1 // SPDX-License-Identifier: GPL-2.0+
5 /* 128-bit long double support routines for Darwin.
13  * Implementations of floating-point long double basic arithmetic
17  * Double-double algorithms are based on the paper "Doubled-Precision
18  * IEEE Standard 754 Floating-Point Arithmetic" by W. Kahan, February 26,
20  * Doubled-Precision Floating-Point Computations", by Seppo Linnainmaa,
21  * ACM TOMS vol 7 no 3, September 1981, pages 272-283.
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.
34  * This code currently assumes big-endian.
41 #define nonfinite(a) unlikely(!isless(fabs(a), inf()))  argument
49 long double __gcc_qadd(double a, double aa, double c, double cc)  in __gcc_qadd()  argument
52 	double z, q, zz, xh;  in __gcc_qadd()  local
54 	z = a + c;  in __gcc_qadd()
56 	if (nonfinite(z)) {  in __gcc_qadd()
57 		z = cc + aa + c + a;  in __gcc_qadd()
58 		if (nonfinite(z))  in __gcc_qadd()
59 			return z;  in __gcc_qadd()
60 		x.dval[0] = z;	/* Will always be DBL_MAX.  */  in __gcc_qadd()
62 		if (fabs(a) > fabs(c))  in __gcc_qadd()
63 			x.dval[1] = a - z + c + zz;  in __gcc_qadd()
65 			x.dval[1] = c - z + a + zz;  in __gcc_qadd()
67 		q = a - z;  in __gcc_qadd()
68 		zz = q + c + (a - (q + z)) + aa + cc;  in __gcc_qadd()
70 		/* Keep -0 result.  */  in __gcc_qadd()
72 			return z;  in __gcc_qadd()
74 		xh = z + zz;  in __gcc_qadd()
79 		x.dval[1] = z - xh + zz;  in __gcc_qadd()
84 long double __gcc_qsub(double a, double b, double c, double d)  in __gcc_qsub()  argument
86 	return __gcc_qadd(a, b, -c, -d);  in __gcc_qsub()
89 long double __gcc_qmul(double a, double b, double c, double d)  in __gcc_qmul()  argument
91 	longDblUnion z;  in __gcc_qmul()  local
94 	t = a * c;		/* Highest order double term.  */  in __gcc_qmul()
96 	if (unlikely(t == 0)	/* Preserve -0.  */  in __gcc_qmul()
102 	/* Use fused multiply-add to get low part of a * c.  */  in __gcc_qmul()
104 	asm("fmsub %0,%1,%2,%3" : "=f"(tau) : "f"(a), "f"(c), "f"(t));  in __gcc_qmul()
106 	tau = fmsub(a, c, t);  in __gcc_qmul()
108 	v = a * d;  in __gcc_qmul()
110 	tau += v + w;		/* Add in other second-order terms.  */  in __gcc_qmul()
116 	z.dval[0] = u;  in __gcc_qmul()
117 	z.dval[1] = (t - u) + tau;  in __gcc_qmul()
118 	return z.ldval;  in __gcc_qmul()