contrib/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 (__NO_FPRS__) || defined (_SOFT_FLOAT)) \
  53      && !defined (__LITTLE_ENDIAN__) \
  54      && (defined (__MACH__) || defined (__powerpc__) || defined (_AIX)))
  55
  56 #define fabs(x) __builtin_fabs(x)
  57 #define isless(x, y) __builtin_isless (x, y)
  58 #define inf() __builtin_inf()
  59
  60 #define unlikely(x) __builtin_expect ((x), 0)
  61
  62 #define nonfinite(a) unlikely (! isless (fabs (a), inf ()))
  63
  64 /* Define ALIASNAME as a strong alias for NAME.  */
  65 # define strong_alias(name, aliasname) _strong_alias(name, aliasname)
  66 # define _strong_alias(name, aliasname) \
  67   extern __typeof (name) aliasname __attribute__ ((alias (#name)));
  68
  69 /* All these routines actually take two long doubles as parameters,
  70    but GCC currently generates poor code when a union is used to turn
  71    a long double into a pair of doubles.  */
  72
  73 long double __gcc_qadd (double, double, double, double);
  74 long double __gcc_qsub (double, double, double, double);
  75 long double __gcc_qmul (double, double, double, double);
  76 long double __gcc_qdiv (double, double, double, double);
  77
  78 #if defined __ELF__ && defined SHARED \
  79     && (defined __powerpc64__ || !(defined __linux__ || defined __gnu_hurd__))
  80 /* Provide definitions of the old symbol names to satisfy apps and
  81    shared libs built against an older libgcc.  To access the _xlq
  82    symbols an explicit version reference is needed, so these won't
  83    satisfy an unadorned reference like _xlqadd.  If dot symbols are
  84    not needed, the assembler will remove the aliases from the symbol
  85    table.  */
  86 __asm__ (".symver __gcc_qadd,_xlqadd@GCC_3.4\n\t"
  87          ".symver __gcc_qsub,_xlqsub@GCC_3.4\n\t"
  88          ".symver __gcc_qmul,_xlqmul@GCC_3.4\n\t"
  89          ".symver __gcc_qdiv,_xlqdiv@GCC_3.4\n\t"
  90          ".symver .__gcc_qadd,._xlqadd@GCC_3.4\n\t"
  91          ".symver .__gcc_qsub,._xlqsub@GCC_3.4\n\t"
  92          ".symver .__gcc_qmul,._xlqmul@GCC_3.4\n\t"
  93          ".symver .__gcc_qdiv,._xlqdiv@GCC_3.4");
  94 #endif
  95
  96 typedef union
  97 {
  98   long double ldval;
  99   double dval[2];
 100 } longDblUnion;
 101
 102 /* Add two 'long double' values and return the result.  */
 103 long double
 104 __gcc_qadd (double a, double aa, double c, double cc)
 105 {
 106   longDblUnion x;
 107   double z, q, zz, xh;
 108
 109   z = a + c;
 110
 111   if (nonfinite (z))
 112     {
 113       z = cc + aa + c + a;
 114       if (nonfinite (z))
 115         return z;
 116       x.dval[0] = z;  /* Will always be DBL_MAX.  */
 117       zz = aa + cc;
 118       if (fabs(a) > fabs(c))
 119         x.dval[1] = a - z + c + zz;
 120       else
 121         x.dval[1] = c - z + a + zz;
 122     }
 123   else
 124     {
 125       q = a - z;
 126       zz = q + c + (a - (q + z)) + aa + cc;
 127
 128       /* Keep -0 result.  */
 129       if (zz == 0.0)
 130         return z;
 131
 132       xh = z + zz;
 133       if (nonfinite (xh))
 134         return xh;
 135
 136       x.dval[0] = xh;
 137       x.dval[1] = z - xh + zz;
 138     }
 139   return x.ldval;
 140 }
 141
 142 long double
 143 __gcc_qsub (double a, double b, double c, double d)
 144 {
 145   return __gcc_qadd (a, b, -c, -d);
 146 }
 147
 148 #ifdef _SOFT_FLOAT
 149 static double fmsub (double, double, double);
 150 #endif
 151
 152 long double
 153 __gcc_qmul (double a, double b, double c, double d)
 154 {
 155   longDblUnion z;
 156   double t, tau, u, v, w;
 157
 158   t = a * c;                    /* Highest order double term.  */
 159
 160   if (unlikely (t == 0)         /* Preserve -0.  */
 161       || nonfinite (t))
 162     return t;
 163
 164   /* Sum terms of two highest orders. */
 165
 166   /* Use fused multiply-add to get low part of a * c.  */
 167 #ifndef _SOFT_FLOAT
 168   asm ("fmsub %0,%1,%2,%3" : "=f"(tau) : "f"(a), "f"(c), "f"(t));
 169 #else
 170   tau = fmsub (a, c, t);
 171 #endif
 172   v = a*d;
 173   w = b*c;
 174   tau += v + w;     /* Add in other second-order terms.  */
 175   u = t + tau;
 176
 177   /* Construct long double result.  */
 178   if (nonfinite (u))
 179     return u;
 180   z.dval[0] = u;
 181   z.dval[1] = (t - u) + tau;
 182   return z.ldval;
 183 }
 184
 185 long double
 186 __gcc_qdiv (double a, double b, double c, double d)
 187 {
 188   longDblUnion z;
 189   double s, sigma, t, tau, u, v, w;
 190
 191   t = a / c;                    /* highest order double term */
 192
 193   if (unlikely (t == 0)         /* Preserve -0.  */
 194       || nonfinite (t))
 195     return t;
 196
 197   /* Finite nonzero result requires corrections to the highest order term.  */
 198
 199   s = c * t;                    /* (s,sigma) = c*t exactly.  */
 200   w = -(-b + d * t);    /* Written to get fnmsub for speed, but not
 201                            numerically necessary.  */
 202
 203   /* Use fused multiply-add to get low part of c * t.    */
 204 #ifndef _SOFT_FLOAT
 205   asm ("fmsub %0,%1,%2,%3" : "=f"(sigma) : "f"(c), "f"(t), "f"(s));
 206 #else
 207   sigma = fmsub (c, t, s);
 208 #endif
 209   v = a - s;
 210
 211   tau = ((v-sigma)+w)/c;   /* Correction to t.  */
 212   u = t + tau;
 213
 214   /* Construct long double result.  */
 215   if (nonfinite (u))
 216     return u;
 217   z.dval[0] = u;
 218   z.dval[1] = (t - u) + tau;
 219   return z.ldval;
 220 }
 221
 222 #if defined (_SOFT_FLOAT) && defined (__LONG_DOUBLE_128__)
 223
 224 long double __gcc_qneg (double, double);
 225 int __gcc_qeq (double, double, double, double);
 226 int __gcc_qne (double, double, double, double);
 227 int __gcc_qge (double, double, double, double);
 228 int __gcc_qle (double, double, double, double);
 229 int __gcc_qunord (double, double, double, double);
 230 long double __gcc_stoq (float);
 231 long double __gcc_dtoq (double);
 232 float __gcc_qtos (double, double);
 233 double __gcc_qtod (double, double);
 234 int __gcc_qtoi (double, double);
 235 unsigned int __gcc_qtou (double, double);
 236 long double __gcc_itoq (int);
 237 long double __gcc_utoq (unsigned int);
 238
 239 extern int __eqdf2 (double, double);
 240 extern int __ledf2 (double, double);
 241 extern int __gedf2 (double, double);
 242 extern int __unorddf2 (double, double);
 243
 244 /* Negate 'long double' value and return the result.    */
 245 long double
 246 __gcc_qneg (double a, double aa)
 247 {
 248   longDblUnion x;
 249
 250   x.dval[0] = -a;
 251   x.dval[1] = -aa;
 252   return x.ldval;
 253 }
 254
 255 /* Compare two 'long double' values for equality.  */
 256 int
 257 __gcc_qeq (double a, double aa, double c, double cc)
 258 {
 259   if (__eqdf2 (a, c) == 0)
 260     return __eqdf2 (aa, cc);
 261   return 1;
 262 }
 263
 264 strong_alias (__gcc_qeq, __gcc_qne);
 265
 266 /* Compare two 'long double' values for less than or equal.  */
 267 int
 268 __gcc_qle (double a, double aa, double c, double cc)
 269 {
 270   if (__eqdf2 (a, c) == 0)
 271     return __ledf2 (aa, cc);
 272   return __ledf2 (a, c);
 273 }
 274
 275 strong_alias (__gcc_qle, __gcc_qlt);
 276
 277 /* Compare two 'long double' values for greater than or equal.  */
 278 int
 279 __gcc_qge (double a, double aa, double c, double cc)
 280 {
 281   if (__eqdf2 (a, c) == 0)
 282     return __gedf2 (aa, cc);
 283   return __gedf2 (a, c);
 284 }
 285
 286 strong_alias (__gcc_qge, __gcc_qgt);
 287
 288 /* Compare two 'long double' values for unordered.  */
 289 int
 290 __gcc_qunord (double a, double aa, double c, double cc)
 291 {
 292   if (__eqdf2 (a, c) == 0)
 293     return __unorddf2 (aa, cc);
 294   return __unorddf2 (a, c);
 295 }
 296
 297 /* Convert single to long double.  */
 298 long double
 299 __gcc_stoq (float a)
 300 {
 301   longDblUnion x;
 302
 303   x.dval[0] = (double) a;
 304   x.dval[1] = 0.0;
 305
 306   return x.ldval;
 307 }
 308
 309 /* Convert double to long double.  */
 310 long double
 311 __gcc_dtoq (double a)
 312 {
 313   longDblUnion x;
 314
 315   x.dval[0] = a;
 316   x.dval[1] = 0.0;
 317
 318   return x.ldval;
 319 }
 320
 321 /* Convert long double to single.  */
 322 float
 323 __gcc_qtos (double a, double aa __attribute__ ((__unused__)))
 324 {
 325   return (float) a;
 326 }
 327
 328 /* Convert long double to double.  */
 329 double
 330 __gcc_qtod (double a, double aa __attribute__ ((__unused__)))
 331 {
 332   return a;
 333 }
 334
 335 /* Convert long double to int.  */
 336 int
 337 __gcc_qtoi (double a, double aa)
 338 {
 339   double z = a + aa;
 340   return (int) z;
 341 }
 342
 343 /* Convert long double to unsigned int.  */
 344 unsigned int
 345 __gcc_qtou (double a, double aa)
 346 {
 347   double z = a + aa;
 348   return (unsigned int) z;
 349 }
 350
 351 /* Convert int to long double.  */
 352 long double
 353 __gcc_itoq (int a)
 354 {
 355   return __gcc_dtoq ((double) a);
 356 }
 357
 358 /* Convert unsigned int to long double.  */
 359 long double
 360 __gcc_utoq (unsigned int a)
 361 {
 362   return __gcc_dtoq ((double) a);
 363 }
 364
 365 #include "config/soft-fp/soft-fp.h"
 366 #include "config/soft-fp/double.h"
 367 #include "config/soft-fp/quad.h"
 368
 369 /* Compute floating point multiply-subtract with higher (quad) precision.  */
 370 static double
 371 fmsub (double a, double b, double c)
 372 {
 373     FP_DECL_EX;
 374     FP_DECL_D(A);
 375     FP_DECL_D(B);
 376     FP_DECL_D(C);
 377     FP_DECL_Q(X);
 378     FP_DECL_Q(Y);
 379     FP_DECL_Q(Z);
 380     FP_DECL_Q(U);
 381     FP_DECL_Q(V);
 382     FP_DECL_D(R);
 383     double r;
 384     long double u, v, x, y, z;
 385
 386     FP_INIT_ROUNDMODE;
 387     FP_UNPACK_RAW_D (A, a);
 388     FP_UNPACK_RAW_D (B, b);
 389     FP_UNPACK_RAW_D (C, c);
 390
 391     /* Extend double to quad.  */
 392 #if (2 * _FP_W_TYPE_SIZE) < _FP_FRACBITS_Q
 393     FP_EXTEND(Q,D,4,2,X,A);
 394     FP_EXTEND(Q,D,4,2,Y,B);
 395     FP_EXTEND(Q,D,4,2,Z,C);
 396 #else
 397     FP_EXTEND(Q,D,2,1,X,A);
 398     FP_EXTEND(Q,D,2,1,Y,B);
 399     FP_EXTEND(Q,D,2,1,Z,C);
 400 #endif
 401     FP_PACK_RAW_Q(x,X);
 402     FP_PACK_RAW_Q(y,Y);
 403     FP_PACK_RAW_Q(z,Z);
 404     FP_HANDLE_EXCEPTIONS;
 405
 406     /* Multiply.  */
 407     FP_INIT_ROUNDMODE;
 408     FP_UNPACK_Q(X,x);
 409     FP_UNPACK_Q(Y,y);
 410     FP_MUL_Q(U,X,Y);
 411     FP_PACK_Q(u,U);
 412     FP_HANDLE_EXCEPTIONS;
 413
 414     /* Subtract.  */
 415     FP_INIT_ROUNDMODE;
 416     FP_UNPACK_SEMIRAW_Q(U,u);
 417     FP_UNPACK_SEMIRAW_Q(Z,z);
 418     FP_SUB_Q(V,U,Z);
 419     FP_PACK_SEMIRAW_Q(v,V);
 420     FP_HANDLE_EXCEPTIONS;
 421
 422     /* Truncate quad to double.  */
 423     FP_INIT_ROUNDMODE;
 424     FP_UNPACK_SEMIRAW_Q(V,v);
 425 #if (2 * _FP_W_TYPE_SIZE) < _FP_FRACBITS_Q
 426     FP_TRUNC(D,Q,2,4,R,V);
 427 #else
 428     FP_TRUNC(D,Q,1,2,R,V);
 429 #endif
 430     FP_PACK_SEMIRAW_D(r,R);
 431     FP_HANDLE_EXCEPTIONS;
 432
 433     return r;
 434 }
 435
 436 #endif
 437
 438 #endif