contrib/gcc/config/rs6000/darwin-ldouble.c

   1 /* 128-bit long double support routines for Darwin.
   2    Copyright (C) 1993, 2003, 2004, 2005 Free Software Foundation, Inc.
   3
   4 This file is part of GCC.
   5
   6 GCC is free software; you can redistribute it and/or modify it under
   7 the terms of the GNU General Public License as published by the Free
   8 Software Foundation; either version 2, or (at your option) any later
   9 version.
  10
  11 In addition to the permissions in the GNU General Public License, the
  12 Free Software Foundation gives you unlimited permission to link the
  13 compiled version of this file into combinations with other programs,
  14 and to distribute those combinations without any restriction coming
  15 from the use of this file.  (The General Public License restrictions
  16 do apply in other respects; for example, they cover modification of
  17 the file, and distribution when not linked into a combine
  18 executable.)
  19
  20 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
  21 WARRANTY; without even the implied warranty of MERCHANTABILITY or
  22 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  23 for more details.
  24
  25 You should have received a copy of the GNU General Public License
  26 along with GCC; see the file COPYING.  If not, write to the Free
  27 Software Foundation, 59 Temple Place - Suite 330, Boston, MA
  28 02111-1307, USA.  */
  29
  30 /* Implementations of floating-point long double basic arithmetic
  31    functions called by the IBM C compiler when generating code for
  32    PowerPC platforms.  In particular, the following functions are
  33    implemented: _xlqadd, _xlqsub, _xlqmul, and _xlqdiv.  Double-double
  34    algorithms are based on the paper "Doubled-Precision IEEE Standard
  35    754 Floating-Point Arithmetic" by W. Kahan, February 26, 1987.  An
  36    alternative published reference is "Software for Doubled-Precision
  37    Floating-Point Computations", by Seppo Linnainmaa, ACM TOMS vol 7
  38    no 3, September 1961, pages 272-283.  */
  39
  40 /* Each long double is made up of two IEEE doubles.  The value of the
  41    long double is the sum of the values of the two parts.  The most
  42    significant part is required to be the value of the long double
  43    rounded to the nearest double, as specified by IEEE.  For Inf
  44    values, the least significant part is required to be one of +0.0 or
  45    -0.0.  No other requirements are made; so, for example, 1.0 may be
  46    represented as (1.0, +0.0) or (1.0, -0.0), and the low part of a
  47    NaN is don't-care.
  48
  49    This code currently assumes big-endian.  */
  50
  51 #if !_SOFT_FLOAT && (defined (__MACH__) || defined (__powerpc64__) || defined (_AIX))
  52
  53 #define fabs(x) __builtin_fabs(x)
  54
  55 #define unlikely(x) __builtin_expect ((x), 0)
  56
  57 /* All these routines actually take two long doubles as parameters,
  58    but GCC currently generates poor code when a union is used to turn
  59    a long double into a pair of doubles.  */
  60
  61 extern long double __gcc_qadd (double, double, double, double);
  62 extern long double __gcc_qsub (double, double, double, double);
  63 extern long double __gcc_qmul (double, double, double, double);
  64 extern long double __gcc_qdiv (double, double, double, double);
  65
  66 #if defined __ELF__ && defined IN_LIBGCC2_S
  67 /* Provide definitions of the old symbol names to statisfy apps and
  68    shared libs built against an older libgcc.  To access the _xlq
  69    symbols an explicit version reference is needed, so these won't
  70    satisfy an unadorned reference like _xlqadd.  If dot symbols are
  71    not needed, the assembler will remove the aliases from the symbol
  72    table.  */
  73 __asm__ (".symver __gcc_qadd,_xlqadd@GCC_3.4\n\t"
  74          ".symver __gcc_qsub,_xlqsub@GCC_3.4\n\t"
  75          ".symver __gcc_qmul,_xlqmul@GCC_3.4\n\t"
  76          ".symver __gcc_qdiv,_xlqdiv@GCC_3.4\n\t"
  77          ".symver .__gcc_qadd,._xlqadd@GCC_3.4\n\t"
  78          ".symver .__gcc_qsub,._xlqsub@GCC_3.4\n\t"
  79          ".symver .__gcc_qmul,._xlqmul@GCC_3.4\n\t"
  80          ".symver .__gcc_qdiv,._xlqdiv@GCC_3.4");
  81 #endif
  82
  83 typedef union
  84 {
  85   long double ldval;
  86   double dval[2];
  87 } longDblUnion;
  88
  89 static const double FPKINF = 1.0/0.0;
  90
  91 /* Add two 'long double' values and return the result.  */
  92 long double
  93 __gcc_qadd (double a, double b, double c, double d)
  94 {
  95   longDblUnion z;
  96   double t, tau, u, FPR_zero, FPR_PosInf;
  97
  98   FPR_zero = 0.0;
  99   FPR_PosInf = FPKINF;
 100
 101   if (unlikely (a != a) || unlikely (c != c))
 102     return a + c;  /* NaN result.  */
 103
 104   /* Ordered operands are arranged in order of their magnitudes.  */
 105
 106   /* Switch inputs if |(c,d)| > |(a,b)|. */
 107   if (fabs (c) > fabs (a))
 108     {
 109       t = a;
 110       tau = b;
 111       a = c;
 112       b = d;
 113       c = t;
 114       d = tau;
 115     }
 116
 117   /* b <- second largest magnitude double. */
 118   if (fabs (c) > fabs (b))
 119     {
 120       t = b;
 121       b = c;
 122       c = t;
 123     }
 124
 125   /* Thanks to commutivity, sum is invariant w.r.t. the next
 126      conditional exchange. */
 127   tau = d + c;
 128
 129   /* Order the smallest magnitude doubles.  */
 130   if (fabs (d) > fabs (c))
 131     {
 132       t = c;
 133       c = d;
 134       d = t;
 135     }
 136
 137   t = (tau + b) + a;         /* Sum values in ascending magnitude order.  */
 138
 139   /* Infinite or zero result.  */
 140   if (unlikely (t == FPR_zero) || unlikely (fabs (t) == FPR_PosInf))
 141     return t;
 142
 143   /* Usual case.  */
 144   tau = (((a-t) + b) + c) + d;
 145   u = t + tau;
 146   z.dval[0] = u;               /* Final fixup for long double result.  */
 147   z.dval[1] = (t - u) + tau;
 148   return z.ldval;
 149 }
 150
 151 long double
 152 __gcc_qsub (double a, double b, double c, double d)
 153 {
 154   return __gcc_qadd (a, b, -c, -d);
 155 }
 156
 157 long double
 158 __gcc_qmul (double a, double b, double c, double d)
 159 {
 160   longDblUnion z;
 161   double t, tau, u, v, w, FPR_zero, FPR_PosInf;
 162
 163   FPR_zero = 0.0;
 164   FPR_PosInf = FPKINF;
 165
 166   t = a * c;                    /* Highest order double term.  */
 167
 168   if (unlikely (t != t) || unlikely (t == FPR_zero)
 169       || unlikely (fabs (t) == FPR_PosInf))
 170     return t;
 171
 172   /* Finite nonzero result requires summing of terms of two highest
 173      orders.    */
 174
 175   /* Use fused multiply-add to get low part of a * c.    */
 176   asm ("fmsub %0,%1,%2,%3" : "=f"(tau) : "f"(a), "f"(c), "f"(t));
 177   v = a*d;
 178   w = b*c;
 179   tau += v + w;     /* Add in other second-order terms.  */
 180   u = t + tau;
 181
 182   /* Construct long double result.  */
 183   z.dval[0] = u;
 184   z.dval[1] = (t - u) + tau;
 185   return z.ldval;
 186 }
 187
 188 long double
 189 __gcc_qdiv (double a, double b, double c, double d)
 190 {
 191   longDblUnion z;
 192   double s, sigma, t, tau, u, v, w, FPR_zero, FPR_PosInf;
 193
 194   FPR_zero = 0.0;
 195   FPR_PosInf = FPKINF;
 196
 197   t = a / c;                    /* highest order double term */
 198
 199   if (unlikely (t != t) || unlikely (t == FPR_zero)
 200       || unlikely (fabs (t) == FPR_PosInf))
 201     return t;
 202
 203   /* Finite nonzero result requires corrections to the highest order term.  */
 204
 205   s = c * t;                    /* (s,sigma) = c*t exactly. */
 206   w = -(-b + d * t);    /* Written to get fnmsub for speed, but not
 207                            numerically necessary.  */
 208
 209   /* Use fused multiply-add to get low part of c * t.    */
 210   asm ("fmsub %0,%1,%2,%3" : "=f"(sigma) : "f"(c), "f"(t), "f"(s));
 211   v = a - s;
 212
 213   tau = ((v-sigma)+w)/c;   /* Correction to t. */
 214   u = t + tau;
 215
 216   /* Construct long double result. */
 217   z.dval[0] = u;
 218   z.dval[1] = (t - u) + tau;
 219   return z.ldval;
 220 }
 221
 222 #endif