1 /* 128-bit long double support routines for Darwin. 2 Copyright (C) 1993-2019 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 455