1 1.1 joerg /* This file is distributed under the University of Illinois Open Source 2 1.1 joerg * License. See LICENSE.TXT for details. 3 1.1 joerg */ 4 1.1 joerg 5 1.1 joerg #include "DD.h" 6 1.1 joerg #include "../int_math.h" 7 1.1 joerg 8 1.1 joerg #if !defined(CRT_INFINITY) && defined(HUGE_VAL) 9 1.1 joerg #define CRT_INFINITY HUGE_VAL 10 1.1 joerg #endif /* CRT_INFINITY */ 11 1.1 joerg 12 1.1 joerg #define makeFinite(x) { \ 13 1.1 joerg (x).s.hi = crt_copysign(crt_isinf((x).s.hi) ? 1.0 : 0.0, (x).s.hi); \ 14 1.1 joerg (x).s.lo = 0.0; \ 15 1.1 joerg } 16 1.1 joerg 17 1.1 joerg long double _Complex 18 1.1 joerg __divtc3(long double a, long double b, long double c, long double d) 19 1.1 joerg { 20 1.1 joerg DD cDD = { .ld = c }; 21 1.1 joerg DD dDD = { .ld = d }; 22 1.1 joerg 23 1.1 joerg int ilogbw = 0; 24 1.1 joerg const double logbw = crt_logb(crt_fmax(crt_fabs(cDD.s.hi), crt_fabs(dDD.s.hi) )); 25 1.1 joerg 26 1.1 joerg if (crt_isfinite(logbw)) 27 1.1 joerg { 28 1.1 joerg ilogbw = (int)logbw; 29 1.1 joerg 30 1.1 joerg cDD.s.hi = crt_scalbn(cDD.s.hi, -ilogbw); 31 1.1 joerg cDD.s.lo = crt_scalbn(cDD.s.lo, -ilogbw); 32 1.1 joerg dDD.s.hi = crt_scalbn(dDD.s.hi, -ilogbw); 33 1.1 joerg dDD.s.lo = crt_scalbn(dDD.s.lo, -ilogbw); 34 1.1 joerg } 35 1.1 joerg 36 1.1 joerg const long double denom = __gcc_qadd(__gcc_qmul(cDD.ld, cDD.ld), __gcc_qmul(dDD.ld, dDD.ld)); 37 1.1 joerg const long double realNumerator = __gcc_qadd(__gcc_qmul(a,cDD.ld), __gcc_qmul(b,dDD.ld)); 38 1.1 joerg const long double imagNumerator = __gcc_qsub(__gcc_qmul(b,cDD.ld), __gcc_qmul(a,dDD.ld)); 39 1.1 joerg 40 1.1 joerg DD real = { .ld = __gcc_qdiv(realNumerator, denom) }; 41 1.1 joerg DD imag = { .ld = __gcc_qdiv(imagNumerator, denom) }; 42 1.1 joerg 43 1.1 joerg real.s.hi = crt_scalbn(real.s.hi, -ilogbw); 44 1.1 joerg real.s.lo = crt_scalbn(real.s.lo, -ilogbw); 45 1.1 joerg imag.s.hi = crt_scalbn(imag.s.hi, -ilogbw); 46 1.1 joerg imag.s.lo = crt_scalbn(imag.s.lo, -ilogbw); 47 1.1 joerg 48 1.1 joerg if (crt_isnan(real.s.hi) && crt_isnan(imag.s.hi)) 49 1.1 joerg { 50 1.1 joerg DD aDD = { .ld = a }; 51 1.1 joerg DD bDD = { .ld = b }; 52 1.1 joerg DD rDD = { .ld = denom }; 53 1.1 joerg 54 1.1 joerg if ((rDD.s.hi == 0.0) && (!crt_isnan(aDD.s.hi) || 55 1.1 joerg !crt_isnan(bDD.s.hi))) 56 1.1 joerg { 57 1.1 joerg real.s.hi = crt_copysign(CRT_INFINITY,cDD.s.hi) * aDD.s.hi; 58 1.1 joerg real.s.lo = 0.0; 59 1.1 joerg imag.s.hi = crt_copysign(CRT_INFINITY,cDD.s.hi) * bDD.s.hi; 60 1.1 joerg imag.s.lo = 0.0; 61 1.1 joerg } 62 1.1 joerg 63 1.1 joerg else if ((crt_isinf(aDD.s.hi) || crt_isinf(bDD.s.hi)) && 64 1.1 joerg crt_isfinite(cDD.s.hi) && crt_isfinite(dDD.s.hi)) 65 1.1 joerg { 66 1.1 joerg makeFinite(aDD); 67 1.1 joerg makeFinite(bDD); 68 1.1 joerg real.s.hi = CRT_INFINITY * (aDD.s.hi*cDD.s.hi + bDD.s.hi*dDD.s.hi); 69 1.1 joerg real.s.lo = 0.0; 70 1.1 joerg imag.s.hi = CRT_INFINITY * (bDD.s.hi*cDD.s.hi - aDD.s.hi*dDD.s.hi); 71 1.1 joerg imag.s.lo = 0.0; 72 1.1 joerg } 73 1.1 joerg 74 1.1 joerg else if ((crt_isinf(cDD.s.hi) || crt_isinf(dDD.s.hi)) && 75 1.1 joerg crt_isfinite(aDD.s.hi) && crt_isfinite(bDD.s.hi)) 76 1.1 joerg { 77 1.1 joerg makeFinite(cDD); 78 1.1 joerg makeFinite(dDD); 79 1.1 joerg real.s.hi = crt_copysign(0.0,(aDD.s.hi*cDD.s.hi + bDD.s.hi*dDD.s.hi)); 80 1.1 joerg real.s.lo = 0.0; 81 1.1 joerg imag.s.hi = crt_copysign(0.0,(bDD.s.hi*cDD.s.hi - aDD.s.hi*dDD.s.hi)); 82 1.1 joerg imag.s.lo = 0.0; 83 1.1 joerg } 84 1.1 joerg } 85 1.1 joerg 86 1.1 joerg long double _Complex z; 87 1.1 joerg __real__ z = real.ld; 88 1.1 joerg __imag__ z = imag.ld; 89 1.1 joerg 90 1.1 joerg return z; 91 1.1 joerg } 92
Indexes created Sun Sep 21 20:09:37 GMT 2025