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