1*0a6a1f1dSLionel Sambuc /* This file is distributed under the University of Illinois Open Source
2*0a6a1f1dSLionel Sambuc * License. See LICENSE.TXT for details.
3*0a6a1f1dSLionel Sambuc */
4*0a6a1f1dSLionel Sambuc
5*0a6a1f1dSLionel Sambuc #include "DD.h"
6*0a6a1f1dSLionel Sambuc #include "../int_math.h"
7*0a6a1f1dSLionel Sambuc
8*0a6a1f1dSLionel Sambuc #define makeFinite(x) { \
9*0a6a1f1dSLionel Sambuc (x).s.hi = crt_copysign(crt_isinf((x).s.hi) ? 1.0 : 0.0, (x).s.hi); \
10*0a6a1f1dSLionel Sambuc (x).s.lo = 0.0; \
11*0a6a1f1dSLionel Sambuc }
12*0a6a1f1dSLionel Sambuc
13*0a6a1f1dSLionel Sambuc #define zeroNaN(x) { \
14*0a6a1f1dSLionel Sambuc if (crt_isnan((x).s.hi)) { \
15*0a6a1f1dSLionel Sambuc (x).s.hi = crt_copysign(0.0, (x).s.hi); \
16*0a6a1f1dSLionel Sambuc (x).s.lo = 0.0; \
17*0a6a1f1dSLionel Sambuc } \
18*0a6a1f1dSLionel Sambuc }
19*0a6a1f1dSLionel Sambuc
20*0a6a1f1dSLionel Sambuc long double __gcc_qadd(long double, long double);
21*0a6a1f1dSLionel Sambuc long double __gcc_qsub(long double, long double);
22*0a6a1f1dSLionel Sambuc long double __gcc_qmul(long double, long double);
23*0a6a1f1dSLionel Sambuc
24*0a6a1f1dSLionel Sambuc long double _Complex
__multc3(long double a,long double b,long double c,long double d)25*0a6a1f1dSLionel Sambuc __multc3(long double a, long double b, long double c, long double d)
26*0a6a1f1dSLionel Sambuc {
27*0a6a1f1dSLionel Sambuc long double ac = __gcc_qmul(a,c);
28*0a6a1f1dSLionel Sambuc long double bd = __gcc_qmul(b,d);
29*0a6a1f1dSLionel Sambuc long double ad = __gcc_qmul(a,d);
30*0a6a1f1dSLionel Sambuc long double bc = __gcc_qmul(b,c);
31*0a6a1f1dSLionel Sambuc
32*0a6a1f1dSLionel Sambuc DD real = { .ld = __gcc_qsub(ac,bd) };
33*0a6a1f1dSLionel Sambuc DD imag = { .ld = __gcc_qadd(ad,bc) };
34*0a6a1f1dSLionel Sambuc
35*0a6a1f1dSLionel Sambuc if (crt_isnan(real.s.hi) && crt_isnan(imag.s.hi))
36*0a6a1f1dSLionel Sambuc {
37*0a6a1f1dSLionel Sambuc int recalc = 0;
38*0a6a1f1dSLionel Sambuc
39*0a6a1f1dSLionel Sambuc DD aDD = { .ld = a };
40*0a6a1f1dSLionel Sambuc DD bDD = { .ld = b };
41*0a6a1f1dSLionel Sambuc DD cDD = { .ld = c };
42*0a6a1f1dSLionel Sambuc DD dDD = { .ld = d };
43*0a6a1f1dSLionel Sambuc
44*0a6a1f1dSLionel Sambuc if (crt_isinf(aDD.s.hi) || crt_isinf(bDD.s.hi))
45*0a6a1f1dSLionel Sambuc {
46*0a6a1f1dSLionel Sambuc makeFinite(aDD);
47*0a6a1f1dSLionel Sambuc makeFinite(bDD);
48*0a6a1f1dSLionel Sambuc zeroNaN(cDD);
49*0a6a1f1dSLionel Sambuc zeroNaN(dDD);
50*0a6a1f1dSLionel Sambuc recalc = 1;
51*0a6a1f1dSLionel Sambuc }
52*0a6a1f1dSLionel Sambuc
53*0a6a1f1dSLionel Sambuc if (crt_isinf(cDD.s.hi) || crt_isinf(dDD.s.hi))
54*0a6a1f1dSLionel Sambuc {
55*0a6a1f1dSLionel Sambuc makeFinite(cDD);
56*0a6a1f1dSLionel Sambuc makeFinite(dDD);
57*0a6a1f1dSLionel Sambuc zeroNaN(aDD);
58*0a6a1f1dSLionel Sambuc zeroNaN(bDD);
59*0a6a1f1dSLionel Sambuc recalc = 1;
60*0a6a1f1dSLionel Sambuc }
61*0a6a1f1dSLionel Sambuc
62*0a6a1f1dSLionel Sambuc if (!recalc)
63*0a6a1f1dSLionel Sambuc {
64*0a6a1f1dSLionel Sambuc DD acDD = { .ld = ac };
65*0a6a1f1dSLionel Sambuc DD bdDD = { .ld = bd };
66*0a6a1f1dSLionel Sambuc DD adDD = { .ld = ad };
67*0a6a1f1dSLionel Sambuc DD bcDD = { .ld = bc };
68*0a6a1f1dSLionel Sambuc
69*0a6a1f1dSLionel Sambuc if (crt_isinf(acDD.s.hi) || crt_isinf(bdDD.s.hi) ||
70*0a6a1f1dSLionel Sambuc crt_isinf(adDD.s.hi) || crt_isinf(bcDD.s.hi))
71*0a6a1f1dSLionel Sambuc {
72*0a6a1f1dSLionel Sambuc zeroNaN(aDD);
73*0a6a1f1dSLionel Sambuc zeroNaN(bDD);
74*0a6a1f1dSLionel Sambuc zeroNaN(cDD);
75*0a6a1f1dSLionel Sambuc zeroNaN(dDD);
76*0a6a1f1dSLionel Sambuc recalc = 1;
77*0a6a1f1dSLionel Sambuc }
78*0a6a1f1dSLionel Sambuc }
79*0a6a1f1dSLionel Sambuc
80*0a6a1f1dSLionel Sambuc if (recalc)
81*0a6a1f1dSLionel Sambuc {
82*0a6a1f1dSLionel Sambuc real.s.hi = CRT_INFINITY * (aDD.s.hi*cDD.s.hi - bDD.s.hi*dDD.s.hi);
83*0a6a1f1dSLionel Sambuc real.s.lo = 0.0;
84*0a6a1f1dSLionel Sambuc imag.s.hi = CRT_INFINITY * (aDD.s.hi*dDD.s.hi + bDD.s.hi*cDD.s.hi);
85*0a6a1f1dSLionel Sambuc imag.s.lo = 0.0;
86*0a6a1f1dSLionel Sambuc }
87*0a6a1f1dSLionel Sambuc }
88*0a6a1f1dSLionel Sambuc
89*0a6a1f1dSLionel Sambuc long double _Complex z;
90*0a6a1f1dSLionel Sambuc __real__ z = real.ld;
91*0a6a1f1dSLionel Sambuc __imag__ z = imag.ld;
92*0a6a1f1dSLionel Sambuc
93*0a6a1f1dSLionel Sambuc return z;
94*0a6a1f1dSLionel Sambuc }
95