libgcc/config/rs6000/ibm-ldouble.c

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