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