1214152Sed/* This file is distributed under the University of Illinois Open Source 2214152Sed * License. See LICENSE.TXT for details. 3214152Sed */ 4214152Sed 5214152Sed/* long double __gcc_qsub(long double x, long double y); 6214152Sed * This file implements the PowerPC 128-bit double-double add operation. 7214152Sed * This implementation is shamelessly cribbed from Apple's DDRT, circa 1993(!) 8214152Sed */ 9214152Sed 10214152Sed#include "DD.h" 11214152Sed 12214152Sedlong double __gcc_qsub(long double x, long double y) 13214152Sed{ 14214152Sed static const uint32_t infinityHi = UINT32_C(0x7ff00000); 15214152Sed 16214152Sed DD dst = { .ld = x }, src = { .ld = y }; 17214152Sed 18214152Sed register double A = dst.s.hi, a = dst.s.lo, 19214152Sed B = -src.s.hi, b = -src.s.lo; 20214152Sed 21214152Sed /* If both operands are zero: */ 22214152Sed if ((A == 0.0) && (B == 0.0)) { 23214152Sed dst.s.hi = A + B; 24214152Sed dst.s.lo = 0.0; 25214152Sed return dst.ld; 26214152Sed } 27214152Sed 28214152Sed /* If either operand is NaN or infinity: */ 29214152Sed const doublebits abits = { .d = A }; 30214152Sed const doublebits bbits = { .d = B }; 31214152Sed if ((((uint32_t)(abits.x >> 32) & infinityHi) == infinityHi) || 32214152Sed (((uint32_t)(bbits.x >> 32) & infinityHi) == infinityHi)) { 33214152Sed dst.s.hi = A + B; 34214152Sed dst.s.lo = 0.0; 35214152Sed return dst.ld; 36214152Sed } 37214152Sed 38214152Sed /* If the computation overflows: */ 39214152Sed /* This may be playing things a little bit fast and loose, but it will do for a start. */ 40214152Sed const double testForOverflow = A + (B + (a + b)); 41214152Sed const doublebits testbits = { .d = testForOverflow }; 42214152Sed if (((uint32_t)(testbits.x >> 32) & infinityHi) == infinityHi) { 43214152Sed dst.s.hi = testForOverflow; 44214152Sed dst.s.lo = 0.0; 45214152Sed return dst.ld; 46214152Sed } 47214152Sed 48214152Sed double H, h; 49214152Sed double T, t; 50214152Sed double W, w; 51214152Sed double Y; 52214152Sed 53214152Sed H = B + (A - (A + B)); 54214152Sed T = b + (a - (a + b)); 55214152Sed h = A + (B - (A + B)); 56214152Sed t = a + (b - (a + b)); 57214152Sed 58214152Sed if (fabs(A) <= fabs(B)) 59214152Sed w = (a + b) + h; 60214152Sed else 61214152Sed w = (a + b) + H; 62214152Sed 63214152Sed W = (A + B) + w; 64214152Sed Y = (A + B) - W; 65214152Sed Y += w; 66214152Sed 67214152Sed if (fabs(a) <= fabs(b)) 68214152Sed w = t + Y; 69214152Sed else 70214152Sed w = T + Y; 71214152Sed 72214152Sed dst.s.hi = Y = W + w; 73214152Sed dst.s.lo = (W - Y) + w; 74214152Sed 75214152Sed return dst.ld; 76214152Sed} 77