| Daniel Dunbar | b3a6901 | 2009-06-26 16:47:03 +0000 | [diff] [blame^] | 1 | // This file is distributed under the University of Illinois Open Source |
| 2 | // License. See LICENSE.TXT for details. |
| 3 | |
| 4 | #include "DD.h" |
| 5 | #include <math.h> |
| 6 | |
| 7 | #define makeFinite(x) { \ |
| 8 | (x).hi = __builtin_copysign(isinf((x).hi) ? 1.0 : 0.0, (x).hi); \ |
| 9 | (x).lo = 0.0; \ |
| 10 | } |
| 11 | |
| 12 | long double __gcc_qadd(long double, long double); |
| 13 | long double __gcc_qsub(long double, long double); |
| 14 | long double __gcc_qmul(long double, long double); |
| 15 | long double __gcc_qdiv(long double, long double); |
| 16 | |
| 17 | long double _Complex |
| 18 | __divtc3(long double a, long double b, long double c, long double d) |
| 19 | { |
| 20 | DD cDD = { .ld = c }; |
| 21 | DD dDD = { .ld = d }; |
| 22 | |
| 23 | int ilogbw = 0; |
| 24 | const double logbw = logb(__builtin_fmax( __builtin_fabs(cDD.hi), __builtin_fabs(dDD.hi) )); |
| 25 | |
| 26 | if (isfinite(logbw)) |
| 27 | { |
| 28 | ilogbw = (int)logbw; |
| 29 | |
| 30 | cDD.hi = scalbn(cDD.hi, -ilogbw); |
| 31 | cDD.lo = scalbn(cDD.lo, -ilogbw); |
| 32 | dDD.hi = scalbn(dDD.hi, -ilogbw); |
| 33 | dDD.lo = scalbn(dDD.lo, -ilogbw); |
| 34 | } |
| 35 | |
| 36 | const long double denom = __gcc_qadd(__gcc_qmul(cDD.ld, cDD.ld), __gcc_qmul(dDD.ld, dDD.ld)); |
| 37 | const long double realNumerator = __gcc_qadd(__gcc_qmul(a,cDD.ld), __gcc_qmul(b,dDD.ld)); |
| 38 | const long double imagNumerator = __gcc_qsub(__gcc_qmul(b,cDD.ld), __gcc_qmul(a,dDD.ld)); |
| 39 | |
| 40 | DD real = { .ld = __gcc_qdiv(realNumerator, denom) }; |
| 41 | DD imag = { .ld = __gcc_qdiv(imagNumerator, denom) }; |
| 42 | |
| 43 | real.hi = scalbn(real.hi, -ilogbw); |
| 44 | real.lo = scalbn(real.lo, -ilogbw); |
| 45 | imag.hi = scalbn(imag.hi, -ilogbw); |
| 46 | imag.lo = scalbn(imag.lo, -ilogbw); |
| 47 | |
| 48 | if (isnan(real.hi) && isnan(imag.hi)) |
| 49 | { |
| 50 | DD aDD = { .ld = a }; |
| 51 | DD bDD = { .ld = b }; |
| 52 | DD rDD = { .ld = denom }; |
| 53 | |
| 54 | if ((rDD.hi == 0.0) && (!isnan(aDD.hi) || !isnan(bDD.hi))) |
| 55 | { |
| 56 | real.hi = __builtin_copysign(INFINITY,cDD.hi) * aDD.hi; |
| 57 | real.lo = 0.0; |
| 58 | imag.hi = __builtin_copysign(INFINITY,cDD.hi) * bDD.hi; |
| 59 | imag.lo = 0.0; |
| 60 | } |
| 61 | |
| 62 | else if ((isinf(aDD.hi) || isinf(bDD.hi)) && isfinite(cDD.hi) && isfinite(dDD.hi)) |
| 63 | { |
| 64 | makeFinite(aDD); |
| 65 | makeFinite(bDD); |
| 66 | real.hi = INFINITY * (aDD.hi*cDD.hi + bDD.hi*dDD.hi); |
| 67 | real.lo = 0.0; |
| 68 | imag.hi = INFINITY * (bDD.hi*cDD.hi - aDD.hi*dDD.hi); |
| 69 | imag.lo = 0.0; |
| 70 | } |
| 71 | |
| 72 | else if ((isinf(cDD.hi) || isinf(dDD.hi)) && isfinite(aDD.hi) && isfinite(bDD.hi)) |
| 73 | { |
| 74 | makeFinite(cDD); |
| 75 | makeFinite(dDD); |
| 76 | real.hi = __builtin_copysign(0.0,(aDD.hi*cDD.hi + bDD.hi*dDD.hi)); |
| 77 | real.lo = 0.0; |
| 78 | imag.hi = __builtin_copysign(0.0,(bDD.hi*cDD.hi - aDD.hi*dDD.hi)); |
| 79 | imag.lo = 0.0; |
| 80 | } |
| 81 | } |
| 82 | |
| 83 | long double _Complex z; |
| 84 | __real__ z = real.ld; |
| 85 | __imag__ z = imag.ld; |
| 86 | |
| 87 | return z; |
| 88 | } |