post/lib_powerpc/fpu/darwin-ldouble.c

   1 /*
   2  * Borrowed from GCC 4.2.2 (which still was GPL v2+)
   3  */
   4 /* 128-bit long double support routines for Darwin.
   5    Copyright (C) 1993, 2003, 2004, 2005, 2006, 2007
   6    Free Software Foundation, Inc.
   7
   8 This file is part of GCC.
   9
  10  * SPDX-License-Identifier:     GPL-2.0+
  11  */
  12
  13 /*
  14  * Implementations of floating-point long double basic arithmetic
  15  * functions called by the IBM C compiler when generating code for
  16  * PowerPC platforms.  In particular, the following functions are
  17  * implemented: __gcc_qadd, __gcc_qsub, __gcc_qmul, and __gcc_qdiv.
  18  * Double-double algorithms are based on the paper "Doubled-Precision
  19  * IEEE Standard 754 Floating-Point Arithmetic" by W. Kahan, February 26,
  20  * 1987.  An alternative published reference is "Software for
  21  * Doubled-Precision Floating-Point Computations", by Seppo Linnainmaa,
  22  * ACM TOMS vol 7 no 3, September 1981, pages 272-283.
  23  */
  24
  25 /*
  26  * Each long double is made up of two IEEE doubles.  The value of the
  27  * long double is the sum of the values of the two parts.  The most
  28  * significant part is required to be the value of the long double
  29  * rounded to the nearest double, as specified by IEEE.  For Inf
  30  * values, the least significant part is required to be one of +0.0 or
  31  * -0.0.  No other requirements are made; so, for example, 1.0 may be
  32  * represented as (1.0, +0.0) or (1.0, -0.0), and the low part of a
  33  * NaN is don't-care.
  34  *
  35  * This code currently assumes big-endian.
  36  */
  37
  38 #define fabs(x) __builtin_fabs(x)
  39 #define isless(x, y) __builtin_isless(x, y)
  40 #define inf() __builtin_inf()
  41 #define unlikely(x) __builtin_expect((x), 0)
  42 #define nonfinite(a) unlikely(!isless(fabs(a), inf()))
  43
  44 typedef union {
  45         long double ldval;
  46         double dval[2];
  47 } longDblUnion;
  48
  49 /* Add two 'long double' values and return the result.  */
  50 long double __gcc_qadd(double a, double aa, double c, double cc)
  51 {
  52         longDblUnion x;
  53         double z, q, zz, xh;
  54
  55         z = a + c;
  56
  57         if (nonfinite(z)) {
  58                 z = cc + aa + c + a;
  59                 if (nonfinite(z))
  60                         return z;
  61                 x.dval[0] = z;  /* Will always be DBL_MAX.  */
  62                 zz = aa + cc;
  63                 if (fabs(a) > fabs(c))
  64                         x.dval[1] = a - z + c + zz;
  65                 else
  66                         x.dval[1] = c - z + a + zz;
  67         } else {
  68                 q = a - z;
  69                 zz = q + c + (a - (q + z)) + aa + cc;
  70
  71                 /* Keep -0 result.  */
  72                 if (zz == 0.0)
  73                         return z;
  74
  75                 xh = z + zz;
  76                 if (nonfinite(xh))
  77                         return xh;
  78
  79                 x.dval[0] = xh;
  80                 x.dval[1] = z - xh + zz;
  81         }
  82         return x.ldval;
  83 }
  84
  85 long double __gcc_qsub(double a, double b, double c, double d)
  86 {
  87         return __gcc_qadd(a, b, -c, -d);
  88 }
  89
  90 long double __gcc_qmul(double a, double b, double c, double d)
  91 {
  92         longDblUnion z;
  93         double t, tau, u, v, w;
  94
  95         t = a * c;              /* Highest order double term.  */
  96
  97         if (unlikely(t == 0)    /* Preserve -0.  */
  98             || nonfinite(t))
  99                 return t;
 100
 101         /* Sum terms of two highest orders. */
 102
 103         /* Use fused multiply-add to get low part of a * c.  */
 104 #ifndef __NO_FPRS__
 105         asm("fmsub %0,%1,%2,%3" : "=f"(tau) : "f"(a), "f"(c), "f"(t));
 106 #else
 107         tau = fmsub(a, c, t);
 108 #endif
 109         v = a * d;
 110         w = b * c;
 111         tau += v + w;           /* Add in other second-order terms.  */
 112         u = t + tau;
 113
 114         /* Construct long double result.  */
 115         if (nonfinite(u))
 116                 return u;
 117         z.dval[0] = u;
 118         z.dval[1] = (t - u) + tau;
 119         return z.ldval;
 120 }