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