1 /* 2 * Borrowed from GCC 4.2.2 (which still was GPL v2+) 3 */ 4 /* 128-bit long double support routines for Darwin. 5 Copyright (C) 1993, 2003, 2004, 2005, 2006, 2007 6 Free Software Foundation, Inc. 7 8 This file is part of GCC. 9 10 * SPDX-License-Identifier: GPL-2.0+ 11 */ 12 13 /* 14 * Implementations of floating-point long double basic arithmetic 15 * functions called by the IBM C compiler when generating code for 16 * PowerPC platforms. In particular, the following functions are 17 * implemented: __gcc_qadd, __gcc_qsub, __gcc_qmul, and __gcc_qdiv. 18 * Double-double algorithms are based on the paper "Doubled-Precision 19 * IEEE Standard 754 Floating-Point Arithmetic" by W. Kahan, February 26, 20 * 1987. An alternative published reference is "Software for 21 * Doubled-Precision Floating-Point Computations", by Seppo Linnainmaa, 22 * ACM TOMS vol 7 no 3, September 1981, pages 272-283. 23 */ 24 25 /* 26 * Each long double is made up of two IEEE doubles. The value of the 27 * long double is the sum of the values of the two parts. The most 28 * significant part is required to be the value of the long double 29 * rounded to the nearest double, as specified by IEEE. For Inf 30 * values, the least significant part is required to be one of +0.0 or 31 * -0.0. No other requirements are made; so, for example, 1.0 may be 32 * represented as (1.0, +0.0) or (1.0, -0.0), and the low part of a 33 * NaN is don't-care. 34 * 35 * This code currently assumes big-endian. 36 */ 37 38 #define fabs(x) __builtin_fabs(x) 39 #define isless(x, y) __builtin_isless(x, y) 40 #define inf() __builtin_inf() 41 #define unlikely(x) __builtin_expect((x), 0) 42 #define nonfinite(a) unlikely(!isless(fabs(a), inf())) 43 44 typedef union { 45 long double ldval; 46 double dval[2]; 47 } longDblUnion; 48 49 /* Add two 'long double' values and return the result. */ 50 long double __gcc_qadd(double a, double aa, double c, double cc) 51 { 52 longDblUnion x; 53 double z, q, zz, xh; 54 55 z = a + c; 56 57 if (nonfinite(z)) { 58 z = cc + aa + c + a; 59 if (nonfinite(z)) 60 return z; 61 x.dval[0] = z; /* Will always be DBL_MAX. */ 62 zz = aa + cc; 63 if (fabs(a) > fabs(c)) 64 x.dval[1] = a - z + c + zz; 65 else 66 x.dval[1] = c - z + a + zz; 67 } else { 68 q = a - z; 69 zz = q + c + (a - (q + z)) + aa + cc; 70 71 /* Keep -0 result. */ 72 if (zz == 0.0) 73 return z; 74 75 xh = z + zz; 76 if (nonfinite(xh)) 77 return xh; 78 79 x.dval[0] = xh; 80 x.dval[1] = z - xh + zz; 81 } 82 return x.ldval; 83 } 84 85 long double __gcc_qsub(double a, double b, double c, double d) 86 { 87 return __gcc_qadd(a, b, -c, -d); 88 } 89 90 long double __gcc_qmul(double a, double b, double c, double d) 91 { 92 longDblUnion z; 93 double t, tau, u, v, w; 94 95 t = a * c; /* Highest order double term. */ 96 97 if (unlikely(t == 0) /* Preserve -0. */ 98 || nonfinite(t)) 99 return t; 100 101 /* Sum terms of two highest orders. */ 102 103 /* Use fused multiply-add to get low part of a * c. */ 104 #ifndef __NO_FPRS__ 105 asm("fmsub %0,%1,%2,%3" : "=f"(tau) : "f"(a), "f"(c), "f"(t)); 106 #else 107 tau = fmsub(a, c, t); 108 #endif 109 v = a * d; 110 w = b * c; 111 tau += v + w; /* Add in other second-order terms. */ 112 u = t + tau; 113 114 /* Construct long double result. */ 115 if (nonfinite(u)) 116 return u; 117 z.dval[0] = u; 118 z.dval[1] = (t - u) + tau; 119 return z.ldval; 120 } 121