libgcc/config/rs6000/ibm-ldouble.c

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