--- a/QDouble.st Fri Nov 29 15:00:42 2019 +0100
+++ b/QDouble.st Fri Nov 29 16:30:23 2019 +0100
@@ -367,7 +367,7 @@
// three_sum2 : d + e = a + b + c, de are nonoverlapping, Bailey
// two_prod : s + e = a * b, s = fl(a * b), e = err(a * b), Verkamp and Dekker
// sqr : s + e = a^2, s = fl(a * a), e = err(a * a)
-// renormalize : renormalization algorithm
+// renorm : renormalization algorithm
// qd_add_s : qd + double
// qd_add_qd : qd + qd (sloppy add)
// s_sub_qd : double - qd
@@ -375,8 +375,9 @@
// s_mul_qd : double * qd
// qd_mul_qd : qd * qd
// qd_div_qd : qd / qd
-// qd_sqr : qd ^ 2
-// qdsqrt_sci : square root (scalar)
+// qd_sqr : qd ^ 2
+// qd_sqrt : square root (scalar)
+// qd_pow : qd ^ n (n integer)
static INLINE void
fast_two_sum(double *s, double *e, double a, double b)
@@ -665,17 +666,18 @@
//
//--------------------------------------------------------------------------------------------
static INLINE void
-qd_add_s(double *c0, double *c1, double *c2, double *c3, double a0, double a1, double a2, double a3, double b)
+qd_add_s(double *o0, double *o1, double *o2, double *o3, double a0, double a1, double a2, double a3, double b)
{
double e,x,y,w,z;
- c0[0] = 0.0; c1[0] = 0.0; c2[0] = 0.0; c3[0] = 0.0;
-
- two_sum(&x, &y, a0, b); c0[0] = x; e = y;
- two_sum(&x, &y, a1, e); c1[0] = x; e = y;
- two_sum(&x, &y, a2, e); c2[0] = x; e = y;
- two_sum(&x, &y, a3, e); c3[0] = x; e = y;
- renorm(&x, &y, &w, &z, c0[0], c1[0], c2[0], c3[0], e);
- c0[0] = x; c1[0] = y; c2[0] = w; c3[0] = z;
+ double c0, c1, c2, c3;
+ c0 = 0.0; c1 = 0.0; c2 = 0.0; c3 = 0.0;
+
+ two_sum(&x, &y, a0, b); c0 = x; e = y;
+ two_sum(&x, &y, a1, e); c1 = x; e = y;
+ two_sum(&x, &y, a2, e); c2 = x; e = y;
+ two_sum(&x, &y, a3, e); c3 = x; e = y;
+ renorm(&x, &y, &w, &z, c0, c1, c2, c3, e);
+ o0[0] = x; o1[0] = y; o2[0] = w; o3[0] = z;
}
//--------------------------------------------------------------------------------------------
@@ -684,21 +686,22 @@
//
//--------------------------------------------------------------------------------------------
static INLINE void
-qd_add_dd(double *c0, double *c1, double *c2, double *c3, double a0, double a1, double a2, double a3, double b0, double b1)
+qd_add_dd(double *o0, double *o1, double *o2, double *o3, double a0, double a1, double a2, double a3, double b0, double b1)
{
double e1,e2,x,y,w,z;
- c0[0] = 0.0; c1[0] = 0.0; c2[0] = 0.0; c3[0] = 0.0;
-
- two_sum(&x, &y, a0, b0); c0[0] = x; e1 = y;
- two_sum(&x, &y, a1, b1); c1[0] = x; e2 = y;
- two_sum(&x, &y, c1[0], e1); c1[0] = x; e1 = y;
- two_sum(&x, &y, a2, e2); c2[0] = x; e2 = y;
- two_sum(&x, &y, c2[0], e1); c2[0] = x; e1 = y;
- two_sum(&x, &y, e1, e2); e1 = x; e2 = y;
- two_sum(&x, &y, a3, e1); c3[0] = x; e1 = y;
+ double c0, c1, c2, c3;
+ c0 = 0.0; c1 = 0.0; c2 = 0.0; c3 = 0.0;
+
+ two_sum(&x, &y, a0, b0); c0 = x; e1 = y;
+ two_sum(&x, &y, a1, b1); c1 = x; e2 = y;
+ two_sum(&x, &y, c1, e1); c1 = x; e1 = y;
+ two_sum(&x, &y, a2, e2); c2 = x; e2 = y;
+ two_sum(&x, &y, c2, e1); c2 = x; e1 = y;
+ two_sum(&x, &y, e1, e2); e1 = x; e2 = y;
+ two_sum(&x, &y, a3, e1); c3 = x; e1 = y;
e1 = e1 + e2;
- renorm(&x, &y, &w, &z, c0[0], c1[0], c2[0], c3[0], e1);
- c0[0] = x; c1[0] = y; c2[0] = w; c3[0] = z;
+ renorm(&x, &y, &w, &z, c0, c1, c2, c3, e1);
+ o0[0] = x; o1[0] = y; o2[0] = w; o3[0] = z;
return;
}
@@ -716,6 +719,7 @@
two_sum(&x, &y, a0, b0); c0[0] = x; e1 = y;
two_sum(&x, &y, a1, b1); c1[0] = x; e2 = y;
two_sum(&x, &y, c1[0], e1); c1[0] = x; e1 = y;
+
two_sum(&x, &y, a2, b2); c2[0] = x; e3 = y;
three_sum(&x, &y, &z, c2[0], e2, e1); c2[0] = x; e1 = y; e2 = z;
two_sum(&x, &y, a3, b3); c3[0] = x; e4 = y;
@@ -731,31 +735,31 @@
//
//--------------------------------------------------------------------------------------------
static INLINE void
-s_sub_qd(double *c0, double *c1, double *c2, double *c3, double a, double b0, double b1, double b2, double b3)
+s_sub_qd(double *o0, double *o1, double *o2, double *o3, double a, double b0, double b1, double b2, double b3)
{
- double e,x,y,w,z;
- e=0.0;
- c0[0] = 0.0; c1[0] = 0.0; c2[0] = 0.0; c3[0] = 0.0;
- b0=-b0; b1=-b1; b2=-b2; b3=-b3;
- two_sum(&x, &y, a, b0);
- c0[0] = x;
- e = y;
- two_sum(&x, &y, b1, e);
- c1[0] = x;
- e = y;
- two_sum(&x, &y, b2, e);
- c2[0] = x;
- e = y;
- two_sum(&x, &y, b3, e);
- c3[0] = x;
- e = y;
- renorm(&x, &y, &w, &z, c0[0], c1[0], c2[0], c3[0], e);
- c0[0] = x;
- c1[0] = y;
- c2[0] = w;
- c3[0] = z;
-
- return;
+ double e,x,y,w,z;
+ double c0, c1, c2, c3;
+
+ e=0.0;
+ c0 = 0.0; c1 = 0.0; c2 = 0.0; c3 = 0.0;
+ b0=-b0; b1=-b1; b2=-b2; b3=-b3;
+ two_sum(&x, &y, a, b0);
+ c0 = x;
+ e = y;
+ two_sum(&x, &y, b1, e);
+ c1 = x;
+ e = y;
+ two_sum(&x, &y, b2, e);
+ c2 = x;
+ e = y;
+ two_sum(&x, &y, b3, e);
+ c3 = x;
+ e = y;
+ renorm(&x, &y, &w, &z, c0, c1, c2, c3, e);
+ o0[0] = x;
+ o1[0] = y;
+ o2[0] = w;
+ o3[0] = z;
}
//--------------------------------------------------------------------------------------------
@@ -788,22 +792,24 @@
//
//--------------------------------------------------------------------------------------------
static INLINE void
-s_mul_qd(double *c0, double *c1, double *c2, double *c3, double b, double a0, double a1, double a2, double a3)
+s_mul_qd(double *o0, double *o1, double *o2, double *o3, double b, double a0, double a1, double a2, double a3)
{
double e0,e1,e2,x,y,w,z;
- c0[0] = 0.0; c1[0] = 0.0; c2[0] = 0.0; c3[0] = 0.0;
-
- two_prod(&x, &y, a0, b); c0[0] = x; e0 = y;
- two_prod(&x, &y, a1, b); c1[0] = x; e1 = y;
- two_sum(&x, &y, c1[0], e0); c1[0] = x; e0 = y;
- two_prod(&x, &y, a2, b); c2[0] = x; e2 = y;
- three_sum(&x, &y, &z, c2[0], e1, e0); c2[0] = x; e0 = y; e1 = z;
- c3[0] = a3*b;
- three_sum2(&x, &y, c3[0], e2, e0); c3[0] = x; e0 = y;
+ double c0, c1, c2, c3;
+
+ c0 = 0.0; c1 = 0.0; c2 = 0.0; c3 = 0.0;
+
+ two_prod(&x, &y, a0, b); c0 = x; e0 = y;
+ two_prod(&x, &y, a1, b); c1 = x; e1 = y;
+ two_sum(&x, &y, c1, e0); c1 = x; e0 = y;
+ two_prod(&x, &y, a2, b); c2 = x; e2 = y;
+ three_sum(&x, &y, &z, c2, e1, e0); c2 = x; e0 = y; e1 = z;
+ c3 = a3*b;
+ three_sum2(&x, &y, c3, e2, e0); c3 = x; e0 = y;
e0 = e0 + e1;
- renorm(&x, &y, &w, &z, c0[0], c1[0], c2[0], c3[0], e0);
- c0[0] = x; c1[0] = y; c2[0] = w; c3[0] = z;
+ renorm(&x, &y, &w, &z, c0, c1, c2, c3, e0);
+ o0[0] = x; o1[0] = y; o2[0] = w; o3[0] = z;
}
//--------------------------------------------------------------------------------------------
@@ -1014,7 +1020,7 @@
}
static INLINE void
-qdsqrt(double *c0, double *c1, double *c2, double *c3, double a0, double a1, double a2, double a3)
+qd_sqrt(double *c0, double *c1, double *c2, double *c3, double a0, double a1, double a2, double a3)
{
double h0,h1,h2,h3,x0,x1,x2,x3,p,q,r,s;
@@ -1051,7 +1057,7 @@
}
static void
-qdpow(double *c0, double *c1, double *c2, double *c3, double a0, double a1, double a2, double a3, int p)
+qd_pow(double *c0, double *c1, double *c2, double *c3, double a0, double a1, double a2, double a3, int p)
{
double r0,r1,r2,r3,x,y,w,z;
int abs_p;
@@ -1663,54 +1669,6 @@
return;
}
-static void
-sin_taylor_qd(double *s0, double *s1, double *s2, double *s3, double x0, double x1, double x2, double x3)
-{
- double eps = 1.21543267145725e-63; // = 2^-209
- double thresh = 0.5*fabs(x0)*eps;
- double fact[15][4] = {
- {1.66666666666666657e-01, 9.25185853854297066e-18, 5.13581318503262866e-34, 2.85094902409834186e-50},
- {4.16666666666666644e-02, 2.31296463463574266e-18, 1.28395329625815716e-34, 7.12737256024585466e-51},
- {8.33333333333333322e-03, 1.15648231731787138e-19, 1.60494162032269652e-36, 2.22730392507682967e-53},
- {1.38888888888888894e-03, -5.30054395437357706e-20, -1.73868675534958776e-36, -1.63335621172300840e-52},
- {1.98412698412698413e-04, 1.72095582934207053e-22, 1.49269123913941271e-40, 1.29470326746002471e-58},
- {2.48015873015873016e-05, 2.15119478667758816e-23, 1.86586404892426588e-41, 1.61837908432503088e-59},
- {2.75573192239858925e-06, -1.85839327404647208e-22, 8.49175460488199287e-39, -5.72661640789429621e-55},
- {2.75573192239858883e-07, 2.37677146222502973e-23, -3.26318890334088294e-40, 1.61435111860404415e-56},
- {2.50521083854417202e-08, -1.44881407093591197e-24, 2.04267351467144546e-41, -8.49632672007163175e-58},
- {2.08767569878681002e-09, -1.20734505911325997e-25, 1.70222792889287100e-42, 1.41609532150396700e-58},
- {1.60590438368216133e-10, 1.25852945887520981e-26, -5.31334602762985031e-43, 3.54021472597605528e-59},
- {1.14707455977297245e-11, 2.06555127528307454e-28, 6.88907923246664603e-45, 5.72920002655109095e-61},
- {7.64716373181981641e-13, 7.03872877733453001e-30, -7.82753927716258345e-48, 1.92138649443790242e-64},
- {4.77947733238738525e-14, 4.39920548583408126e-31, -4.89221204822661465e-49, 1.20086655902368901e-65},
- {2.81145725434552060e-15, 1.65088427308614326e-31, -2.87777179307447918e-50, 4.27110689256293549e-67}
- };
- double y0,y1,y2,y3,r0,r1,r2,r3,t0,t1,t2,t3;
- int i;
-
- if(x0==0.0) {
- s0[0]=0.0; s1[0]=0.0; s2[0]=0.0; s3[0]=0.0;
- return;
- }
-
- i=0;
- qd_mul_qd(&y0,&y1,&y2,&y3,x0,x1,x2,x3,x0,x1,x2,x3);
- y0 = -y0; y1 = -y1; y2 = -y2; y3 = -y3;
- s0[0]=x0; s1[0]=x1; s2[0]=x2; s3[0]=x3;
- r0=x0; r1=x1; r2=x2; r3=x3;
-
- qd_mul_qd(&r0,&r1,&r2,&r3,r0,r1,r2,r3,y0,y1,y2,y3);
- qd_mul_qd(&t0,&t1,&t2,&t3,r0,r1,r2,r3,fact[i][0],fact[i][1],fact[i][2],fact[i][3]);
- qd_add_qd(&s0[0],&s1[0],&s2[0],&s3[0],s0[0],s1[0],s2[0],s3[0],t0,t1,t2,t3);
- i=i+2;
- while ((i<=15)||(fabs(t0)>thresh)) {
- qd_mul_qd(&r0,&r1,&r2,&r3,r0,r1,r2,r3,y0,y1,y2,y3);
- qd_mul_qd(&t0,&t1,&t2,&t3,r0,r1,r2,r3,fact[i][0],fact[i][1],fact[i][2],fact[i][3]);
- qd_add_qd(&s0[0],&s1[0],&s2[0],&s3[0],s0[0],s1[0],s2[0],s3[0],t0,t1,t2,t3);
- i=i+2;
- }
-}
-
static double inv_fact[15][4] = {
{1.66666666666666657e-01, 9.25185853854297066e-18, 5.13581318503262866e-34, 2.85094902409834186e-50},
{4.16666666666666644e-02, 2.31296463463574266e-18, 1.28395329625815716e-34, 7.12737256024585466e-51},
@@ -1730,6 +1688,37 @@
};
static void
+sin_taylor_qd(double *s0, double *s1, double *s2, double *s3, double x0, double x1, double x2, double x3)
+{
+ double eps = 1.21543267145725e-63; // = 2^-209
+ double thresh = 0.5*fabs(x0)*eps;
+ double y0,y1,y2,y3,r0,r1,r2,r3,t0,t1,t2,t3;
+ int i;
+
+ if(x0==0.0) {
+ s0[0]=0.0; s1[0]=0.0; s2[0]=0.0; s3[0]=0.0;
+ return;
+ }
+
+ i=0;
+ qd_mul_qd(&y0,&y1,&y2,&y3,x0,x1,x2,x3,x0,x1,x2,x3);
+ y0 = -y0; y1 = -y1; y2 = -y2; y3 = -y3;
+ s0[0]=x0; s1[0]=x1; s2[0]=x2; s3[0]=x3;
+ r0=x0; r1=x1; r2=x2; r3=x3;
+
+ qd_mul_qd(&r0,&r1,&r2,&r3,r0,r1,r2,r3,y0,y1,y2,y3);
+ qd_mul_qd(&t0,&t1,&t2,&t3,r0,r1,r2,r3,inv_fact[i][0],inv_fact[i][1],inv_fact[i][2],inv_fact[i][3]);
+ qd_add_qd(&s0[0],&s1[0],&s2[0],&s3[0],s0[0],s1[0],s2[0],s3[0],t0,t1,t2,t3);
+ i=i+2;
+ while ((i<=15)||(fabs(t0)>thresh)) {
+ qd_mul_qd(&r0,&r1,&r2,&r3,r0,r1,r2,r3,y0,y1,y2,y3);
+ qd_mul_qd(&t0,&t1,&t2,&t3,r0,r1,r2,r3,inv_fact[i][0],inv_fact[i][1],inv_fact[i][2],inv_fact[i][3]);
+ qd_add_qd(&s0[0],&s1[0],&s2[0],&s3[0],s0[0],s1[0],s2[0],s3[0],t0,t1,t2,t3);
+ i=i+2;
+ }
+}
+
+static void
cos_taylor_qd(double *c0, double *c1, double *c2, double *c3, double x0, double x1, double x2, double x3)
{
double eps = 1.21543267145725e-63; // = 2^-209
@@ -1794,7 +1783,7 @@
}
qd_mul_qd(&p0,&p1,&p2,&p3,s0[0],s1[0],s2[0],s3[0],s0[0],s1[0],s2[0],s3[0]); // Modified,2012/01/16 Saito
s_sub_qd(&q0,&q1,&q2,&q3,1.0,p0,p1,p2,p3);
- qdsqrt(&c0[0],&c1[0],&c2[0],&c3[0],q0,q1,q2,q3);
+ qd_sqrt(&c0[0],&c1[0],&c2[0],&c3[0],q0,q1,q2,q3);
}
//--------------------------------------------------------------------------------------------
@@ -1809,7 +1798,7 @@
//
//--------------------------------------------------------------------------------------------
static void
-qdsin(double *s0, double *s1, double *s2, double *s3, double *a0, double *a1, double *a2, double *a3)
+qd_sin(double *s0, double *s1, double *s2, double *s3, double *a0, double *a1, double *a2, double *a3)
{
double p0,p1,p2,p3, q0,q1,q2,q3, z0,z1,z2,z3, r0,r1,r2,r3, j, t0,t1,t2,t3, k,abs_k;
double _2pi[4] = {6.283185307179586232e+00,2.449293598294706414e-16,-5.989539619436679332e-33,2.224908441726730563e-49}; // 2*pi
@@ -1940,7 +1929,7 @@
//
//--------------------------------------------------------------------------------------------
static void
-qdcos(double *c0, double *c1, double *c2, double *c3, double *a0, double *a1, double *a2, double *a3)
+qd_cos(double *c0, double *c1, double *c2, double *c3, double *a0, double *a1, double *a2, double *a3)
{
double p0,p1,p2,p3, q0,q1,q2,q3, z0,z1,z2,z3, r0,r1,r2,r3, j, t0,t1,t2,t3, k,abs_k;
double _2pi[4] = {6.283185307179586232e+00,2.449293598294706414e-16,-5.989539619436679332e-33,2.224908441726730563e-49}; // 2*pi
@@ -2245,7 +2234,7 @@
//
//--------------------------------------------------------------------------------------------
static void
-qdtan(double *t0, double *t1, double *t2, double *t3, double *a0, double *a1, double *a2, double *a3)
+qd_tan(double *t0, double *t1, double *t2, double *t3, double *a0, double *a1, double *a2, double *a3)
{
double sin0,sin1,sin2,sin3,cos0,cos1,cos2,cos3;
sincos_qd(&sin0,&sin1,&sin2,&sin3,&cos0,&cos1,&cos2,&cos3,a0[0],a1[0],a2[0],a3[0]);
@@ -2264,7 +2253,7 @@
//
//--------------------------------------------------------------------------------------------
static void
-qdexp(double *e0, double *e1, double *e2, double *e3, double *x0, double *x1, double *x2, double *x3) {
+qd_exp(double *e0, double *e1, double *e2, double *e3, double *x0, double *x1, double *x2, double *x3) {
double k = ldexp(1.0,16);
double inv_k = 1.0 / k;
@@ -2353,7 +2342,7 @@
qd_add_s(&s0, &s1, &s2, &s3, s0, s1, s2, s3, 1.0);
for(i=0; i<m; i++) {
- t=t*2;
+ t=t*2;
}
e0[0] = s0*t; e1[0] = s1*t; e2[0] = s2*t; e3[0] = s3*t;
@@ -2838,25 +2827,25 @@
if (__isSmallInteger(anInteger)) {
INT iVal = __smallIntegerVal(anInteger);
- double *__d__;
+ double *d;
__qNew(newQD, sizeof(struct __qDoubleStruct));
__stx_setClass(newQD, QDouble);
- __d__ = __QDoubleInstPtr(newQD)->d_qDoubleValue;
- __d__[1] = 0.0;
- __d__[2] = 0.0;
- __d__[3] = 0.0;
+ d = __QDoubleInstPtr(newQD)->d_qDoubleValue;
+ d[1] = 0.0;
+ d[2] = 0.0;
+ d[3] = 0.0;
// need more than 52bits?
if ((sizeof(INT) > 52)
&& ((iVal > 0xFFFFFFFFFFFFF) || (iVal < -0xFFFFFFFFFFFFF))) {
- __d__[0] = (double)(iVal & ~0xFFFFFFFF);
- __d__[1] = (double)(iVal & 0xFFFFFFFF);
- renorm(&(__d__[0]), &(__d__[1]), &(__d__[2]), &(__d__[3]), __d__[0], __d__[0], __d__[0], __d__[0], 0.0);
+ d[0] = (double)(iVal & ~0xFFFFFFFF);
+ d[1] = (double)(iVal & 0xFFFFFFFF);
+ renorm(&(d[0]), &(d[1]), &(d[2]), &(d[3]), d[0], d[1], 0.0, 0.0, 0.0);
// renorm4(&(a[0]), &(a[1]), &(a[2]), &(a[3]));
} else {
- __d__[0] = (double)iVal;
+ d[0] = (double)iVal;
}
RETURN (newQD);
}
@@ -4021,7 +4010,7 @@
int savedCV;
fpu_fix_start(&savedCV);
- qdcos(&q0, &q1, &q2, &q3, &a[0], &a[1], &a[2], &a[3]);
+ qd_cos(&q0, &q1, &q2, &q3, &a[0], &a[1], &a[2], &a[3]);
fpu_fix_end(&savedCV);
__qNew_qdReal(newQD, q0, q1, q2, q3);
RETURN( newQD );
@@ -4042,7 +4031,7 @@
int savedCV;
fpu_fix_start(&savedCV);
- qdexp(&q0, &q1, &q2, &q3, &a[0], &a[1], &a[2], &a[3]);
+ qd_exp(&q0, &q1, &q2, &q3, &a[0], &a[1], &a[2], &a[3]);
fpu_fix_end(&savedCV);
__qNew_qdReal(newQD, q0, q1, q2, q3);
RETURN( newQD );
@@ -4179,7 +4168,7 @@
int savedCV;
fpu_fix_start(&savedCV);
- qdpow(&q0, &q1, &q2, &q3, a[0], a[1], a[2], a[3], __intVal(n));
+ qd_pow(&q0, &q1, &q2, &q3, a[0], a[1], a[2], a[3], __intVal(n));
fpu_fix_end(&savedCV);
__qNew_qdReal(newQD, q0, q1, q2, q3);
RETURN( newQD );
@@ -4204,7 +4193,7 @@
int savedCV;
fpu_fix_start(&savedCV);
- qdsin(&q0, &q1, &q2, &q3, &a[0], &a[1], &a[2], &a[3]);
+ qd_sin(&q0, &q1, &q2, &q3, &a[0], &a[1], &a[2], &a[3]);
fpu_fix_end(&savedCV);
__qNew_qdReal(newQD, q0, q1, q2, q3);
RETURN( newQD );
@@ -4226,7 +4215,7 @@
int savedCV;
fpu_fix_start(&savedCV);
- qdsqrt(&q0, &q1, &q2, &q3, a[0], a[1], a[2], a[3]);
+ qd_sqrt(&q0, &q1, &q2, &q3, a[0], a[1], a[2], a[3]);
fpu_fix_end(&savedCV);
__qNew_qdReal(newQD, q0, q1, q2, q3);
RETURN( newQD );
@@ -4273,7 +4262,7 @@
int savedCV;
fpu_fix_start(&savedCV);
- qdtan(&q0, &q1, &q2, &q3, &a[0], &a[1], &a[2], &a[3]);
+ qd_tan(&q0, &q1, &q2, &q3, &a[0], &a[1], &a[2], &a[3]);
fpu_fix_end(&savedCV);
__qNew_qdReal(newQD, q0, q1, q2, q3);
RETURN( newQD );
@@ -4757,7 +4746,9 @@
x0 -= 1.0;
}
}
- m_renorm4(x0, x1, x2, x3);
+ renorm(&x0, &x1, &x2, &x3, x0, x1, x2, x3, 0.0);
+ // m_renorm4(x0, x1, x2, x3);
+
__qNew_qdReal(newQD, x0, x1, x2, x3);
RETURN( newQD );
%}.