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