Tom Rini | 83d290c | 2018-05-06 17:58:06 -0400 | [diff] [blame] | 1 | // SPDX-License-Identifier: GPL-2.0+ |
Wolfgang Denk | 139e187 | 2011-12-22 04:29:41 +0000 | [diff] [blame] | 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. |
Wolfgang Denk | 1a45966 | 2013-07-08 09:37:19 +0200 | [diff] [blame] | 10 | */ |
Wolfgang Denk | 139e187 | 2011-12-22 04:29:41 +0000 | [diff] [blame] | 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 | } |