blob: 16105c837bea526d730cf68b6855a7c2f9a47c92 [file] [log] [blame]
Tom Rini83d290c2018-05-06 17:58:06 -04001// SPDX-License-Identifier: GPL-2.0+
Wolfgang Denk139e1872011-12-22 04:29:41 +00002/*
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
9This file is part of GCC.
Wolfgang Denk1a459662013-07-08 09:37:19 +020010 */
Wolfgang Denk139e1872011-12-22 04:29:41 +000011
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
43typedef union {
44 long double ldval;
45 double dval[2];
46} longDblUnion;
47
48/* Add two 'long double' values and return the result. */
49long 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
84long double __gcc_qsub(double a, double b, double c, double d)
85{
86 return __gcc_qadd(a, b, -c, -d);
87}
88
89long 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}