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 445