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