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#if !defined(CRT_INFINITY) && defined(HUGE_VAL) 9229135Sed#define CRT_INFINITY HUGE_VAL 10229135Sed#endif /* CRT_INFINITY */ 11214152Sed 12229135Sed#define makeFinite(x) { \ 13229135Sed (x).s.hi = crt_copysign(crt_isinf((x).s.hi) ? 1.0 : 0.0, (x).s.hi); \ 14229135Sed (x).s.lo = 0.0; \ 15229135Sed } 16214152Sed 17214152Sedlong double __gcc_qadd(long double, long double); 18214152Sedlong double __gcc_qsub(long double, long double); 19214152Sedlong double __gcc_qmul(long double, long double); 20214152Sedlong double __gcc_qdiv(long double, long double); 21214152Sed 22214152Sedlong double _Complex 23214152Sed__divtc3(long double a, long double b, long double c, long double d) 24214152Sed{ 25214152Sed DD cDD = { .ld = c }; 26214152Sed DD dDD = { .ld = d }; 27214152Sed 28214152Sed int ilogbw = 0; 29229135Sed const double logbw = crt_logb(crt_fmax(crt_fabs(cDD.s.hi), crt_fabs(dDD.s.hi) )); 30214152Sed 31229135Sed if (crt_isfinite(logbw)) 32214152Sed { 33214152Sed ilogbw = (int)logbw; 34214152Sed 35229135Sed cDD.s.hi = crt_scalbn(cDD.s.hi, -ilogbw); 36229135Sed cDD.s.lo = crt_scalbn(cDD.s.lo, -ilogbw); 37229135Sed dDD.s.hi = crt_scalbn(dDD.s.hi, -ilogbw); 38229135Sed dDD.s.lo = crt_scalbn(dDD.s.lo, -ilogbw); 39214152Sed } 40214152Sed 41214152Sed const long double denom = __gcc_qadd(__gcc_qmul(cDD.ld, cDD.ld), __gcc_qmul(dDD.ld, dDD.ld)); 42214152Sed const long double realNumerator = __gcc_qadd(__gcc_qmul(a,cDD.ld), __gcc_qmul(b,dDD.ld)); 43214152Sed const long double imagNumerator = __gcc_qsub(__gcc_qmul(b,cDD.ld), __gcc_qmul(a,dDD.ld)); 44214152Sed 45214152Sed DD real = { .ld = __gcc_qdiv(realNumerator, denom) }; 46214152Sed DD imag = { .ld = __gcc_qdiv(imagNumerator, denom) }; 47214152Sed 48229135Sed real.s.hi = crt_scalbn(real.s.hi, -ilogbw); 49229135Sed real.s.lo = crt_scalbn(real.s.lo, -ilogbw); 50229135Sed imag.s.hi = crt_scalbn(imag.s.hi, -ilogbw); 51229135Sed imag.s.lo = crt_scalbn(imag.s.lo, -ilogbw); 52214152Sed 53229135Sed if (crt_isnan(real.s.hi) && crt_isnan(imag.s.hi)) 54214152Sed { 55214152Sed DD aDD = { .ld = a }; 56214152Sed DD bDD = { .ld = b }; 57214152Sed DD rDD = { .ld = denom }; 58214152Sed 59229135Sed if ((rDD.s.hi == 0.0) && (!crt_isnan(aDD.s.hi) || 60229135Sed !crt_isnan(bDD.s.hi))) 61214152Sed { 62229135Sed real.s.hi = crt_copysign(CRT_INFINITY,cDD.s.hi) * aDD.s.hi; 63214152Sed real.s.lo = 0.0; 64229135Sed imag.s.hi = crt_copysign(CRT_INFINITY,cDD.s.hi) * bDD.s.hi; 65214152Sed imag.s.lo = 0.0; 66214152Sed } 67214152Sed 68229135Sed else if ((crt_isinf(aDD.s.hi) || crt_isinf(bDD.s.hi)) && 69229135Sed crt_isfinite(cDD.s.hi) && crt_isfinite(dDD.s.hi)) 70214152Sed { 71214152Sed makeFinite(aDD); 72214152Sed makeFinite(bDD); 73229135Sed real.s.hi = CRT_INFINITY * (aDD.s.hi*cDD.s.hi + bDD.s.hi*dDD.s.hi); 74214152Sed real.s.lo = 0.0; 75229135Sed imag.s.hi = CRT_INFINITY * (bDD.s.hi*cDD.s.hi - aDD.s.hi*dDD.s.hi); 76214152Sed imag.s.lo = 0.0; 77214152Sed } 78214152Sed 79229135Sed else if ((crt_isinf(cDD.s.hi) || crt_isinf(dDD.s.hi)) && 80229135Sed crt_isfinite(aDD.s.hi) && crt_isfinite(bDD.s.hi)) 81214152Sed { 82214152Sed makeFinite(cDD); 83214152Sed makeFinite(dDD); 84229135Sed real.s.hi = crt_copysign(0.0,(aDD.s.hi*cDD.s.hi + bDD.s.hi*dDD.s.hi)); 85214152Sed real.s.lo = 0.0; 86229135Sed imag.s.hi = crt_copysign(0.0,(bDD.s.hi*cDD.s.hi - aDD.s.hi*dDD.s.hi)); 87214152Sed imag.s.lo = 0.0; 88214152Sed } 89214152Sed } 90214152Sed 91214152Sed long double _Complex z; 92214152Sed __real__ z = real.ld; 93214152Sed __imag__ z = imag.ld; 94214152Sed 95214152Sed return z; 96214152Sed} 97