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