static inline double scaleDouble_1(double x, int n) { double t; {union {double d; unsigned long i;} _bitsy; _bitsy.i = (((long)n + 1023) << 52); t = _bitsy.d;}; return x*t; } static inline double scaleDouble_2(double x, int n) { double t1, t2; int n1, n2; n1 = n / 2; n2 = n - n1; {union {double d; unsigned long i;} _bitsy; _bitsy.i = (((long)n1 + 1023) << 52); t1 = _bitsy.d;}; {union {double d; unsigned long i;} _bitsy; _bitsy.i = (((long)n2 + 1023) << 52); t2 = _bitsy.d;}; return (x*t1)*t2; } double __fma(double a, double b, double sum) { double ha, ta, hb, tb, z, zz, r, s, az, asum; int ua, ub, usum; int scaled, expover, expunder, scaleexp; unsigned long u; {union {double d; unsigned long i;} _bitsy; _bitsy.d = (a); u = _bitsy.i;}; ua = (int)((u & 0x7ff0000000000000) >> 52) - 1023; {union {double d; unsigned long i;} _bitsy; _bitsy.d = (b); u = _bitsy.i;}; ub = (int)((u & 0x7ff0000000000000) >> 52) - 1023; {union {double d; unsigned long i;} _bitsy; _bitsy.d = (sum); u = _bitsy.i;}; usum = (int)((u & 0x7ff0000000000000) >> 52) - 1023; if (ua == 1023 + 1 || ub == 1023 + 1 || usum == 1023 + 1) { return a * b + sum; } else if (ua + ub > usum + 2 * 53) { return a*b; } else if (usum > ua + ub + 53) { return sum; } expover = 1023 - 2; expunder = -1022 + 53; scaleexp = 0; if (ua + ub > expover || usum > expover) { scaled = 1; scaleexp = expover / 2; a = scaleDouble_1(a, -scaleexp); b = scaleDouble_1(b, -scaleexp); sum = scaleDouble_2(sum, -2*scaleexp); } else if (ua + ub < expunder) { scaled = 1; scaleexp = expunder / 2; a = scaleDouble_1(a, -scaleexp); b = scaleDouble_1(b, -scaleexp); sum = scaleDouble_2(sum, -2*scaleexp); } else scaled = 0; ha = a; {union {double d; unsigned long i;} _bitsy; _bitsy.d = (ha); u = _bitsy.i;}; u &= 0xfffffffff8000000; {union {double d; unsigned long i;} _bitsy; _bitsy.i = (u); ha = _bitsy.d;}; ta = a - ha; hb = b; {union {double d; unsigned long i;} _bitsy; _bitsy.d = (hb); u = _bitsy.i;}; u &= 0xfffffffff8000000; {union {double d; unsigned long i;} _bitsy; _bitsy.i = (u); hb = _bitsy.d;}; tb = b - hb; z = a * b; zz = (((ha * hb - z) + ha * tb) + ta * hb) + ta * tb; {union {double d; unsigned long i;} _bitsy; _bitsy.d = (z); u = _bitsy.i;}; u &= ~0x8000000000000000; {union {double d; unsigned long i;} _bitsy; _bitsy.i = (u); az = _bitsy.d;}; {union {double d; unsigned long i;} _bitsy; _bitsy.d = (sum); u = _bitsy.i;}; u &= ~0x8000000000000000; {union {double d; unsigned long i;} _bitsy; _bitsy.i = (u); asum = _bitsy.d;}; r = z + sum; if (az > asum) s = ((z - r) + sum) + zz; else s = ((sum - r) + z) + zz; if (scaled) return scaleDouble_1(r + s, 2*scaleexp); else return r + s; }