1/* 128-bit long double support routines for Darwin. 2 Copyright (C) 1993-2020 Free Software Foundation, Inc. 3 4This file is part of GCC. 5 6GCC is free software; you can redistribute it and/or modify it under 7the terms of the GNU General Public License as published by the Free 8Software Foundation; either version 3, or (at your option) any later 9version. 10 11GCC is distributed in the hope that it will be useful, but WITHOUT ANY 12WARRANTY; without even the implied warranty of MERCHANTABILITY or 13FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 14for more details. 15 16Under Section 7 of GPL version 3, you are granted additional 17permissions described in the GCC Runtime Library Exception, version 183.1, as published by the Free Software Foundation. 19 20You should have received a copy of the GNU General Public License and 21a copy of the GCC Runtime Library Exception along with this program; 22see 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 78IBM128_TYPE __gcc_qadd (double, double, double, double); 79IBM128_TYPE __gcc_qsub (double, double, double, double); 80IBM128_TYPE __gcc_qmul (double, double, double, double); 81IBM128_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. */ 102static inline IBM128_TYPE 103pack_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. */ 121IBM128_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 160IBM128_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__ 167static double fmsub (double, double, double); 168#endif 169 170IBM128_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 202IBM128_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 249IBM128_TYPE __gcc_qneg (double, double); 250int __gcc_qeq (double, double, double, double); 251int __gcc_qne (double, double, double, double); 252int __gcc_qge (double, double, double, double); 253int __gcc_qle (double, double, double, double); 254IBM128_TYPE __gcc_stoq (float); 255IBM128_TYPE __gcc_dtoq (double); 256float __gcc_qtos (double, double); 257double __gcc_qtod (double, double); 258int __gcc_qtoi (double, double); 259unsigned int __gcc_qtou (double, double); 260IBM128_TYPE __gcc_itoq (int); 261IBM128_TYPE __gcc_utoq (unsigned int); 262 263extern int __eqdf2 (double, double); 264extern int __ledf2 (double, double); 265extern int __gedf2 (double, double); 266 267/* Negate 'IBM128_TYPE' value and return the result. */ 268IBM128_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. */ 275int 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 283strong_alias (__gcc_qeq, __gcc_qne); 284 285/* Compare two 'IBM128_TYPE' values for less than or equal. */ 286int 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 294strong_alias (__gcc_qle, __gcc_qlt); 295 296/* Compare two 'IBM128_TYPE' values for greater than or equal. */ 297int 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 305strong_alias (__gcc_qge, __gcc_qgt); 306 307/* Convert single to IBM128_TYPE. */ 308IBM128_TYPE 309__gcc_stoq (float a) 310{ 311 return pack_ldouble ((double) a, 0.0); 312} 313 314/* Convert double to IBM128_TYPE. */ 315IBM128_TYPE 316__gcc_dtoq (double a) 317{ 318 return pack_ldouble (a, 0.0); 319} 320 321/* Convert IBM128_TYPE to single. */ 322float 323__gcc_qtos (double a, double aa __attribute__ ((__unused__))) 324{ 325 return (float) a; 326} 327 328/* Convert IBM128_TYPE to double. */ 329double 330__gcc_qtod (double a, double aa __attribute__ ((__unused__))) 331{ 332 return a; 333} 334 335/* Convert IBM128_TYPE to int. */ 336int 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. */ 344unsigned 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. */ 352IBM128_TYPE 353__gcc_itoq (int a) 354{ 355 return __gcc_dtoq ((double) a); 356} 357 358/* Convert unsigned int to IBM128_TYPE. */ 359IBM128_TYPE 360__gcc_utoq (unsigned int a) 361{ 362 return __gcc_dtoq ((double) a); 363} 364 365#endif 366 367#ifdef __NO_FPRS__ 368 369int __gcc_qunord (double, double, double, double); 370 371extern int __eqdf2 (double, double); 372extern int __unorddf2 (double, double); 373 374/* Compare two 'IBM128_TYPE' values for unordered. */ 375int 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. */ 388static double 389fmsub (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 _FP_W_TYPE_SIZE < 64 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 _FP_W_TYPE_SIZE < 64 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