gcc/config/rs6000/darwin-ldouble.c

   1 /* 128-bit long double support routines for Darwin.
   2    Copyright (C) 1993, 2003, 2004, 2005, 2006, 2007
   3    Free Software Foundation, Inc.
   4
   5 This file is part of GCC.
   6
   7 GCC is free software; you can redistribute it and/or modify it under
   8 the terms of the GNU General Public License as published by the Free
   9 Software Foundation; either version 2, or (at your option) any later
  10 version.
  11
  12 In addition to the permissions in the GNU General Public License, the
  13 Free Software Foundation gives you unlimited permission to link the
  14 compiled version of this file into combinations with other programs,
  15 and to distribute those combinations without any restriction coming
  16 from the use of this file.  (The General Public License restrictions
  17 do apply in other respects; for example, they cover modification of
  18 the file, and distribution when not linked into a combine
  19 executable.)
  20
  21 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
  22 WARRANTY; without even the implied warranty of MERCHANTABILITY or
  23 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  24 for more details.
  25
  26 You should have received a copy of the GNU General Public License
  27 along with GCC; see the file COPYING.  If not, write to the Free
  28 Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA
  29 02110-1301, USA.  */
  30
  31 /* Implementations of floating-point long double basic arithmetic
  32    functions called by the IBM C compiler when generating code for
  33    PowerPC platforms.  In particular, the following functions are
  34    implemented: __gcc_qadd, __gcc_qsub, __gcc_qmul, and __gcc_qdiv.
  35    Double-double algorithms are based on the paper "Doubled-Precision
  36    IEEE Standard 754 Floating-Point Arithmetic" by W. Kahan, February 26,
  37    1987.  An alternative published reference is "Software for
  38    Doubled-Precision Floating-Point Computations", by Seppo Linnainmaa,
  39    ACM TOMS vol 7 no 3, September 1981, pages 272-283.  */
  40
  41 /* Each long double is made up of two IEEE doubles.  The value of the
  42    long double is the sum of the values of the two parts.  The most
  43    significant part is required to be the value of the long double
  44    rounded to the nearest double, as specified by IEEE.  For Inf
  45    values, the least significant part is required to be one of +0.0 or
  46    -0.0.  No other requirements are made; so, for example, 1.0 may be
  47    represented as (1.0, +0.0) or (1.0, -0.0), and the low part of a
  48    NaN is don't-care.
  49
  50    This code currently assumes big-endian.  */
  51
  52 #if (!defined (__LITTLE_ENDIAN__) \
  53      && (defined (__MACH__) || defined (__powerpc__) || defined (_AIX)))
  54
  55 #define fabs(x) __builtin_fabs(x)
  56 #define isless(x, y) __builtin_isless (x, y)
  57 #define inf() __builtin_inf()
  58
  59 #define unlikely(x) __builtin_expect ((x), 0)
  60
  61 #define nonfinite(a) unlikely (! isless (fabs (a), inf ()))
  62
  63 /* Define ALIASNAME as a strong alias for NAME.  */
  64 # define strong_alias(name, aliasname) _strong_alias(name, aliasname)
  65 # define _strong_alias(name, aliasname) \
  66   extern __typeof (name) aliasname __attribute__ ((alias (#name)));
  67
  68 /* All these routines actually take two long doubles as parameters,
  69    but GCC currently generates poor code when a union is used to turn
  70    a long double into a pair of doubles.  */
  71
  72 long double __gcc_qadd (double, double, double, double);
  73 long double __gcc_qsub (double, double, double, double);
  74 long double __gcc_qmul (double, double, double, double);
  75 long double __gcc_qdiv (double, double, double, double);
  76
  77 #if defined __ELF__ && defined SHARED \
  78     && (defined __powerpc64__ || !(defined __linux__ || defined __gnu_hurd__))
  79 /* Provide definitions of the old symbol names to satisfy apps and
  80    shared libs built against an older libgcc.  To access the _xlq
  81    symbols an explicit version reference is needed, so these won't
  82    satisfy an unadorned reference like _xlqadd.  If dot symbols are
  83    not needed, the assembler will remove the aliases from the symbol
  84    table.  */
  85 __asm__ (".symver __gcc_qadd,_xlqadd@GCC_3.4\n\t"
  86          ".symver __gcc_qsub,_xlqsub@GCC_3.4\n\t"
  87          ".symver __gcc_qmul,_xlqmul@GCC_3.4\n\t"
  88          ".symver __gcc_qdiv,_xlqdiv@GCC_3.4\n\t"
  89          ".symver .__gcc_qadd,._xlqadd@GCC_3.4\n\t"
  90          ".symver .__gcc_qsub,._xlqsub@GCC_3.4\n\t"
  91          ".symver .__gcc_qmul,._xlqmul@GCC_3.4\n\t"
  92          ".symver .__gcc_qdiv,._xlqdiv@GCC_3.4");
  93 #endif
  94
  95 typedef union
  96 {
  97   long double ldval;
  98   double dval[2];
  99 } longDblUnion;
 100
 101 /* Add two 'long double' values and return the result.  */
 102 long double
 103 __gcc_qadd (double a, double aa, double c, double cc)
 104 {
 105   longDblUnion x;
 106   double z, q, zz, xh;
 107
 108   z = a + c;
 109
 110   if (nonfinite (z))
 111     {
 112       z = cc + aa + c + a;
 113       if (nonfinite (z))
 114         return z;
 115       x.dval[0] = z;  /* Will always be DBL_MAX.  */
 116       zz = aa + cc;
 117       if (fabs(a) > fabs(c))
 118         x.dval[1] = a - z + c + zz;
 119       else
 120         x.dval[1] = c - z + a + zz;
 121     }
 122   else
 123     {
 124       q = a - z;
 125       zz = q + c + (a - (q + z)) + aa + cc;
 126
 127       /* Keep -0 result.  */
 128       if (zz == 0.0)
 129         return z;
 130
 131       xh = z + zz;
 132       if (nonfinite (xh))
 133         return xh;
 134
 135       x.dval[0] = xh;
 136       x.dval[1] = z - xh + zz;
 137     }
 138   return x.ldval;
 139 }
 140
 141 long double
 142 __gcc_qsub (double a, double b, double c, double d)
 143 {
 144   return __gcc_qadd (a, b, -c, -d);
 145 }
 146
 147 #ifdef __NO_FPRS__
 148 static double fmsub (double, double, double);
 149 #endif
 150
 151 long double
 152 __gcc_qmul (double a, double b, double c, double d)
 153 {
 154   longDblUnion z;
 155   double t, tau, u, v, w;
 156
 157   t = a * c;                    /* Highest order double term.  */
 158
 159   if (unlikely (t == 0)         /* Preserve -0.  */
 160       || nonfinite (t))
 161     return t;
 162
 163   /* Sum terms of two highest orders. */
 164
 165   /* Use fused multiply-add to get low part of a * c.  */
 166 #ifndef __NO_FPRS__
 167   asm ("fmsub %0,%1,%2,%3" : "=f"(tau) : "f"(a), "f"(c), "f"(t));
 168 #else
 169   tau = fmsub (a, c, t);
 170 #endif
 171   v = a*d;
 172   w = b*c;
 173   tau += v + w;     /* Add in other second-order terms.  */
 174   u = t + tau;
 175
 176   /* Construct long double result.  */
 177   if (nonfinite (u))
 178     return u;
 179   z.dval[0] = u;
 180   z.dval[1] = (t - u) + tau;
 181   return z.ldval;
 182 }
 183
 184 long double
 185 __gcc_qdiv (double a, double b, double c, double d)
 186 {
 187   longDblUnion z;
 188   double s, sigma, t, tau, u, v, w;
 189
 190   t = a / c;                    /* highest order double term */
 191
 192   if (unlikely (t == 0)         /* Preserve -0.  */
 193       || nonfinite (t))
 194     return t;
 195
 196   /* Finite nonzero result requires corrections to the highest order term.  */
 197
 198   s = c * t;                    /* (s,sigma) = c*t exactly.  */
 199   w = -(-b + d * t);    /* Written to get fnmsub for speed, but not
 200                            numerically necessary.  */
 201
 202   /* Use fused multiply-add to get low part of c * t.    */
 203 #ifndef __NO_FPRS__
 204   asm ("fmsub %0,%1,%2,%3" : "=f"(sigma) : "f"(c), "f"(t), "f"(s));
 205 #else
 206   sigma = fmsub (c, t, s);
 207 #endif
 208   v = a - s;
 209
 210   tau = ((v-sigma)+w)/c;   /* Correction to t.  */
 211   u = t + tau;
 212
 213   /* Construct long double result.  */
 214   if (nonfinite (u))
 215     return u;
 216   z.dval[0] = u;
 217   z.dval[1] = (t - u) + tau;
 218   return z.ldval;
 219 }
 220
 221 #if defined (_SOFT_DOUBLE) && defined (__LONG_DOUBLE_128__)
 222
 223 long double __gcc_qneg (double, double);
 224 int __gcc_qeq (double, double, double, double);
 225 int __gcc_qne (double, double, double, double);
 226 int __gcc_qge (double, double, double, double);
 227 int __gcc_qle (double, double, double, double);
 228 long double __gcc_stoq (float);
 229 long double __gcc_dtoq (double);
 230 float __gcc_qtos (double, double);
 231 double __gcc_qtod (double, double);
 232 int __gcc_qtoi (double, double);
 233 unsigned int __gcc_qtou (double, double);
 234 long double __gcc_itoq (int);
 235 long double __gcc_utoq (unsigned int);
 236
 237 extern int __eqdf2 (double, double);
 238 extern int __ledf2 (double, double);
 239 extern int __gedf2 (double, double);
 240
 241 /* Negate 'long double' value and return the result.    */
 242 long double
 243 __gcc_qneg (double a, double aa)
 244 {
 245   longDblUnion x;
 246
 247   x.dval[0] = -a;
 248   x.dval[1] = -aa;
 249   return x.ldval;
 250 }
 251
 252 /* Compare two 'long double' values for equality.  */
 253 int
 254 __gcc_qeq (double a, double aa, double c, double cc)
 255 {
 256   if (__eqdf2 (a, c) == 0)
 257     return __eqdf2 (aa, cc);
 258   return 1;
 259 }
 260
 261 strong_alias (__gcc_qeq, __gcc_qne);
 262
 263 /* Compare two 'long double' values for less than or equal.  */
 264 int
 265 __gcc_qle (double a, double aa, double c, double cc)
 266 {
 267   if (__eqdf2 (a, c) == 0)
 268     return __ledf2 (aa, cc);
 269   return __ledf2 (a, c);
 270 }
 271
 272 strong_alias (__gcc_qle, __gcc_qlt);
 273
 274 /* Compare two 'long double' values for greater than or equal.  */
 275 int
 276 __gcc_qge (double a, double aa, double c, double cc)
 277 {
 278   if (__eqdf2 (a, c) == 0)
 279     return __gedf2 (aa, cc);
 280   return __gedf2 (a, c);
 281 }
 282
 283 strong_alias (__gcc_qge, __gcc_qgt);
 284
 285 /* Convert single to long double.  */
 286 long double
 287 __gcc_stoq (float a)
 288 {
 289   longDblUnion x;
 290
 291   x.dval[0] = (double) a;
 292   x.dval[1] = 0.0;
 293
 294   return x.ldval;
 295 }
 296
 297 /* Convert double to long double.  */
 298 long double
 299 __gcc_dtoq (double a)
 300 {
 301   longDblUnion x;
 302
 303   x.dval[0] = a;
 304   x.dval[1] = 0.0;
 305
 306   return x.ldval;
 307 }
 308
 309 /* Convert long double to single.  */
 310 float
 311 __gcc_qtos (double a, double aa __attribute__ ((__unused__)))
 312 {
 313   return (float) a;
 314 }
 315
 316 /* Convert long double to double.  */
 317 double
 318 __gcc_qtod (double a, double aa __attribute__ ((__unused__)))
 319 {
 320   return a;
 321 }
 322
 323 /* Convert long double to int.  */
 324 int
 325 __gcc_qtoi (double a, double aa)
 326 {
 327   double z = a + aa;
 328   return (int) z;
 329 }
 330
 331 /* Convert long double to unsigned int.  */
 332 unsigned int
 333 __gcc_qtou (double a, double aa)
 334 {
 335   double z = a + aa;
 336   return (unsigned int) z;
 337 }
 338
 339 /* Convert int to long double.  */
 340 long double
 341 __gcc_itoq (int a)
 342 {
 343   return __gcc_dtoq ((double) a);
 344 }
 345
 346 /* Convert unsigned int to long double.  */
 347 long double
 348 __gcc_utoq (unsigned int a)
 349 {
 350   return __gcc_dtoq ((double) a);
 351 }
 352
 353 #endif
 354
 355 #ifdef __NO_FPRS__
 356
 357 int __gcc_qunord (double, double, double, double);
 358
 359 extern int __eqdf2 (double, double);
 360 extern int __unorddf2 (double, double);
 361
 362 /* Compare two 'long double' values for unordered.  */
 363 int
 364 __gcc_qunord (double a, double aa, double c, double cc)
 365 {
 366   if (__eqdf2 (a, c) == 0)
 367     return __unorddf2 (aa, cc);
 368   return __unorddf2 (a, c);
 369 }
 370
 371 #include "config/soft-fp/soft-fp.h"
 372 #include "config/soft-fp/double.h"
 373 #include "config/soft-fp/quad.h"
 374
 375 /* Compute floating point multiply-subtract with higher (quad) precision.  */
 376 static double
 377 fmsub (double a, double b, double c)
 378 {
 379     FP_DECL_EX;
 380     FP_DECL_D(A);
 381     FP_DECL_D(B);
 382     FP_DECL_D(C);
 383     FP_DECL_Q(X);
 384     FP_DECL_Q(Y);
 385     FP_DECL_Q(Z);
 386     FP_DECL_Q(U);
 387     FP_DECL_Q(V);
 388     FP_DECL_D(R);
 389     double r;
 390     long double u, v, x, y, z;
 391
 392     FP_INIT_ROUNDMODE;
 393     FP_UNPACK_RAW_D (A, a);
 394     FP_UNPACK_RAW_D (B, b);
 395     FP_UNPACK_RAW_D (C, c);
 396
 397     /* Extend double to quad.  */
 398 #if (2 * _FP_W_TYPE_SIZE) < _FP_FRACBITS_Q
 399     FP_EXTEND(Q,D,4,2,X,A);
 400     FP_EXTEND(Q,D,4,2,Y,B);
 401     FP_EXTEND(Q,D,4,2,Z,C);
 402 #else
 403     FP_EXTEND(Q,D,2,1,X,A);
 404     FP_EXTEND(Q,D,2,1,Y,B);
 405     FP_EXTEND(Q,D,2,1,Z,C);
 406 #endif
 407     FP_PACK_RAW_Q(x,X);
 408     FP_PACK_RAW_Q(y,Y);
 409     FP_PACK_RAW_Q(z,Z);
 410     FP_HANDLE_EXCEPTIONS;
 411
 412     /* Multiply.  */
 413     FP_INIT_ROUNDMODE;
 414     FP_UNPACK_Q(X,x);
 415     FP_UNPACK_Q(Y,y);
 416     FP_MUL_Q(U,X,Y);
 417     FP_PACK_Q(u,U);
 418     FP_HANDLE_EXCEPTIONS;
 419
 420     /* Subtract.  */
 421     FP_INIT_ROUNDMODE;
 422     FP_UNPACK_SEMIRAW_Q(U,u);
 423     FP_UNPACK_SEMIRAW_Q(Z,z);
 424     FP_SUB_Q(V,U,Z);
 425
 426     /* Truncate quad to double.  */
 427 #if (2 * _FP_W_TYPE_SIZE) < _FP_FRACBITS_Q
 428     V_f[3] &= 0x0007ffff;
 429     FP_TRUNC(D,Q,2,4,R,V);
 430 #else
 431     V_f1 &= 0x0007ffffffffffffL;
 432     FP_TRUNC(D,Q,1,2,R,V);
 433 #endif
 434     FP_PACK_SEMIRAW_D(r,R);
 435     FP_HANDLE_EXCEPTIONS;
 436
 437     return r;
 438 }
 439
 440 #endif
 441
 442 #endif