• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
1 /* This file is distributed under the University of Illinois Open Source
2  *  License. See LICENSE.TXT for details.
3  */
4 
5 /* long double __gcc_qadd(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 
10 #include "DD.h"
11 
__gcc_qadd(long double x,long double y)12 long double __gcc_qadd(long double x, long double y)
13 {
14 	static const uint32_t infinityHi = UINT32_C(0x7ff00000);
15 
16 	DD dst = { .ld = x }, src = { .ld = y };
17 
18 	register double A = dst.s.hi, a = dst.s.lo,
19 					B = src.s.hi, b = src.s.lo;
20 
21 	/* If both operands are zero: */
22 	if ((A == 0.0) && (B == 0.0)) {
23 		dst.s.hi = A + B;
24 		dst.s.lo = 0.0;
25 		return dst.ld;
26 	}
27 
28 	/* If either operand is NaN or infinity: */
29 	const doublebits abits = { .d = A };
30 	const doublebits bbits = { .d = B };
31 	if ((((uint32_t)(abits.x >> 32) & infinityHi) == infinityHi) ||
32 		(((uint32_t)(bbits.x >> 32) & infinityHi) == infinityHi)) {
33 		dst.s.hi = A + B;
34 		dst.s.lo = 0.0;
35 		return dst.ld;
36 	}
37 
38 	/* If the computation overflows: */
39 	/* This may be playing things a little bit fast and loose, but it will do for a start. */
40 	const double testForOverflow = A + (B + (a + b));
41 	const doublebits testbits = { .d = testForOverflow };
42 	if (((uint32_t)(testbits.x >> 32) & infinityHi) == infinityHi) {
43 		dst.s.hi = testForOverflow;
44 		dst.s.lo = 0.0;
45 		return dst.ld;
46 	}
47 
48 	double H, h;
49 	double T, t;
50 	double W, w;
51 	double Y;
52 
53 	H = B + (A - (A + B));
54 	T = b + (a - (a + b));
55 	h = A + (B - (A + B));
56 	t = a + (b - (a + b));
57 
58 	if (local_fabs(A) <= local_fabs(B))
59 		w = (a + b) + h;
60 	else
61 		w = (a + b) + H;
62 
63 	W = (A + B) + w;
64 	Y = (A + B) - W;
65 	Y += w;
66 
67 	if (local_fabs(a) <= local_fabs(b))
68 		w = t + Y;
69 	else
70 		w = T + Y;
71 
72 	dst.s.hi = Y = W + w;
73 	dst.s.lo = (W - Y) + w;
74 
75 	return dst.ld;
76 }
77