compiler-rt/lib/builtins/ppc/gcc_qsub.c

   1 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
   2 // See https://llvm.org/LICENSE.txt for license information.
   3 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
   4
   5 // long double __gcc_qsub(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 #include "DD.h"
  10
  11 long double __gcc_qsub(long double x, long double y) {
  12   static const uint32_t infinityHi = UINT32_C(0x7ff00000);
  13
  14   DD dst = {.ld = x}, src = {.ld = y};
  15
  16   register double A = dst.s.hi, a = dst.s.lo, B = -src.s.hi, b = -src.s.lo;
  17
  18   // If both operands are zero:
  19   if ((A == 0.0) && (B == 0.0)) {
  20     dst.s.hi = A + B;
  21     dst.s.lo = 0.0;
  22     return dst.ld;
  23   }
  24
  25   // If either operand is NaN or infinity:
  26   const doublebits abits = {.d = A};
  27   const doublebits bbits = {.d = B};
  28   if ((((uint32_t)(abits.x >> 32) & infinityHi) == infinityHi) ||
  29       (((uint32_t)(bbits.x >> 32) & infinityHi) == infinityHi)) {
  30     dst.s.hi = A + B;
  31     dst.s.lo = 0.0;
  32     return dst.ld;
  33   }
  34
  35   // If the computation overflows:
  36   // This may be playing things a little bit fast and loose, but it will do for
  37   // a start.
  38   const double testForOverflow = A + (B + (a + b));
  39   const doublebits testbits = {.d = testForOverflow};
  40   if (((uint32_t)(testbits.x >> 32) & infinityHi) == infinityHi) {
  41     dst.s.hi = testForOverflow;
  42     dst.s.lo = 0.0;
  43     return dst.ld;
  44   }
  45
  46   double H, h;
  47   double T, t;
  48   double W, w;
  49   double Y;
  50
  51   H = B + (A - (A + B));
  52   T = b + (a - (a + b));
  53   h = A + (B - (A + B));
  54   t = a + (b - (a + b));
  55
  56   if (local_fabs(A) <= local_fabs(B))
  57     w = (a + b) + h;
  58   else
  59     w = (a + b) + H;
  60
  61   W = (A + B) + w;
  62   Y = (A + B) - W;
  63   Y += w;
  64
  65   if (local_fabs(a) <= local_fabs(b))
  66     w = t + Y;
  67   else
  68     w = T + Y;
  69
  70   dst.s.hi = Y = W + w;
  71   dst.s.lo = (W - Y) + w;
  72
  73   return dst.ld;
  74 }