18bf1e094893cb24796137b47ee0d46d18d299996Edward O'Callaghan/* This file is distributed under the University of Illinois Open Source 28bf1e094893cb24796137b47ee0d46d18d299996Edward O'Callaghan * License. See LICENSE.TXT for details. 38bf1e094893cb24796137b47ee0d46d18d299996Edward O'Callaghan */ 4b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar 58bf1e094893cb24796137b47ee0d46d18d299996Edward O'Callaghan/* long double __gcc_qadd(long double x, long double y); 68bf1e094893cb24796137b47ee0d46d18d299996Edward O'Callaghan * This file implements the PowerPC 128-bit double-double add operation. 78bf1e094893cb24796137b47ee0d46d18d299996Edward O'Callaghan * This implementation is shamelessly cribbed from Apple's DDRT, circa 1993(!) 88bf1e094893cb24796137b47ee0d46d18d299996Edward O'Callaghan */ 9b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar 10b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar#include "DD.h" 11b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar 12b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbarlong double __gcc_qadd(long double x, long double y) 13b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar{ 14b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar static const uint32_t infinityHi = UINT32_C(0x7ff00000); 15b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar 16b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar DD dst = { .ld = x }, src = { .ld = y }; 17b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar 188bf1e094893cb24796137b47ee0d46d18d299996Edward O'Callaghan register double A = dst.s.hi, a = dst.s.lo, 198bf1e094893cb24796137b47ee0d46d18d299996Edward O'Callaghan B = src.s.hi, b = src.s.lo; 20b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar 218bf1e094893cb24796137b47ee0d46d18d299996Edward O'Callaghan /* If both operands are zero: */ 22b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar if ((A == 0.0) && (B == 0.0)) { 238bf1e094893cb24796137b47ee0d46d18d299996Edward O'Callaghan dst.s.hi = A + B; 248bf1e094893cb24796137b47ee0d46d18d299996Edward O'Callaghan dst.s.lo = 0.0; 25b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar return dst.ld; 26b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar } 27b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar 288bf1e094893cb24796137b47ee0d46d18d299996Edward O'Callaghan /* If either operand is NaN or infinity: */ 29b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar const doublebits abits = { .d = A }; 30b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar const doublebits bbits = { .d = B }; 31b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar if ((((uint32_t)(abits.x >> 32) & infinityHi) == infinityHi) || 32b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar (((uint32_t)(bbits.x >> 32) & infinityHi) == infinityHi)) { 338bf1e094893cb24796137b47ee0d46d18d299996Edward O'Callaghan dst.s.hi = A + B; 348bf1e094893cb24796137b47ee0d46d18d299996Edward O'Callaghan dst.s.lo = 0.0; 35b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar return dst.ld; 36b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar } 37b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar 388bf1e094893cb24796137b47ee0d46d18d299996Edward O'Callaghan /* If the computation overflows: */ 398bf1e094893cb24796137b47ee0d46d18d299996Edward O'Callaghan /* This may be playing things a little bit fast and loose, but it will do for a start. */ 40b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar const double testForOverflow = A + (B + (a + b)); 41b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar const doublebits testbits = { .d = testForOverflow }; 42b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar if (((uint32_t)(testbits.x >> 32) & infinityHi) == infinityHi) { 438bf1e094893cb24796137b47ee0d46d18d299996Edward O'Callaghan dst.s.hi = testForOverflow; 448bf1e094893cb24796137b47ee0d46d18d299996Edward O'Callaghan dst.s.lo = 0.0; 45b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar return dst.ld; 46b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar } 47b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar 48b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar double H, h; 49b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar double T, t; 50b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar double W, w; 51b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar double Y; 52b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar 53b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar H = B + (A - (A + B)); 54b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar T = b + (a - (a + b)); 55b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar h = A + (B - (A + B)); 56b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar t = a + (b - (a + b)); 57b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar 58b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar if (fabs(A) <= fabs(B)) 59b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar w = (a + b) + h; 60b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar else 61b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar w = (a + b) + H; 62b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar 63b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar W = (A + B) + w; 64b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar Y = (A + B) - W; 65b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar Y += w; 66b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar 67b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar if (fabs(a) <= fabs(b)) 68b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar w = t + Y; 69b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar else 70b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar w = T + Y; 71b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar 728bf1e094893cb24796137b47ee0d46d18d299996Edward O'Callaghan dst.s.hi = Y = W + w; 738bf1e094893cb24796137b47ee0d46d18d299996Edward O'Callaghan dst.s.lo = (W - Y) + w; 74b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar 75b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar return dst.ld; 76b3a6901e66f55b35aa9e01bcb24134e6a65ea004Daniel Dunbar} 77