#include typedef struct { double hi; double lo; } DoubleDouble; inline DoubleDouble dd_new(double hi, double lo) { DoubleDouble num; num.hi = hi; num.lo = lo; return num; } unsigned int dd_get_ui(DoubleDouble num) { return num.hi; } inline DoubleDouble dd_sqrt(DoubleDouble num) { double a, b, c, d, e; d = 1 / sqrt(num.hi); e = num.hi * d; a = 0x08000001 * e; a += e - a; b = e - a; c = e * e; b = ((a * a - c) + a * b * 2) + b * b; a = num.hi - c; c = num.hi - a; c = (a + ((((num.hi - (c + a)) + (c - c)) + num.lo) - b)) * d * 0.5; a = e + c; b = e - a; return dd_new(a, (e - (b + a)) + (b + c)); } inline DoubleDouble dd_div(DoubleDouble num1, DoubleDouble num2) { double a, b, c, d, e, f, g; f = num1.hi / num2.hi; a = 0x08000001 * num2.hi; a += num2.hi - a; b = num2.hi - a; c = 0x08000001 * f; c += f - c; d = f - c; e = num2.hi * f; c = (((a * c - e) + (a * d + b * c)) + b * d) + num2.lo * f; b = num1.lo - c; d = num1.lo - b; a = num1.hi - e; e = (num1.hi - ((num1.hi - a) + a)) + b; g = a + e; e += (a - g) + ((num1.lo - (d + b)) + (d - c)); a = g + e; b = a / num2.hi; f += (e + (g - a)) / num2.hi; a = f + b; return dd_new(a, b + (f - a)); } inline DoubleDouble dd_ui_div(unsigned int num1, DoubleDouble num2) { return dd_div(dd_new(num1, 0), num2); } inline DoubleDouble dd_div_ui(DoubleDouble num1, unsigned int num2) { return dd_div(num1, dd_new(num2, 0)); } inline DoubleDouble dd_si_div(signed int num1, DoubleDouble num2) { return dd_div(dd_new(num1, 0), num2); } inline DoubleDouble dd_div_si(DoubleDouble num1, signed int num2) { return dd_div(num1, dd_new(num2, 0)); } inline DoubleDouble dd_add(DoubleDouble num1, DoubleDouble num2) { double a, b, c, d, e, f; e = num1.hi + num2.hi; d = num1.hi - e; a = num1.lo + num2.lo; f = num1.lo - a; d = ((num1.hi - (d + e)) + (d + num2.hi)) + a; b = e + d; c = ((num1.lo - (f + a)) + (f + num2.lo)) + (d + (e - b)); a = b + c; return dd_new(a, c + (b - a)); } inline DoubleDouble dd_mul(DoubleDouble num1, DoubleDouble num2) { double a, b, c, d, e; a = 0x08000001 * num1.hi; a += num1.hi - a; b = num1.hi - a; c = 0x08000001 * num2.hi; c += num2.hi - c; d = num2.hi - c; e = num1.hi * num2.hi; c = (((a * c - e) + (a * d + b * c)) + b * d) + (num1.lo * num2.hi + num1.hi * num2.lo); a = e + c; return dd_new(a, c + (e - a)); } inline DoubleDouble dd_mul2(DoubleDouble num1, DoubleDouble num2) { double a, b, c, d, e; a = 0x08000001 * num1.hi; a += num1.hi - a; b = num1.hi - a; c = 0x08000001 * num2.hi; c += num2.hi - c; d = num2.hi - c; e = num1.hi * num2.hi; c = 2*((((a * c - e) + (a * d + b * c)) + b * d) + (num1.lo * num2.hi + num1.hi * num2.lo)); a = 2*e + c; return dd_new(a, c + (2*e - a)); } inline DoubleDouble dd_mul_ui(DoubleDouble num1, unsigned int num2) { return dd_mul(num1, dd_new(num2, 0)); } inline DoubleDouble dd_mul_d(DoubleDouble num1, double num2) { return dd_mul(num1, dd_new(num2, 0)); } inline DoubleDouble dd_sub(DoubleDouble num1, DoubleDouble num2) { double a, b, c, d, e, f, g; g = num1.lo - num2.lo; f = num1.lo - g; e = num1.hi - num2.hi; d = num1.hi - e; d = ((num1.hi - (d + e)) + (d - num2.hi)) + g; b = e + d; c = (d + (e - b)) + ((num1.lo - (f + g)) + (f - num2.lo)); a = b + c; return dd_new(a, c + (b - a)); } inline DoubleDouble dd_sqr(DoubleDouble num) { double a, b, c; a = 0x08000001 * num.hi; a += num.hi - a; b = num.hi - a; c = num.hi * num.hi; b = ((((a * a - c) + a * b * 2) + b * b) + num.hi * num.lo * 2) + num.lo * num.lo; a = b + c; return dd_new(a, b + (c - a)); }