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