1 | // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. |
2 | // See https://llvm.org/LICENSE.txt for license information. |
3 | // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception |
4 | |
5 | // long double __gcc_qsub(long double x, long double y); |
6 | // This file implements the PowerPC 128-bit double-double add operation. |
7 | // This implementation is shamelessly cribbed from Apple's DDRT, circa 1993(!) |
8 | |
9 | #include "DD.h" |
10 | |
11 | long double __gcc_qsub(long double x, long double y) { |
12 | static const uint32_t infinityHi = UINT32_C(0x7ff00000); |
13 | |
14 | DD dst = {.ld = x}, src = {.ld = y}; |
15 | |
16 | register double A = dst.s.hi, a = dst.s.lo, B = -src.s.hi, b = -src.s.lo; |
17 | |
18 | // If both operands are zero: |
19 | if ((A == 0.0) && (B == 0.0)) { |
20 | dst.s.hi = A + B; |
21 | dst.s.lo = 0.0; |
22 | return dst.ld; |
23 | } |
24 | |
25 | // If either operand is NaN or infinity: |
26 | const doublebits abits = {.d = A}; |
27 | const doublebits bbits = {.d = B}; |
28 | if ((((uint32_t)(abits.x >> 32) & infinityHi) == infinityHi) || |
29 | (((uint32_t)(bbits.x >> 32) & infinityHi) == infinityHi)) { |
30 | dst.s.hi = A + B; |
31 | dst.s.lo = 0.0; |
32 | return dst.ld; |
33 | } |
34 | |
35 | // If the computation overflows: |
36 | // This may be playing things a little bit fast and loose, but it will do for |
37 | // a start. |
38 | const double testForOverflow = A + (B + (a + b)); |
39 | const doublebits testbits = {.d = testForOverflow}; |
40 | if (((uint32_t)(testbits.x >> 32) & infinityHi) == infinityHi) { |
41 | dst.s.hi = testForOverflow; |
42 | dst.s.lo = 0.0; |
43 | return dst.ld; |
44 | } |
45 | |
46 | double H, h; |
47 | double T, t; |
48 | double W, w; |
49 | double Y; |
50 | |
51 | H = B + (A - (A + B)); |
52 | T = b + (a - (a + b)); |
53 | h = A + (B - (A + B)); |
54 | t = a + (b - (a + b)); |
55 | |
56 | if (local_fabs(x: A) <= local_fabs(x: B)) |
57 | w = (a + b) + h; |
58 | else |
59 | w = (a + b) + H; |
60 | |
61 | W = (A + B) + w; |
62 | Y = (A + B) - W; |
63 | Y += w; |
64 | |
65 | if (local_fabs(x: a) <= local_fabs(x: b)) |
66 | w = t + Y; |
67 | else |
68 | w = T + Y; |
69 | |
70 | dst.s.hi = Y = W + w; |
71 | dst.s.lo = (W - Y) + w; |
72 | |
73 | return dst.ld; |
74 | } |
75 | |