lib/ppc/gcc_qadd.c

   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_qadd(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_qadd(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 }