libgcc/config/rs6000/ibm-ldouble.c

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