Home | History | Annotate | Download | only in ppc
      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_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 
     10 #include "DD.h"
     11 
     12 long double __gcc_qsub(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 (fabs(A) <= 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 (fabs(a) <= 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