1 /* 128-bit long double support routines for Darwin. 2 Copyright (C) 1993-2022 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 static inline IBM128_TYPE 122 ldouble_qadd_internal (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_qadd (double a, double aa, double c, double cc) 162 { 163 return ldouble_qadd_internal (a, aa, c, cc); 164 } 165 166 IBM128_TYPE 167 __gcc_qsub (double a, double aa, double c, double cc) 168 { 169 return ldouble_qadd_internal (a, aa, -c, -cc); 170 } 171 172 #ifdef __NO_FPRS__ 173 static double fmsub (double, double, double); 174 #endif 175 176 IBM128_TYPE 177 __gcc_qmul (double a, double b, double c, double d) 178 { 179 double xh, xl, t, tau, u, v, w; 180 181 t = a * c; /* Highest order double term. */ 182 183 if (unlikely (t == 0) /* Preserve -0. */ 184 || nonfinite (t)) 185 return t; 186 187 /* Sum terms of two highest orders. */ 188 189 /* Use fused multiply-add to get low part of a * c. */ 190 #ifndef __NO_FPRS__ 191 asm ("fmsub %0,%1,%2,%3" : "=f"(tau) : "f"(a), "f"(c), "f"(t)); 192 #else 193 tau = fmsub (a, c, t); 194 #endif 195 v = a*d; 196 w = b*c; 197 tau += v + w; /* Add in other second-order terms. */ 198 u = t + tau; 199 200 /* Construct IBM128_TYPE result. */ 201 if (nonfinite (u)) 202 return u; 203 xh = u; 204 xl = (t - u) + tau; 205 return pack_ldouble (xh, xl); 206 } 207 208 IBM128_TYPE 209 __gcc_qdiv (double a, double b, double c, double d) 210 { 211 double xh, xl, s, sigma, t, tau, u, v, w; 212 213 t = a / c; /* highest order double term */ 214 215 if (unlikely (t == 0) /* Preserve -0. */ 216 || nonfinite (t)) 217 return t; 218 219 /* Finite nonzero result requires corrections to the highest order 220 term. These corrections require the low part of c * t to be 221 exactly represented in double. */ 222 if (fabs (a) <= 0x1p-969) 223 { 224 a *= 0x1p106; 225 b *= 0x1p106; 226 c *= 0x1p106; 227 d *= 0x1p106; 228 } 229 230 s = c * t; /* (s,sigma) = c*t exactly. */ 231 w = -(-b + d * t); /* Written to get fnmsub for speed, but not 232 numerically necessary. */ 233 234 /* Use fused multiply-add to get low part of c * t. */ 235 #ifndef __NO_FPRS__ 236 asm ("fmsub %0,%1,%2,%3" : "=f"(sigma) : "f"(c), "f"(t), "f"(s)); 237 #else 238 sigma = fmsub (c, t, s); 239 #endif 240 v = a - s; 241 242 tau = ((v-sigma)+w)/c; /* Correction to t. */ 243 u = t + tau; 244 245 /* Construct IBM128_TYPE result. */ 246 if (nonfinite (u)) 247 return u; 248 xh = u; 249 xl = (t - u) + tau; 250 return pack_ldouble (xh, xl); 251 } 252 253 #if defined (_SOFT_DOUBLE) && defined (__LONG_DOUBLE_128__) 254 255 IBM128_TYPE __gcc_qneg (double, double); 256 int __gcc_qeq (double, double, double, double); 257 int __gcc_qne (double, double, double, double); 258 int __gcc_qge (double, double, double, double); 259 int __gcc_qle (double, double, double, double); 260 IBM128_TYPE __gcc_stoq (float); 261 IBM128_TYPE __gcc_dtoq (double); 262 float __gcc_qtos (double, double); 263 double __gcc_qtod (double, double); 264 int __gcc_qtoi (double, double); 265 unsigned int __gcc_qtou (double, double); 266 IBM128_TYPE __gcc_itoq (int); 267 IBM128_TYPE __gcc_utoq (unsigned int); 268 269 extern int __eqdf2 (double, double); 270 extern int __ledf2 (double, double); 271 extern int __gedf2 (double, double); 272 273 /* Negate 'IBM128_TYPE' value and return the result. */ 274 IBM128_TYPE 275 __gcc_qneg (double a, double aa) 276 { 277 return pack_ldouble (-a, -aa); 278 } 279 280 /* Compare two 'IBM128_TYPE' values for equality. */ 281 int 282 __gcc_qeq (double a, double aa, double c, double cc) 283 { 284 if (__eqdf2 (a, c) == 0) 285 return __eqdf2 (aa, cc); 286 return 1; 287 } 288 289 strong_alias (__gcc_qeq, __gcc_qne); 290 291 /* Compare two 'IBM128_TYPE' values for less than or equal. */ 292 int 293 __gcc_qle (double a, double aa, double c, double cc) 294 { 295 if (__eqdf2 (a, c) == 0) 296 return __ledf2 (aa, cc); 297 return __ledf2 (a, c); 298 } 299 300 strong_alias (__gcc_qle, __gcc_qlt); 301 302 /* Compare two 'IBM128_TYPE' values for greater than or equal. */ 303 int 304 __gcc_qge (double a, double aa, double c, double cc) 305 { 306 if (__eqdf2 (a, c) == 0) 307 return __gedf2 (aa, cc); 308 return __gedf2 (a, c); 309 } 310 311 strong_alias (__gcc_qge, __gcc_qgt); 312 313 /* Convert single to IBM128_TYPE. */ 314 IBM128_TYPE 315 __gcc_stoq (float a) 316 { 317 return pack_ldouble ((double) a, 0.0); 318 } 319 320 /* Convert double to IBM128_TYPE. */ 321 IBM128_TYPE 322 __gcc_dtoq (double a) 323 { 324 return pack_ldouble (a, 0.0); 325 } 326 327 /* Convert IBM128_TYPE to single. */ 328 float 329 __gcc_qtos (double a, double aa __attribute__ ((__unused__))) 330 { 331 return (float) a; 332 } 333 334 /* Convert IBM128_TYPE to double. */ 335 double 336 __gcc_qtod (double a, double aa __attribute__ ((__unused__))) 337 { 338 return a; 339 } 340 341 /* Convert IBM128_TYPE to int. */ 342 int 343 __gcc_qtoi (double a, double aa) 344 { 345 double z = a + aa; 346 return (int) z; 347 } 348 349 /* Convert IBM128_TYPE to unsigned int. */ 350 unsigned int 351 __gcc_qtou (double a, double aa) 352 { 353 double z = a + aa; 354 return (unsigned int) z; 355 } 356 357 /* Convert int to IBM128_TYPE. */ 358 IBM128_TYPE 359 __gcc_itoq (int a) 360 { 361 return __gcc_dtoq ((double) a); 362 } 363 364 /* Convert unsigned int to IBM128_TYPE. */ 365 IBM128_TYPE 366 __gcc_utoq (unsigned int a) 367 { 368 return __gcc_dtoq ((double) a); 369 } 370 371 #endif 372 373 #ifdef __NO_FPRS__ 374 375 int __gcc_qunord (double, double, double, double); 376 377 extern int __eqdf2 (double, double); 378 extern int __unorddf2 (double, double); 379 380 /* Compare two 'IBM128_TYPE' values for unordered. */ 381 int 382 __gcc_qunord (double a, double aa, double c, double cc) 383 { 384 if (__eqdf2 (a, c) == 0) 385 return __unorddf2 (aa, cc); 386 return __unorddf2 (a, c); 387 } 388 389 #include "soft-fp/soft-fp.h" 390 #include "soft-fp/double.h" 391 #include "soft-fp/quad.h" 392 393 /* Compute floating point multiply-subtract with higher (quad) precision. */ 394 static double 395 fmsub (double a, double b, double c) 396 { 397 FP_DECL_EX; 398 FP_DECL_D(A); 399 FP_DECL_D(B); 400 FP_DECL_D(C); 401 FP_DECL_Q(X); 402 FP_DECL_Q(Y); 403 FP_DECL_Q(Z); 404 FP_DECL_Q(U); 405 FP_DECL_Q(V); 406 FP_DECL_D(R); 407 double r; 408 IBM128_TYPE u, x, y, z; 409 410 FP_INIT_ROUNDMODE; 411 FP_UNPACK_RAW_D (A, a); 412 FP_UNPACK_RAW_D (B, b); 413 FP_UNPACK_RAW_D (C, c); 414 415 /* Extend double to quad. */ 416 #if _FP_W_TYPE_SIZE < 64 417 FP_EXTEND(Q,D,4,2,X,A); 418 FP_EXTEND(Q,D,4,2,Y,B); 419 FP_EXTEND(Q,D,4,2,Z,C); 420 #else 421 FP_EXTEND(Q,D,2,1,X,A); 422 FP_EXTEND(Q,D,2,1,Y,B); 423 FP_EXTEND(Q,D,2,1,Z,C); 424 #endif 425 FP_PACK_RAW_Q(x,X); 426 FP_PACK_RAW_Q(y,Y); 427 FP_PACK_RAW_Q(z,Z); 428 FP_HANDLE_EXCEPTIONS; 429 430 /* Multiply. */ 431 FP_INIT_ROUNDMODE; 432 FP_UNPACK_Q(X,x); 433 FP_UNPACK_Q(Y,y); 434 FP_MUL_Q(U,X,Y); 435 FP_PACK_Q(u,U); 436 FP_HANDLE_EXCEPTIONS; 437 438 /* Subtract. */ 439 FP_INIT_ROUNDMODE; 440 FP_UNPACK_SEMIRAW_Q(U,u); 441 FP_UNPACK_SEMIRAW_Q(Z,z); 442 FP_SUB_Q(V,U,Z); 443 444 /* Truncate quad to double. */ 445 #if _FP_W_TYPE_SIZE < 64 446 V_f[3] &= 0x0007ffff; 447 FP_TRUNC(D,Q,2,4,R,V); 448 #else 449 V_f1 &= 0x0007ffffffffffffL; 450 FP_TRUNC(D,Q,1,2,R,V); 451 #endif 452 FP_PACK_SEMIRAW_D(r,R); 453 FP_HANDLE_EXCEPTIONS; 454 455 return r; 456 } 457 458 #endif 459 460 #endif 461