1214152Sed/* This file is distributed under the University of Illinois Open Source 2214152Sed * License. See LICENSE.TXT for details. 3214152Sed */ 4214152Sed 5214152Sed#include "DD.h" 6229135Sed#include "../int_math.h" 7214152Sed 8229135Sed#define makeFinite(x) { \ 9229135Sed (x).s.hi = crt_copysign(crt_isinf((x).s.hi) ? 1.0 : 0.0, (x).s.hi); \ 10229135Sed (x).s.lo = 0.0; \ 11229135Sed } 12214152Sed 13229135Sed#define zeroNaN(x) { \ 14229135Sed if (crt_isnan((x).s.hi)) { \ 15229135Sed (x).s.hi = crt_copysign(0.0, (x).s.hi); \ 16229135Sed (x).s.lo = 0.0; \ 17229135Sed } \ 18229135Sed } 19214152Sed 20214152Sedlong double __gcc_qadd(long double, long double); 21214152Sedlong double __gcc_qsub(long double, long double); 22214152Sedlong double __gcc_qmul(long double, long double); 23214152Sed 24214152Sedlong double _Complex 25214152Sed__multc3(long double a, long double b, long double c, long double d) 26214152Sed{ 27214152Sed long double ac = __gcc_qmul(a,c); 28214152Sed long double bd = __gcc_qmul(b,d); 29214152Sed long double ad = __gcc_qmul(a,d); 30214152Sed long double bc = __gcc_qmul(b,c); 31214152Sed 32214152Sed DD real = { .ld = __gcc_qsub(ac,bd) }; 33214152Sed DD imag = { .ld = __gcc_qadd(ad,bc) }; 34214152Sed 35229135Sed if (crt_isnan(real.s.hi) && crt_isnan(imag.s.hi)) 36214152Sed { 37214152Sed int recalc = 0; 38214152Sed 39214152Sed DD aDD = { .ld = a }; 40214152Sed DD bDD = { .ld = b }; 41214152Sed DD cDD = { .ld = c }; 42214152Sed DD dDD = { .ld = d }; 43214152Sed 44229135Sed if (crt_isinf(aDD.s.hi) || crt_isinf(bDD.s.hi)) 45214152Sed { 46214152Sed makeFinite(aDD); 47214152Sed makeFinite(bDD); 48214152Sed zeroNaN(cDD); 49214152Sed zeroNaN(dDD); 50214152Sed recalc = 1; 51214152Sed } 52214152Sed 53229135Sed if (crt_isinf(cDD.s.hi) || crt_isinf(dDD.s.hi)) 54214152Sed { 55214152Sed makeFinite(cDD); 56214152Sed makeFinite(dDD); 57214152Sed zeroNaN(aDD); 58214152Sed zeroNaN(bDD); 59214152Sed recalc = 1; 60214152Sed } 61214152Sed 62214152Sed if (!recalc) 63214152Sed { 64214152Sed DD acDD = { .ld = ac }; 65214152Sed DD bdDD = { .ld = bd }; 66214152Sed DD adDD = { .ld = ad }; 67214152Sed DD bcDD = { .ld = bc }; 68214152Sed 69229135Sed if (crt_isinf(acDD.s.hi) || crt_isinf(bdDD.s.hi) || 70229135Sed crt_isinf(adDD.s.hi) || crt_isinf(bcDD.s.hi)) 71214152Sed { 72214152Sed zeroNaN(aDD); 73214152Sed zeroNaN(bDD); 74214152Sed zeroNaN(cDD); 75214152Sed zeroNaN(dDD); 76214152Sed recalc = 1; 77214152Sed } 78214152Sed } 79214152Sed 80214152Sed if (recalc) 81214152Sed { 82229135Sed real.s.hi = CRT_INFINITY * (aDD.s.hi*cDD.s.hi - bDD.s.hi*dDD.s.hi); 83214152Sed real.s.lo = 0.0; 84229135Sed imag.s.hi = CRT_INFINITY * (aDD.s.hi*dDD.s.hi + bDD.s.hi*cDD.s.hi); 85214152Sed imag.s.lo = 0.0; 86214152Sed } 87214152Sed } 88214152Sed 89214152Sed long double _Complex z; 90214152Sed __real__ z = real.ld; 91214152Sed __imag__ z = imag.ld; 92214152Sed 93214152Sed return z; 94214152Sed} 95