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_qdiv(long double x, long double y); 6 * This file implements the PowerPC 128-bit double-double division operation. 7 * This implementation is shamelessly cribbed from Apple's DDRT, circa 1993(!) 8 */ 9 10 #include "DD.h" 11 12 long double __gcc_qdiv(long double a, long double b) 13 { 14 static const uint32_t infinityHi = UINT32_C(0x7ff00000); 15 DD dst = { .ld = a }, src = { .ld = b }; 16 17 register double x = dst.s.hi, x1 = dst.s.lo, 18 y = src.s.hi, y1 = src.s.lo; 19 20 double yHi, yLo, qHi, qLo; 21 double yq, tmp, q; 22 23 q = x / y; 24 25 /* Detect special cases */ 26 if (q == 0.0) { 27 dst.s.hi = q; 28 dst.s.lo = 0.0; 29 return dst.ld; 30 } 31 32 const doublebits qBits = { .d = q }; 33 if (((uint32_t)(qBits.x >> 32) & infinityHi) == infinityHi) { 34 dst.s.hi = q; 35 dst.s.lo = 0.0; 36 return dst.ld; 37 } 38 39 yHi = high26bits(y); 40 qHi = high26bits(q); 41 42 yq = y * q; 43 yLo = y - yHi; 44 qLo = q - qHi; 45 46 tmp = LOWORDER(yq, yHi, yLo, qHi, qLo); 47 tmp = (x - yq) - tmp; 48 tmp = ((tmp + x1) - y1 * q) / y; 49 x = q + tmp; 50 51 dst.s.lo = (q - x) + tmp; 52 dst.s.hi = x; 53 54 return dst.ld; 55 } 56
