gcc/config/rs6000/darwin-ldouble.c

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