1276789Sdim/* This file is distributed under the University of Illinois Open Source 2276789Sdim * License. See LICENSE.TXT for details. 3276789Sdim */ 4276789Sdim 5276789Sdim/* long double __gcc_qadd(long double x, long double y); 6276789Sdim * This file implements the PowerPC 128-bit double-double add operation. 7276789Sdim * This implementation is shamelessly cribbed from Apple's DDRT, circa 1993(!) 8276789Sdim */ 9276789Sdim 10276789Sdim#include "DD.h" 11276789Sdim 12276789Sdimlong double __gcc_qadd(long double x, long double y) 13276789Sdim{ 14276789Sdim static const uint32_t infinityHi = UINT32_C(0x7ff00000); 15276789Sdim 16276789Sdim DD dst = { .ld = x }, src = { .ld = y }; 17276789Sdim 18276789Sdim register double A = dst.s.hi, a = dst.s.lo, 19276789Sdim B = src.s.hi, b = src.s.lo; 20276789Sdim 21276789Sdim /* If both operands are zero: */ 22276789Sdim if ((A == 0.0) && (B == 0.0)) { 23276789Sdim dst.s.hi = A + B; 24276789Sdim dst.s.lo = 0.0; 25276789Sdim return dst.ld; 26276789Sdim } 27276789Sdim 28276789Sdim /* If either operand is NaN or infinity: */ 29276789Sdim const doublebits abits = { .d = A }; 30276789Sdim const doublebits bbits = { .d = B }; 31276789Sdim if ((((uint32_t)(abits.x >> 32) & infinityHi) == infinityHi) || 32276789Sdim (((uint32_t)(bbits.x >> 32) & infinityHi) == infinityHi)) { 33276789Sdim dst.s.hi = A + B; 34276789Sdim dst.s.lo = 0.0; 35276789Sdim return dst.ld; 36276789Sdim } 37276789Sdim 38276789Sdim /* If the computation overflows: */ 39276789Sdim /* This may be playing things a little bit fast and loose, but it will do for a start. */ 40276789Sdim const double testForOverflow = A + (B + (a + b)); 41276789Sdim const doublebits testbits = { .d = testForOverflow }; 42276789Sdim if (((uint32_t)(testbits.x >> 32) & infinityHi) == infinityHi) { 43276789Sdim dst.s.hi = testForOverflow; 44276789Sdim dst.s.lo = 0.0; 45276789Sdim return dst.ld; 46276789Sdim } 47276789Sdim 48276789Sdim double H, h; 49276789Sdim double T, t; 50276789Sdim double W, w; 51276789Sdim double Y; 52276789Sdim 53276789Sdim H = B + (A - (A + B)); 54276789Sdim T = b + (a - (a + b)); 55276789Sdim h = A + (B - (A + B)); 56276789Sdim t = a + (b - (a + b)); 57276789Sdim 58276789Sdim if (local_fabs(A) <= local_fabs(B)) 59276789Sdim w = (a + b) + h; 60276789Sdim else 61276789Sdim w = (a + b) + H; 62276789Sdim 63276789Sdim W = (A + B) + w; 64276789Sdim Y = (A + B) - W; 65276789Sdim Y += w; 66276789Sdim 67276789Sdim if (local_fabs(a) <= local_fabs(b)) 68276789Sdim w = t + Y; 69276789Sdim else 70276789Sdim w = T + Y; 71276789Sdim 72276789Sdim dst.s.hi = Y = W + w; 73276789Sdim dst.s.lo = (W - Y) + w; 74276789Sdim 75276789Sdim return dst.ld; 76276789Sdim} 77