1276789Sdim/* This file is distributed under the University of Illinois Open Source 2276789Sdim * License. See LICENSE.TXT for details. 3276789Sdim */ 4276789Sdim 5276789Sdim#include "DD.h" 6276789Sdim#include "../int_math.h" 7276789Sdim 8276789Sdim#define makeFinite(x) { \ 9276789Sdim (x).s.hi = crt_copysign(crt_isinf((x).s.hi) ? 1.0 : 0.0, (x).s.hi); \ 10276789Sdim (x).s.lo = 0.0; \ 11276789Sdim } 12276789Sdim 13276789Sdim#define zeroNaN(x) { \ 14276789Sdim if (crt_isnan((x).s.hi)) { \ 15276789Sdim (x).s.hi = crt_copysign(0.0, (x).s.hi); \ 16276789Sdim (x).s.lo = 0.0; \ 17276789Sdim } \ 18276789Sdim } 19276789Sdim 20276789Sdimlong double _Complex 21276789Sdim__multc3(long double a, long double b, long double c, long double d) 22276789Sdim{ 23276789Sdim long double ac = __gcc_qmul(a,c); 24276789Sdim long double bd = __gcc_qmul(b,d); 25276789Sdim long double ad = __gcc_qmul(a,d); 26276789Sdim long double bc = __gcc_qmul(b,c); 27276789Sdim 28276789Sdim DD real = { .ld = __gcc_qsub(ac,bd) }; 29276789Sdim DD imag = { .ld = __gcc_qadd(ad,bc) }; 30276789Sdim 31276789Sdim if (crt_isnan(real.s.hi) && crt_isnan(imag.s.hi)) 32276789Sdim { 33276789Sdim int recalc = 0; 34276789Sdim 35276789Sdim DD aDD = { .ld = a }; 36276789Sdim DD bDD = { .ld = b }; 37276789Sdim DD cDD = { .ld = c }; 38276789Sdim DD dDD = { .ld = d }; 39276789Sdim 40276789Sdim if (crt_isinf(aDD.s.hi) || crt_isinf(bDD.s.hi)) 41276789Sdim { 42276789Sdim makeFinite(aDD); 43276789Sdim makeFinite(bDD); 44276789Sdim zeroNaN(cDD); 45276789Sdim zeroNaN(dDD); 46276789Sdim recalc = 1; 47276789Sdim } 48276789Sdim 49276789Sdim if (crt_isinf(cDD.s.hi) || crt_isinf(dDD.s.hi)) 50276789Sdim { 51276789Sdim makeFinite(cDD); 52276789Sdim makeFinite(dDD); 53276789Sdim zeroNaN(aDD); 54276789Sdim zeroNaN(bDD); 55276789Sdim recalc = 1; 56276789Sdim } 57276789Sdim 58276789Sdim if (!recalc) 59276789Sdim { 60276789Sdim DD acDD = { .ld = ac }; 61276789Sdim DD bdDD = { .ld = bd }; 62276789Sdim DD adDD = { .ld = ad }; 63276789Sdim DD bcDD = { .ld = bc }; 64276789Sdim 65276789Sdim if (crt_isinf(acDD.s.hi) || crt_isinf(bdDD.s.hi) || 66276789Sdim crt_isinf(adDD.s.hi) || crt_isinf(bcDD.s.hi)) 67276789Sdim { 68276789Sdim zeroNaN(aDD); 69276789Sdim zeroNaN(bDD); 70276789Sdim zeroNaN(cDD); 71276789Sdim zeroNaN(dDD); 72276789Sdim recalc = 1; 73276789Sdim } 74276789Sdim } 75276789Sdim 76276789Sdim if (recalc) 77276789Sdim { 78276789Sdim real.s.hi = CRT_INFINITY * (aDD.s.hi*cDD.s.hi - bDD.s.hi*dDD.s.hi); 79276789Sdim real.s.lo = 0.0; 80276789Sdim imag.s.hi = CRT_INFINITY * (aDD.s.hi*dDD.s.hi + bDD.s.hi*cDD.s.hi); 81276789Sdim imag.s.lo = 0.0; 82276789Sdim } 83276789Sdim } 84276789Sdim 85276789Sdim long double _Complex z; 86276789Sdim __real__ z = real.ld; 87276789Sdim __imag__ z = imag.ld; 88276789Sdim 89276789Sdim return z; 90276789Sdim} 91