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