#include "math.h"
#include "fp_private.h"
#if (__WANT_LONG_DOUBLE_FORMAT__ - 0L == 128L)
#include "DD.h"
static const hexdouble infinity = {{0x7ff00000, 0x00000000}};
static const hexdouble tiny = {{0x00100000, 0x00000000}};
void _d3div(double a, double b, double c, double d, double e, double f,
double *high, double *mid, double *low)
{
double p, q, r, s, t, u, v, w, x, y, z;
double g, h, o, rmid, rlow;
p = 1.0/d; if (fabs(d) < tiny.d) {
r=a/d; if (fabs(r)==infinity.d) { *high = r;
*mid = 0.0e0;
*low = 0.0e0; return;
}
}
else {
q = a*p;
r = q + (a - q*d)*p; }
s = a - r*d; t = __FMUL( r, e); u = (r*e - t) + r*f; w = __FSUB( b, t ); y = __FADD( s, w );
if (fabs(b) > fabs(t)) x = b - w - t; else
x = b - (w + t);
if (fabs(s) > fabs(w))
z = s - y + w;
else
z = w - y + s;
g = c - u + x + z; if (fabs(d) < tiny.d) { o = y/d;
v = (y - __FMUL( o, d ) + g - (e*o + (f*o))) / d;
}
else {
h = y*p;
o = h + (y - h*d)*p; v = (y - o*d + g - (e*o + (f*o)))*p; }
rmid = o + v; rlow = o - rmid + v;
*high = __FADD( r, rmid );
*mid = (r - *high + rmid) + rlow;
*low = (r - *high + rmid) - *mid + rlow;
return;
}
void _d3mul(double a, double b, double c, double d, double e, double f,
double *high, double *mid, double *low)
{
double s, t, u, v, w, x, y, z, tlow, resbump;
u = d*c + e*b + f*a; v = __FMUL( d, b ); w = d*b - v;
x = __FMUL( e, a );
y = e*a - x;
s = w + y + u;
z = __FMUL( d, a);
w = d*a - z;
u = __FADD( v, x );
if (fabs(x) > fabs(v))
s = x - u + v + s;
else
s = v - u + x + s;
t = __FADD( u, w );
if (fabs(w) > fabs(u))
tlow = w - t + u + s;
else
tlow = u - t + w + s;
resbump = __FADD( t, tlow );
tlow = t - resbump + tlow;
*high = __FADD( z, resbump );
*mid = (z - *high + resbump) + tlow;
*low = (z - *high + resbump) - *mid + tlow;
return;
}
void _d3add(double a, double b, double c, double d, double e, double f,
double *high, double *mid, double *low)
{
double p, q, r, s, t, u, v, w, x, y, z;
double res, carry, z1, z2, resmid, reslow;
p = __FADD( a, d );
q = __FADD( b, e );
r = __FADD( c, f );
if (fabs(a) > fabs(d)) s = a - p + d; else s = d - p + a;
if (fabs(b) > fabs(e))
t = b - q + e; else
t = e - q + b;
if (fabs(c) > fabs(f))
u = c - r + f; else
u = f - r + c;
v = __FADD( s, q );
x = __FADD( t, r );
if (fabs(s) > fabs(q))
w = s - v + q; else
w = q - v + s;
if (fabs(t) > fabs(r))
y = t - x + r; else
y = r - x + t;
z = __FADD( x, w );
res = __FADD( p, v );
if (fabs(p) > fabs(v))
carry = p - res + v; else
carry = v - res + p;
if (fabs(x) > fabs(w))
z1 = x - z + w; else
z1 = w - z + x;
z2 = u + y + z1; resmid = __FADD( carry, z );
if (fabs(carry) > fabs(z))
reslow = carry - resmid + z; else
reslow = z - resmid + carry;
*high = __FADD( res, resmid );
*mid = (res - *high + resmid) + reslow;
*low = (res - *high + resmid) - *mid + reslow + z2;
return;
}
#endif
#if 0 // set to 1 to include code for self-contained sanity test
#include <stdio.h>
main()
{
double a1,a2,a3, b1,b2,b3, c1,c2,c3, d1,d2,d3, e1,e2,e3;
_d3div( 2, 0, 0, 3, 0, 0, &a1, &a2, &a3); _d3div( 5, 0, 0, 7, 0, 0, &b1, &b2, &b3); _d3mul(a1, a2, a3, b1, b2, b3, &c1, &c2, &c3); _d3div(21, 0, 0, 10, 0, 0, &d1, &d2, &d3); _d3mul(c1, c2, c3, d1, d2, d3, &e1, &e2, &e3);
printf("e: %24.18e %24.18e %24.18e\n", e1, e2, e3);
_d3add(a1, a2, a3, b1, b2, b3, &c1, &c2, &c3); _d3div( 8, 0, 0, 21, 0, 0, &d1, &d2, &d3); _d3add(c1, c2, c3,-d1,-d2,-d3, &e1, &e2, &e3);
printf("e: %24.18e %24.18e %24.18e\n", e1, e2, e3);
}
#endif