--- a/QDouble.st Wed Nov 27 11:16:57 2019 +0100
+++ b/QDouble.st Fri Nov 29 14:21:34 2019 +0100
@@ -116,47 +116,23 @@
#endif
-
-struct qd_real {
- double x[4]; /* The Components. */
-};
-
-struct __quadDoubleStruct {
- STX_OBJ_HEADER
-#ifdef __NEED_DOUBLE_ALIGN
- __FILLTYPE_DOUBLE f_filler;
-#endif
- double d_quadDoubleValue[4];
-};
-
-#define __QuadDoubleInstPtr(obj) ((struct __quadDoubleStruct *)(__objPtr(obj)))
-
-#ifndef __isQuadDouble
-# define __isQuadDouble(o) \
- (__Class(o) == @global(QuadDouble))
-#endif
-
-#ifndef __qIsQuadDouble
-# define __qIsQuadDouble(o) \
- (__qClass(o) == @global(QuadDouble))
-#endif
-
#define __qNew_qdReal(newQD, d0,d1,d2,d3) { \
- __qNew(newQD, sizeof(struct __quadDoubleStruct)); \
- __stx_setClass(newQD, QDouble); \
- __QuadDoubleInstPtr(newQD)->d_quadDoubleValue[0] = d0; \
- __QuadDoubleInstPtr(newQD)->d_quadDoubleValue[1] = d1; \
- __QuadDoubleInstPtr(newQD)->d_quadDoubleValue[2] = d2; \
- __QuadDoubleInstPtr(newQD)->d_quadDoubleValue[3] = d3; \
+ double* __d__; \
+ __qNew(newQD, sizeof(struct __qDoubleStruct)); \
+ __stx_setClass(newQD, QDouble); \
+ __d__ = __QDoubleInstPtr(newQD)->d_qDoubleValue; \
+ __d__[0] = d0; \
+ __d__[1] = d1; \
+ __d__[2] = d2; \
+ __d__[3] = d3; \
}
+
// sigh: not all compilers (borland) support inline functions;
// therefore we have to use macros...
// sigh2: c-macros are unhygienic - to avoid catching/hiding variable bindings,
// use different names in each macros (i.e. a_xxx)
-// qd_real(c0, c1, c2, c3);
-
#define _QD_SPLITTER 134217729.0 // = 2^27 + 1
#define _QD_SPLIT_THRESH 6.69692879491417e+299 // = 2^996
@@ -378,11 +354,2017 @@
!QDouble primitiveFunctions!
%{
+#ifdef __BORLANDC__#
+# define INLINE /* */
+#else
+# define INLINE inline
+#endif
+
+// routines
+// fast_two_sum : s + e = a + b, s = fl(a + b), e = err(a + b), assumption |a|>|b|, Dekker('71)
+// two_sum : s + e = a + b, s = fl(a + b), e = err(a + b), Knuth('69)
+// three_sum : d + e + f = a + b + c, def are nonoverlapping, Bailey
+// 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
+// qd_add_s : qd + double
+// qd_add_qd : qd + qd (sloppy add)
+// s_sub_qd : double - qd
+// qd_sub_qd : qd - qd
+// s_mul_qd : double * qd
+// qd_mul_qd : qd * qd
+// qd_div_qd : qd / qd
+// qd_sqr : qd ^ 2
+// qdsqrt_sci : square root (scalar)
+
+static INLINE void
+fast_two_sum(double *s, double *e, double a, double b)
+{
+ double v;
+ s[0] = 0.0; e[0] = 0.0;
+
+ s[0] = a + b;
+ v = s[0] - a;
+ e[0] = b - v;
+}
+
+static INLINE void
+two_sum(double *s, double *e, double a, double b)
+{
+ double v;
+ s[0] = 0.0; e[0] = 0.0;
+
+ s[0] = a + b;
+ v = s[0] - a;
+ e[0] = (a - (s[0] - v)) + (b - v);
+}
+
+static INLINE void
+three_sum(double *d, double *e, double *f, double a, double b, double c)
+{
+ double t1,t2,t3,v;
+
+ d[0] = 0.0; e[0] = 0.0; f[0] = 0.0;
+
+ t1= a + b;
+ v = t1 - a;
+ t2= (a - (t1 - v))+(b - v);
+
+ d[0] = t1 + c;
+ v = d[0] - t1;
+ t3= (t1 - (d[0] - v))+(c - v);
+
+ e[0] = t2 + t3;
+ v = e[0] - t2;
+ f[0] = (t2 - (e[0] - v))+(t3 - v);
+}
+
+static INLINE void
+three_sum2(double *d, double *e, double a, double b, double c)
+{
+ double t1,t2,t3,v;
+ d[0] = 0.0; e[0] = 0.0;
+
+ t1= a + b;
+ v = t1 - a;
+ t2= (a - (t1 - v))+(b - v);
+
+ d[0] = t1 + c;
+ v = d[0] - t1;
+ t3= (t1 - (d[0] - v))+(c - v);
+
+ e[0] = t2 + t3;
+}
+
+static INLINE void
+two_prod(double *p, double *e, double a, double b)
+{
+ double t,ah,al,bh,bl;
+
+ p[0] = a * b;
+
+ t = 134217729 * a; // splitter: 2^27 + 1
+ ah = t -(t - a);
+ al = a - ah;
+ t = 134217729 * b; // splitter: 2^27 + 1
+ bh = t -(t - b);
+ bl = b - bh;
+
+ e[0] = ((ah*bh - p[0]) + ah*bl + al*bh) + al*bl;
+}
+
+static INLINE void
+sqr(double *p, double *e, double a)
+{
+ double t,ah,al;
+
+ p[0] = a * a;
+ t = 134217729 * a; // splitter: 2^27 + 1
+ ah = t -(t - a);
+ al = a - ah;
+ e[0] = ((ah*ah - p[0]) + (ah*al)*2.0) + al*al;
+}
+
+static INLINE void
+renorm(double *b0, double *b1, double *b2, double *b3, double a0, double a1, double a2, double a3, double a4)
+{
+ double t0,t1,t2,t3,t4,s,ss;
+ s = 0.0; ss = 0.0;
+ t0 = 0.0;t1 = 0.0;t2 = 0.0;t3 = 0.0;t4 = 0.0;
+
+// fast_two_sum(&x, &y, a3, a4);
+// s = x;
+// t4 = y;
+// fast_two_sum(&x, &y, a2, s);
+// s = x;
+// t3 = y;
+// fast_two_sum(&x, &y, a1, s);
+// s = x;
+// t2 = y;
+// fast_two_sum(&x, &y, a0, s);
+// t0 = x;
+// t1 = y;
+// if(t1 != 0.0) {
+// fast_two_sum(&x, &y, t1, t2);
+// t1 = x;
+// t2 = y;
+// if(t2 != 0.0) {
+// fast_two_sum(&x, &y,t2, t3);
+// t2 = x;
+// t3 = y;
+// if(t3 != 0.0) {
+// t3 += t4;
+// } else {
+// t2 += t4;
+// }
+// } else {
+// fast_two_sum(&x, &y, t1, t3);
+// t1 = x;
+// t2 = y;
+// if(t2 != 0.0) {
+// fast_two_sum(&x, &y, t2, t4);
+// t2 = x;
+// t3 = y;
+// } else {
+// fast_two_sum(&x, &y, t1, t4);
+// t1 = x;
+// t2 = y;
+// }
+// }
+// } else {
+// fast_two_sum(&x, &y, t0, t2);
+// t0 = x;
+// t1 = y;
+// if(t1 != 0.0) {
+// fast_two_sum(&x, &y, t1, t3);
+// t1 = x;
+// t2 = y;
+// if(t2 != 0.0) {
+// fast_two_sum(&x, &y, t2, t4);
+// t2 = x;
+// t3 = y;
+// } else {
+// fast_two_sum(&x, &y, t1, t4);
+// t1 = x;
+// t2 = y;
+// }
+// } else {
+// fast_two_sum(&x, &y, t0, t3);
+// t0 = x;
+// t1 = y;
+// if(t1 != 0.0) {
+// fast_two_sum(&x, &y, t1, t4);
+// t1 = x;
+// t2 = y;
+// } else {
+// fast_two_sum(&x, &y, t0, t4);
+// t0 = x;
+// t1 = y;
+// }
+// }
+// }
+
+
+ //[s,t4] = fast_two_sum(a4,a5);
+ s = a3 + a4;
+ t3 = a4 - (s - a3);
+ //[ss,t3] = fast_two_sum(a3,s);
+ ss = a2 + s;
+ t2 = s - (ss - a2);
+ //[s,t2] = fast_two_sum(a2,ss);
+ s = a1 + ss;
+ t1 = ss - (s - a1);
+ //[b1,t1] = fast_two_sum(a1,s);
+ b0[0] = a0 + s;
+ t0 = s - (b0[0] - a0);
+ //[s,t3] = fast_two_sum(t3,t4);
+ s = t2 + t3;
+ t2 = t3 - (s - t2);
+ //[ss,t2] = fast_two_sum(t2,s);
+ ss = t1 + s;
+ t1 = s - (ss - t1);
+ //[b2,t1] = fast_two_sum(t1,ss);
+ b1[0] = t0 + ss;
+ t0 = ss - (b1[0] - t0);
+ //[s,t2] = fast_two_sum(t2,t3);
+ s = t1 + t2;
+ t1 = t2 - (s -t1);
+ //[b3,t1] = fast_two_sum(t1,s);
+ b2[0] = t0 + s;
+ t0 = s - (b2[0] - t0);
+ b3[0] = t0 + t1;
+}
+
+static INLINE void
+renorm4(double *c0Ptr, double *c1Ptr, double *c2Ptr, double *c3Ptr) {
+ double s0, s1, s2 = 0.0, s3 = 0.0;
+ double c0 = *c0Ptr;
+
+ if (isinf(c0)) return;
+ fast_two_sum(&s0, c3Ptr, *c2Ptr, *c3Ptr); // s0 = fast_two_sum(*c2Ptr, *c3Ptr, c3Ptr);
+ fast_two_sum(&s0, c2Ptr, *c1Ptr, s0); // s0 = quick_two_sum(*c1Ptr, s0, c2Ptr);
+ fast_two_sum(&c0, c1Ptr, c0, s0); // c0 = quick_two_sum(c0, s0, c1Ptr);
+
+ s0 = c0;
+ s1 = *c1Ptr;
+ if (s1 != 0.0) {
+ fast_two_sum(&s1, &s2, s1, *c2Ptr); // s1 = quick_two_sum(s1, *c2Ptr, &s2);
+ if (s2 != 0.0)
+ fast_two_sum(&s2, &s3, s2, *c3Ptr); // s2 = quick_two_sum(s2, *c3Ptr, &s3);
+ else
+ fast_two_sum(&s1, &s2, s1, *c3Ptr); // s1 = quick_two_sum(s1, *c3Ptr, &s2);
+ } else {
+ fast_two_sum(&s0, &s1, s0, *c2Ptr); // s0 = quick_two_sum(s0, *c2Ptr, &s1);
+ if (s1 != 0.0)
+ fast_two_sum(&s1, &s2, s1, *c3Ptr); // s1 = quick_two_sum(s1, *c3Ptr, &s2);
+ else
+ fast_two_sum(&s0, &s1, s0, *c3Ptr); // s0 = quick_two_sum(s0, *c3Ptr, &s1);
+ }
+
+ *c0Ptr = s0;
+ *c1Ptr = s1;
+ *c2Ptr = s2;
+ *c3Ptr = s3;
+}
+
+//--------------------------------------------------------------------------------------------
+//
+// quad-double square
+//
+//--------------------------------------------------------------------------------------------
+static INLINE void
+qd_sqr(double *c0, double *c1, double *c2, double *c3, double a0, double a1, double a2, double a3)
+{
+ double p01,p11,p02,p03,p12,e00,e01,e11,e02,x,y,w,z,s0,s1,t0,t1;
+ c0[0] = 0.0; c1[0] = 0.0; c2[0] = 0.0; c3[0] = 0.0;
+
+ //O(1) term
+ sqr(&x, &y, a0); c0[0] = x; e00 = y; //O(1) term ok
+
+ //O(eps) term
+ two_prod(&x, &y, a0, a1); p01 = 2.0*x; e01 = 2.0*y;
+
+ //O(eps^2) terms
+ two_prod(&x, &y, a0, a2); p02 = 2.0*x; e02 = 2.0*y;
+ sqr(&x, &y, a1); p11 = x; e11 = y;
+
+ two_sum(&x, &y, p01, e00); c1[0] = x; e00 = y; //O(eps) term ok
+ two_sum(&x, &y, e00, e01); e00 = x; e01 = y;
+ two_sum(&x, &y, p02, p11); p02 = x; p11 = y;
+ two_sum(&x, &y, e00, p02); s0 = x; t0 = y;
+ two_sum(&x, &y, e01, p11); s1 = x; t1 = y;
+ two_sum(&x, &y, s1, t0); s1 = x; t0 = y;
+
+ t0 = t0 + t1;
+
+ fast_two_sum(&x, &y, s1, t0); s1 = x; t0 = y;
+ fast_two_sum(&x, &y, s0, s1); c2[0] = x; t1 = y; //O(eps^2) term ok
+ fast_two_sum(&x, &y, t1, t0); p11 = x; e00 = y;
+
+ //O(eps^3) terms
+ p03 = 2.0 * a0 * a3;
+ p12 = 2.0 * a1 * a2;
+
+ two_sum(&x, &y, p03, p12); p03 = x; p12 = y;
+ two_sum(&x, &y, e02, e11); e02 = x; e11 = y;
+ two_sum(&x, &y, p03, e02); t0 = x; t1 = y;
+
+ t1 = t1 + p12 + e11;
+
+ two_sum(&x, &y, p11, t0); c3[0] = x; p03 = y; //O(eps^3) term ok
+ p03 = p03 + e00 + t1; //O(eps^4) term ok
+
+ renorm(&x, &y, &w, &z, c0[0], c1[0], c2[0], c3[0], p03);
+ c0[0] = x; c1[0] = y; c2[0] = w; c3[0] = z;
+}
+
+//--------------------------------------------------------------------------------------------
+//
+// addition quad-double + double
+//
+//--------------------------------------------------------------------------------------------
+static INLINE void
+qd_add_s(double *c0, double *c1, double *c2, double *c3, 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;
+}
+
+//--------------------------------------------------------------------------------------------
+//
+// addition quad-double + double-double
+//
+//--------------------------------------------------------------------------------------------
+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)
+{
+ 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;
+ 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;
+ return;
+}
+
+//--------------------------------------------------------------------------------------------
+//
+// addition quad-double + quad-double
+//
+//--------------------------------------------------------------------------------------------
+static INLINE void
+qd_add_qd(double *c0, double *c1, double *c2, double *c3, double a0, double a1, double a2, double a3, double b0, double b1, double b2, double b3)
+{
+ double e1,e2,e3,e4,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, 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;
+ three_sum2(&x, &y, c3[0], e3, e1); c3[0] = x; e1 = y;
+ e1 = e1 + e2 + e4;
+ 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;
+}
+
+//--------------------------------------------------------------------------------------------
+//
+// subtraction double - quad-double
+//
+//--------------------------------------------------------------------------------------------
+static INLINE void
+s_sub_qd(double *c0, double *c1, double *c2, double *c3, 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;
+}
+
+//--------------------------------------------------------------------------------------------
+//
+// subtraction quad-double - quad-double
+//
+//--------------------------------------------------------------------------------------------
+static INLINE void
+qd_sub_qd(double *c0, double *c1, double *c2, double *c3, double a0, double a1, double a2, double a3, double b0, double b1, double b2, double b3)
+{
+ double e1,e2,e3,e4,x,y,w,z;
+ b0 = -b0; b1 = -b1; b2 = -b2; b3 = -b3;
+ 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, 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;
+ three_sum2(&x, &y, c3[0], e3, e1); c3[0] = x; e1 = y;
+ e1 = e1 + e2 + e4;
+ 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;
+}
+
+//--------------------------------------------------------------------------------------------
+//
+// multiplication double * quad-double
+//
+//--------------------------------------------------------------------------------------------
+static INLINE void
+s_mul_qd(double *c0, double *c1, double *c2, double *c3, 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;
+ 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;
+}
+
+//--------------------------------------------------------------------------------------------
+//
+// multiplication quad-double * quad-double
+//
+//--------------------------------------------------------------------------------------------
+static INLINE void
+qd_mul_qd(double *c0, double *c1, double *c2, double *c3, double a0, double a1, double a2, double a3, double b0, double b1, double b2, double b3)
+{
+ double p10,p01,p11,p20,p02,e00,e10,e01,e11,e20,e02,x,y,w,z;
+ c0[0] = 0.0; c1[0] = 0.0; c2[0] = 0.0; c3[0] = 0.0;
+
+ //O(1) terms
+ two_prod(&x, &y, a0, b0); c0[0] = x; e00 = y;
+
+ //O(eps) terms
+ two_prod(&x, &y, a0, b1); p01 = x; e01 = y;
+ two_prod(&x, &y, a1, b0); p10 = x; e10 = y;
+ three_sum(&x, &y, &z, p01, p10, e00);
+ c1[0] = x; //O(eps)
+ p10 = y; //O(eps^2)
+ p01 = z; //O(eps^3)
+
+ //O(eps^2) terms
+ two_prod(&x, &y, a0, b2); p02 = x; e02 = y;
+ two_prod(&x, &y, a1, b1); p11 = x; e11 = y;
+ two_prod(&x, &y, a2, b0); p20 = x; e20 = y;
+
+ //six three sum for p10, e01, e10, p02, p11, p20
+ three_sum(&x, &y, &z, p10, e01, e10); p10 = x; e01 = y; e10 = z;
+ three_sum(&x, &y, &z, p02, p11, p20); p02 = x; p11 = y; p20 = z;
+ two_sum(&x, &y, p02, p10); c2[0] = x; p10 = y;
+ two_sum(&x, &y, p11, e01); p11 = x; e01 = y;
+ two_sum(&x, &y, p10, p11); p10 = x; p11 = y;
+
+ e10 = e10 + p20 + e01 + p11; //O(eps^4) terms
+
+ //O(eps^3) terms
+ c3[0] = p10 + a0*b3 + a1*b2 + a2*b1 + a3*b0 + e02 + e11 + e20;
+
+ renorm(&x, &y, &w, &z, c0[0], c1[0], c2[0], c3[0], e10);
+ c0[0] = x;
+ c1[0] = y;
+ c2[0] = w;
+ c3[0] = z;
+}
+
+//--------------------------------------------------------------------------------------------
+//
+// division quad-double / double
+//
+//--------------------------------------------------------------------------------------------
+static INLINE void
+
+qd_div_s(double *c0, double *c1, double *c2, double *c3, double a0, double a1, double a2, double a3, double b)
+{
+ double x,y,w,z,t0,t1,r0,r1,r2,r3,e;
+ c0[0] = 0.0; c1[0] = 0.0; c2[0] = 0.0; c3[0] = 0.0;
+
+ c0[0] = a0/b;
+ // reminder a - c_0*b
+ two_prod(&x, &y, c0[0], b);
+ t0 = -x;
+ t1 = -y;
+ //qd subtruction (a - t)
+ qd_add_dd(&x, &y, &w, &z, a0, a1, a2, a3, t0, t1);
+ r0 = x; r1 = y; r2 = w; r3 = z;
+
+ c1[0] = r0/b;
+ // reminder r - c_1*b
+ two_prod(&x, &y, c1[0], b);
+ t0 = -x;
+ t1 = -y;
+ //qd subtruction (r - t)
+ qd_add_dd(&x, &y, &w, &z, r0, r1, r2, r3, t0, t1);
+ r0 = x; r1 = y; r2 = w; r3 = z;
+
+ c2[0] = r0/b;
+ // reminder r - c_2*b
+ two_prod(&x, &y, c2[0], b);
+ t0 = -x;
+ t1 = -y;
+ //qd subtruction (r - t)
+ qd_add_dd(&x, &y, &w, &z, r0, r1, r2, r3, t0, t1);
+ r0 = x; r1 = y; r2 = w; r3 = z;
+
+ c3[0] = r0/b;
+ // reminder r - c_3*b
+ two_prod(&x, &y, c3[0], b);
+ t0 = -x;
+ t1 = -y;
+ //qd subtruction (r - t)
+ qd_add_dd(&x, &y, &w, &z, r0, r1, r2, r3, t0, t1);
+ r0 = x; r1 = y; r2 = w; r3 = z;
+
+ e = r0/b;
+ 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;
+}
+
+//--------------------------------------------------------------------------------------------
+//
+// division quad-double / quad-double
+//
+//--------------------------------------------------------------------------------------------
+static INLINE void
+qd_div_qd(double *c0, double *c1, double *c2, double *c3, double a0, double a1, double a2, double a3, double b0, double b1, double b2, double b3)
+{
+ double x,y,w,z,t0,t1,t2,t3,r0,r1,r2,r3,e;
+
+ c0[0] = 0.0; c1[0] = 0.0; c2[0] = 0.0; c3[0] = 0.0;
+
+ c0[0] = a0/b0;
+ // reminder a - c_0*b
+ //multiplication
+ s_mul_qd(&x, &y, &w, &z, c0[0], b0, b1, b2, b3);
+ t0 = -x; t1 = -y; t2 = -w; t3 = -z;
+ //qd subtruction (a - t)
+ qd_add_qd(&x, &y, &w, &z, a0, a1, a2, a3, t0, t1, t2, t3);
+ r0 = x; r1 = y; r2 = w; r3 = z;
+
+ c1[0] = r0/b0;
+ // reminder r - c_1*b
+ //multiplication
+ s_mul_qd(&x, &y, &w, &z, c1[0], b0, b1, b2, b3);
+ t0 = -x; t1 = -y; t2 = -w; t3 = -z;
+ //qd subtruction (r - t)
+ qd_add_qd(&x, &y, &w, &z, r0, r1, r2, r3, t0, t1, t2, t3);
+ r0 = x; r1 = y; r2 = w; r3 = z;
+
+ c2[0] = r0/b0;
+ // reminder r - c_2*b
+ //multiplication
+ s_mul_qd(&x, &y, &w, &z, c2[0], b0, b1, b2, b3);
+ t0 = -x; t1 = -y; t2 = -w; t3 = -z;
+ //qd subtruction (r - t)
+ qd_add_qd(&x, &y, &w, &z, r0, r1, r2, r3, t0, t1, t2, t3);
+ r0 = x; r1 = y; r2 = w; r3 = z;
+
+ c3[0] = r0/b0;
+ // reminder r - c_3*b
+ //multiplication
+ s_mul_qd(&x, &y, &w, &z, c3[0], b0, b1, b2, b3);
+ t0 = -x; t1 = -y; t2 = -w; t3 = -z;
+ //qd subtruction (r - t)
+ qd_add_qd(&x, &y, &w, &z, r0, r1, r2, r3, t0, t1, t2, t3);
+ r0 = x; r1 = y; r2 = w; r3 = z;
+
+ e = r0/b0;
+ 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;
+}
+
+//--------------------------------------------------------------------------------------------
+//
+// division double / quad-double sloppy
+//
+//--------------------------------------------------------------------------------------------
+static INLINE void
+s_div_qd(double *c0, double *c1, double *c2, double *c3, double a, double b0, double b1, double b2, double b3)
+{
+ double x,y,w,z,t0,t1,t2,t3,r0,r1,r2,r3;
+
+ c0[0] = 0.0; c1[0] = 0.0; c2[0] = 0.0; c3[0] = 0.0;
+
+ c0[0] = a/b0;
+ // reminder a - c_0*b
+ //multiplication
+ s_mul_qd(&x, &y, &w, &z, c0[0], b0, b1, b2, b3);
+ t0 = -x; t1 = -y; t2 = -w; t3 = -z;
+ //qd subtruction (a - t)
+ qd_add_s(&x, &y, &w, &z, t0, t1, t2, t3, a);
+
+ r0 = x; r1 = y; r2 = w; r3 = z;
+
+
+ c1[0] = r0/b0;
+ // reminder r - c_1*b
+ s_mul_qd(&x, &y, &w, &z, c1[0], b0, b1, b2, b3);
+ t0 = -x; t1 = -y; t2 = -w; t3 = -z;
+ //qd subtruction (r - t)
+ qd_add_qd(&x, &y, &w, &z, r0, r1, r2, r3, t0, t1, t2, t3);
+ r0 = x; r1 = y; r2 = w; r3 = z;
+
+ c2[0] = r0/b0;
+ // reminder r - c_2*b
+ s_mul_qd(&x, &y, &w, &z, c2[0], b0, b1, b2, b3);
+ t0 = -x; t1 = -y; t2 = -w; t3 = -z;
+ //qd subtruction (r - t)
+ qd_add_qd(&x, &y, &w, &z, r0, r1, r2, r3, t0, t1, t2, t3);
+ r0 = x; r1 = y; r2 = w; r3 = z;
+
+ c3[0] = r0/b0;
+
+ renorm(&x, &y, &w, &z, c0[0], c1[0], c2[0], c3[0], 0.0);
+ c0[0] = x;
+ c1[0] = y;
+ c2[0] = w;
+ c3[0] = z;
+}
+
+static INLINE void
+qdsqrt(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;
+
+ c0[0] = 0.0; c1[0] = 0.0; c2[0] = 0.0; c3[0] = 0.0;
+
+ c0[0] = 1.0/sqrt(a0);
+ h0 = 0.5*a0; h1 = 0.5*a1; h2 = 0.5*a2; h3 = 0.5*a3;
+
+ qd_sqr(&x0, &x1, &x2, &x3, c0[0], c1[0], c2[0], c3[0]); //x_k^2
+ qd_mul_qd(&p, &q, &r, &s, h0, h1, h2, h3, x0, x1, x2, x3); //0.5 * a * x_k^2
+ x0 = -p; x1 = -q; x2 = -r; x3 = -s;
+ qd_add_s(&p, &q, &r, &s, x0, x1, x2, x3, 0.5); //0.5 - 0.5 * a * x_k^2
+ qd_mul_qd(&x0, &x1, &x2, &x3, p, q, r, s, c0[0], c1[0], c2[0], c3[0]); //(0.5 - 0.5 * a * x_k^2)*x_k
+ qd_add_qd(&p, &q, &r, &s, c0[0], c1[0], c2[0], c3[0], x0, x1, x2, x3); //x_k+1 = x_k + (0.5 - 0.5 * a * x_k^2)*x_k
+ c0[0] = p; c1[0] = q; c2[0] = r; c3[0] = s;
+
+ qd_sqr(&x0, &x1, &x2, &x3, c0[0], c1[0], c2[0], c3[0]); //x_k^2
+ qd_mul_qd(&p, &q, &r, &s, h0, h1, h2, h3, x0, x1, x2, x3); //0.5 * a * x_k^2
+ x0 = -p; x1 = -q; x2 = -r; x3 = -s;
+ qd_add_s(&p, &q, &r, &s, x0, x1, x2, x3, 0.5); //0.5 - 0.5 * a * x_k^2
+ qd_mul_qd(&x0, &x1, &x2, &x3, p, q, r, s, c0[0], c1[0], c2[0], c3[0]); //(0.5 - 0.5 * a * x_k^2)*x_k
+ qd_add_qd(&p, &q, &r, &s, c0[0], c1[0], c2[0], c3[0], x0, x1, x2, x3); //x_k+1 = x_k + (0.5 - 0.5 * a * x_k^2)*x_k
+ c0[0] = p; c1[0] = q; c2[0] = r; c3[0] = s;
+
+ qd_sqr(&x0, &x1, &x2, &x3, c0[0], c1[0], c2[0], c3[0]); //x_k^2
+ qd_mul_qd(&p, &q, &r, &s, h0, h1, h2, h3, x0, x1, x2, x3); //0.5 * a * x_k^2
+ x0 = -p; x1 = -q; x2 = -r; x3 = -s;
+ qd_add_s(&p, &q, &r, &s, x0, x1, x2, x3, 0.5); //0.5 - 0.5 * a * x_k^2
+ qd_mul_qd(&x0, &x1, &x2, &x3, p, q, r, s, c0[0], c1[0], c2[0], c3[0]); //(0.5 - 0.5 * a * x_k^2)*x_k
+ qd_add_qd(&p, &q, &r, &s, c0[0], c1[0], c2[0], c3[0], x0, x1, x2, x3); //x_k+1 = x_k + (0.5 - 0.5 * a * x_k^2)*x_k
+
+ qd_mul_qd(&x0, &x1, &x2, &x3, a0, a1, a2, a3, p, q, r, s); //(0.5 - 0.5 * a * x_k^2)*x_k*a
+ c0[0] = x0; c1[0] = x1; c2[0] = x2; c3[0] = x3;
+}
+
+static void
+qdpow(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;
+
+ c0[0] = 0.0; c1[0] = 0.0; c2[0] = 0.0; c3[0] = 0.0;
+
+ if (p == 0) {
+ c0[0] = 1.0;
+ } else {
+ r0 = a0; r1 = a1; r2 = a2; r3 = a3;
+ c0[0] = 1.0;
+ abs_p = abs(p);
+
+ if (abs_p > 1) {
+ while (abs_p > 0) {
+ if ((abs_p % 2)==1) {
+ qd_mul_qd(&x, &y, &w, &z, c0[0], c1[0], c2[0], c3[0], r0, r1, r2, r3);
+ c0[0] = x; c1[0] = y; c2[0] = w; c3[0] = z;
+ }
+ abs_p = abs_p / 2;
+ if (abs_p > 0) {
+ qd_sqr(&x, &y, &w, &z, r0, r1, r2, r3);
+ r0 = x; r1 = y; r2 = w; r3 = z;
+ }
+ }
+ } else {
+ c0[0] = r0; c1[0] = r1; c2[0] = r2; c3[0] = r3;
+ }
+ if (p < 0) {
+ s_div_qd(&x, &y, &w, &z, 1.0, c0[0], c1[0], c2[0], c3[0]);
+ c0[0] = x; c1[0] = y; c2[0] = w; c3[0] = z;
+ }
+ }
+}
+
+// round to nearest integer
+#define round(x) (floor((x)+0.5))
+
+static INLINE void
+nint_qd(double *x0, double *x1, double *x2, double *x3, double a0, double a1, double a2, double a3)
+{
+ x0[0]=round(a0);
+ x1[0]=0.0; x2[0]=0.0; x3[0]=0.0;
+
+ if(x0[0]==a0) {
+ x1[0]=round(a1);
+ if(x1[0]==a1) {
+ x2[0]=round(a2);
+ if(x2[0]==a2) {
+ x3[0]=round(a3);
+ }
+ else {
+ if(((int)fabs(x2[0]-a2)==0.5) && (a3<0.0)) {
+ x2[0]=x2[0]-1.0;
+ }
+ }
+ }
+ else {
+ if(((int)fabs(x1[0]-a1)==0.5) && (a2<0.0)) {
+ x1[0]=x1[0]-1.0;
+ }
+ }
+ }
+ else {
+ if(((int)fabs(x0[0]-a0)==0.5) && (a1<0.0)) {
+ x0[0]=x0[0]-1.0;
+ }
+ }
+ renorm(&x0[0],&x1[0],&x2[0],&x3[0],x0[0],x1[0],x2[0],x3[0],0.0);
+ return;
+}
+
+static double s_table[256][4]= {
+ {3.0679567629659761e-03, 1.2690279085455925e-19, 5.2879464245328389e-36, -1.7820334081955298e-52},
+ {6.1358846491544753e-03, 9.0545257482474933e-20, 1.6260113133745320e-37, -9.7492001208767410e-55},
+ {9.2037547820598194e-03, -1.2136591693535934e-19, 5.5696903949425567e-36, 1.2505635791936951e-52},
+ {1.2271538285719925e-02, 6.9197907640283170e-19, -4.0203726713435555e-36, -2.0688703606952816e-52},
+ {1.5339206284988102e-02, -8.4462578865401696e-19, 4.6535897505058629e-35, -1.3923682978570467e-51},
+ {1.8406729905804820e-02, 7.4195533812833160e-19, 3.9068476486787607e-35, 3.6393321292898614e-52},
+ {2.1474080275469508e-02, -4.5407960207688566e-19, -2.2031770119723005e-35, 1.2709814654833741e-52},
+ {2.4541228522912288e-02, -9.1868490125778782e-20, 4.8706148704467061e-36, -2.8153947855469224e-52},
+ {2.7608145778965743e-02, -1.5932358831389269e-18, -7.0475416242776030e-35, -2.7518494176602744e-51},
+ {3.0674803176636626e-02, -1.6936054844107918e-20, -2.0039543064442544e-36, -1.6267505108658196e-52},
+ {3.3741171851377587e-02, -2.0096074292368340e-18, -1.3548237016537134e-34, 6.5554881875899973e-51},
+ {3.6807222941358832e-02, 6.1060088803529842e-19, -4.0448721259852727e-35, -2.1111056765671495e-51},
+ {3.9872927587739811e-02, 4.6657453481183289e-19, 3.4119333562288684e-35, 2.4007534726187511e-51},
+ {4.2938256934940820e-02, 2.8351940588660907e-18, 1.6991309601186475e-34, 6.8026536098672629e-51},
+ {4.6003182130914630e-02, -1.1182813940157788e-18, 7.5235020270378946e-35, 4.1187304955493722e-52},
+ {4.9067674327418015e-02, -6.7961037205182801e-19, -4.4318868124718325e-35, -9.9376628132525316e-52},
+ {5.2131704680283324e-02, -2.4243695291953779e-18, -1.3675405320092298e-34, -8.3938137621145070e-51},
+ {5.5195244349689941e-02, -1.3340299860891103e-18, -3.4359574125665608e-35, 1.1911462755409369e-51},
+ {5.8258264500435759e-02, 2.3299905496077492e-19, 1.9376108990628660e-36, -5.1273775710095301e-53},
+ {6.1320736302208578e-02, -5.1181134064638108e-19, -4.2726335866706313e-35, 2.6368495557440691e-51},
+ {6.4382630929857465e-02, -4.2325997000052705e-18, 3.3260117711855937e-35, 1.4736267706718352e-51},
+ {6.7443919563664065e-02, -6.9221796556983636e-18, 1.5909286358911040e-34, -7.8828946891835218e-51},
+ {7.0504573389613870e-02, -6.8552791107342883e-18, -1.9961177630841580e-34, 2.0127129580485300e-50},
+ {7.3564563599667426e-02, -2.7784941506273593e-18, -9.1240375489852821e-35, -1.9589752023546795e-51},
+ {7.6623861392031492e-02, 2.3253700287958801e-19, -1.3186083921213440e-36, -4.9927872608099673e-53},
+ {7.9682437971430126e-02, -4.4867664311373041e-18, 2.8540789143650264e-34, 2.8491348583262741e-51},
+ {8.2740264549375692e-02, 1.4735983530877760e-18, 3.7284093452233713e-35, 2.9024430036724088e-52},
+ {8.5797312344439894e-02, -3.3881893830684029e-18, -1.6135529531508258e-34, 7.7294651620588049e-51},
+ {8.8853552582524600e-02, -3.7501775830290691e-18, 3.7543606373911573e-34, 2.2233701854451859e-50},
+ {9.1908956497132724e-02, 4.7631594854274564e-18, 1.5722874642939344e-34, -4.8464145447831456e-51},
+ {9.4963495329639006e-02, -6.5885886400417564e-18, -2.1371116991641965e-34, 1.3819370559249300e-50},
+ {9.8017140329560604e-02, -1.6345823622442560e-18, -1.3209238810006454e-35, -3.5691060049117942e-52},
+ {1.0106986275482782e-01, 3.3164325719308656e-18, -1.2004224885132282e-34, 7.2028828495418631e-51},
+ {1.0412163387205457e-01, 6.5760254085385100e-18, 1.7066246171219214e-34, -4.9499340996893514e-51},
+ {1.0717242495680884e-01, 6.4424044279026198e-18, -8.3956976499698139e-35, -4.0667730213318321e-51},
+ {1.1022220729388306e-01, -5.6789503537823233e-19, 1.0380274792383233e-35, 1.5213997918456695e-52},
+ {1.1327095217756435e-01, 2.7100481012132900e-18, 1.5323292999491619e-35, 4.9564432810360879e-52},
+ {1.1631863091190477e-01, 1.0294914877509705e-18, -9.3975734948993038e-35, 1.3534827323719708e-52},
+ {1.1936521481099137e-01, -3.9500089391898506e-18, 3.5317349978227311e-34, 1.8856046807012275e-51},
+ {1.2241067519921620e-01, 2.8354501489965335e-18, 1.8151655751493305e-34, -2.8716592177915192e-51},
+ {1.2545498341154623e-01, 4.8686751763148235e-18, 5.9878105258097936e-35, -3.3534629098722107e-51},
+ {1.2849811079379317e-01, 3.8198603954988802e-18, -1.8627501455947798e-34, -2.4308161133527791e-51},
+ {1.3154002870288312e-01, -5.0039708262213813e-18, -1.2983004159245552e-34, -4.6872034915794122e-51},
+ {1.3458070850712620e-01, -9.1670359171480699e-18, 1.5916493007073973e-34, 4.0237002484366833e-51},
+ {1.3762012158648604e-01, 6.6253255866774482e-18, -2.3746583031401459e-34, -9.3703876173093250e-52},
+ {1.4065823933284924e-01, -7.9193932965524741e-18, 6.0972464202108397e-34, 2.4566623241035797e-50},
+ {1.4369503315029444e-01, 1.1472723016618666e-17, -5.1884954557576435e-35, -4.2220684832186607e-51},
+ {1.4673047445536175e-01, 3.7269471470465677e-18, 3.7352398151250827e-34, -4.0881822289508634e-51},
+ {1.4976453467732151e-01, 8.0812114131285151e-18, 1.2979142554917325e-34, 9.9380667487736254e-51},
+ {1.5279718525844344e-01, -7.6313573938416838e-18, 5.7714690450284125e-34, -3.7731132582986687e-50},
+ {1.5582839765426523e-01, 3.0351307187678221e-18, -1.0976942315176184e-34, 7.8734647685257867e-51},
+ {1.5885814333386145e-01, -4.0163200573859079e-18, -9.2840580257628812e-35, -2.8567420029274875e-51},
+ {1.6188639378011183e-01, 1.1850519643573528e-17, -5.0440990519162957e-34, 3.0510028707928009e-50},
+ {1.6491312048996992e-01, -7.0405288319166738e-19, 3.3211107491245527e-35, 8.6663299254686031e-52},
+ {1.6793829497473117e-01, 5.4284533721558139e-18, -3.3263339336181369e-34, -1.8536367335123848e-50},
+ {1.7096188876030122e-01, 9.1919980181759094e-18, -6.7688743940982606e-34, -1.0377711384318389e-50},
+ {1.7398387338746382e-01, 5.8151994618107928e-18, -1.6751014298301606e-34, -6.6982259797164963e-51},
+ {1.7700422041214875e-01, 6.7329300635408167e-18, 2.8042736644246623e-34, 3.6786888232793599e-51},
+ {1.8002290140569951e-01, 7.9701826047392143e-18, -7.0765920110524977e-34, 1.9622512608461784e-50},
+ {1.8303988795514095e-01, 7.7349918688637383e-18, -4.4803769968145083e-34, 1.1201148793328890e-50},
+ {1.8605515166344666e-01, -1.2564893007679552e-17, 7.5953844248530810e-34, -3.8471695132415039e-51},
+ {1.8906866414980622e-01, -7.6208955803527778e-18, -4.4792298656662981e-34, -4.4136824096645007e-50},
+ {1.9208039704989244e-01, 4.3348343941174903e-18, -2.3404121848139937e-34, 1.5789970962611856e-50},
+ {1.9509032201612828e-01, -7.9910790684617313e-18, 6.1846270024220713e-34, -3.5840270918032937e-50},
+ {1.9809841071795359e-01, -1.8434411800689445e-18, 1.4139031318237285e-34, 1.0542811125343809e-50},
+ {2.0110463484209190e-01, 1.1010032669300739e-17, -3.9123576757413791e-34, 2.4084852500063531e-51},
+ {2.0410896609281687e-01, 6.0941297773957752e-18, -2.8275409970449641e-34, 4.6101008563532989e-51},
+ {2.0711137619221856e-01, -1.0613362528971356e-17, 2.2456805112690884e-34, 1.3483736125280904e-50},
+ {2.1011183688046961e-01, 1.1561548476512844e-17, 6.0355905610401254e-34, 3.3329909618405675e-50},
+ {2.1311031991609136e-01, 1.2031873821063860e-17, -3.4142699719695635e-34, -1.2436262780241778e-50},
+ {2.1610679707621952e-01, -1.0111196082609117e-17, 7.2789545335189643e-34, -2.9347540365258610e-50},
+ {2.1910124015686980e-01, -3.6513812299150776e-19, -2.3359499418606442e-35, 3.1785298198458653e-52},
+ {2.2209362097320354e-01, -3.0337210995812162e-18, 6.6654668033632998e-35, 2.0110862322656942e-51},
+ {2.2508391135979283e-01, 3.9507040822556510e-18, 2.4287993958305375e-35, 5.6662797513020322e-52},
+ {2.2807208317088573e-01, 8.2361837339258012e-18, 6.9786781316397937e-34, -6.4122962482639504e-51},
+ {2.3105810828067111e-01, 1.0129787149761869e-17, -6.9359234615816044e-34, -2.8877355604883782e-50},
+ {2.3404195858354343e-01, -6.9922402696101173e-18, -5.7323031922750280e-34, 5.3092579966872727e-51},
+ {2.3702360599436720e-01, 8.8544852285039918e-18, 1.3588480826354134e-34, 1.0381022520213867e-50},
+ {2.4000302244874150e-01, -1.2137758975632164e-17, -2.6448807731703891e-34, -1.9929733800670473e-51},
+ {2.4298017990326390e-01, -8.7514315297196632e-18, -6.5723260373079431e-34, -1.0333158083172177e-50},
+ {2.4595505033579462e-01, -1.1129044052741832e-17, 4.3805998202883397e-34, 1.2219399554686291e-50},
+ {2.4892760574572018e-01, -8.1783436100020990e-18, 5.5666875261111840e-34, 3.8080473058748167e-50},
+ {2.5189781815421697e-01, -1.7591436032517039e-17, -1.0959681232525285e-33, 5.6209426020232456e-50},
+ {2.5486565960451457e-01, -1.3602299806901461e-19, -6.0073844642762535e-36, -3.0072751311893878e-52},
+ {2.5783110216215899e-01, 1.8480038630879957e-17, 3.3201664714047599e-34, -5.5547819290576764e-51},
+ {2.6079411791527551e-01, 4.2721420983550075e-18, 5.6782126934777920e-35, 3.1428338084365397e-51},
+ {2.6375467897483140e-01, -1.8837947680038700e-17, 1.3720129045754794e-33, -8.2763406665966033e-50},
+ {2.6671275747489837e-01, 2.0941222578826688e-17, -1.1303466524727989e-33, 1.9954224050508963e-50},
+ {2.6966832557291509e-01, 1.5765657618133259e-17, -6.9696142173370086e-34, -4.0455346879146776e-50},
+ {2.7262135544994898e-01, 7.8697166076387850e-18, 6.6179388602933372e-35, -2.7642903696386267e-51},
+ {2.7557181931095814e-01, 1.9320328962556582e-17, 1.3932094180100280e-33, 1.3617253920018116e-50},
+ {2.7851968938505312e-01, -1.0030273719543544e-17, 7.2592115325689254e-34, -1.0068516296655851e-50},
+ {2.8146493792575800e-01, -1.2322299641274009e-17, -1.0564788706386435e-34, 7.5137424251265885e-51},
+ {2.8440753721127182e-01, 2.2209268510661475e-17, -9.1823095629523708e-34, -5.2192875308892218e-50},
+ {2.8734745954472951e-01, 1.5461117367645717e-17, -6.3263973663444076e-34, -2.2982538416476214e-50},
+ {2.9028467725446239e-01, -1.8927978707774251e-17, 1.1522953157142315e-33, 7.4738655654716596e-50},
+ {2.9321916269425863e-01, 2.2385430811901833e-17, 1.3662484646539680e-33, -4.2451325253996938e-50},
+ {2.9615088824362384e-01, -2.0220736360876938e-17, -7.9252212533920413e-35, -2.8990577729572470e-51},
+ {2.9907982630804048e-01, 1.6701181609219447e-18, 8.6091151117316292e-35, 3.9931286230012102e-52},
+ {3.0200594931922808e-01, -1.7167666235262474e-17, 2.3336182149008069e-34, 8.3025334555220004e-51},
+ {3.0492922973540243e-01, -2.2989033898191262e-17, -1.4598901099661133e-34, 3.7760487693121827e-51},
+ {3.0784964004153487e-01, 2.7074088527245185e-17, 1.2568858206899284e-33, 7.2931815105901645e-50},
+ {3.1076715274961147e-01, 2.0887076364048513e-17, -3.0130590791065942e-34, 1.3876739009935179e-51},
+ {3.1368174039889146e-01, 1.4560447299968912e-17, 3.6564186898011595e-34, 1.1654264734999375e-50},
+ {3.1659337555616585e-01, 2.1435292512726283e-17, 1.2338169231377316e-33, 3.3963542100989293e-50},
+ {3.1950203081601569e-01, -1.3981562491096626e-17, 8.1730000697411350e-34, -7.7671096270210952e-50},
+ {3.2240767880106985e-01, -4.0519039937959398e-18, 3.7438302780296796e-34, 8.7936731046639195e-51},
+ {3.2531029216226293e-01, 7.9171249463765892e-18, -6.7576622068146391e-35, 2.3021655066929538e-51},
+ {3.2820984357909255e-01, -2.6693140719641896e-17, 7.8928851447534788e-34, 2.5525163821987809e-51},
+ {3.3110630575987643e-01, -2.7469465474778694e-17, -1.3401245916610206e-33, 6.5531762489976163e-50},
+ {3.3399965144200938e-01, 2.2598986806288142e-17, 7.8063057192586115e-34, 2.0427600895486683e-50},
+ {3.3688985339222005e-01, -4.2000940033475092e-19, -2.9178652969985438e-36, -1.1597376437036749e-52},
+ {3.3977688440682685e-01, 6.6028679499418282e-18, 1.2575009988669683e-34, 2.5569067699008304e-51},
+ {3.4266071731199438e-01, 1.9261518449306319e-17, -9.2754189135990867e-34, 8.5439996687390166e-50},
+ {3.4554132496398904e-01, 2.7251143672916123e-17, 7.0138163601941737e-34, -1.4176292197454015e-50},
+ {3.4841868024943456e-01, 3.6974420514204918e-18, 3.5532146878499996e-34, 1.9565462544501322e-50},
+ {3.5129275608556715e-01, -2.2670712098795844e-17, -1.6994216673139631e-34, -1.2271556077284517e-50},
+ {3.5416352542049040e-01, -1.6951763305764860e-17, 1.2772331777814617e-33, -3.3703785435843310e-50},
+ {3.5703096123343003e-01, -4.8218191137919166e-19, -4.1672436994492361e-35, -7.1531167149364352e-52},
+ {3.5989503653498817e-01, -1.7601687123839282e-17, 1.3375125473046791e-33, 7.9467815593584340e-50},
+ {3.6275572436739723e-01, -9.1668352663749849e-18, -7.4317843956936735e-34, -2.0199582511804564e-50},
+ {3.6561299780477385e-01, 1.6217898770457546e-17, 1.1286970151961055e-33, -7.1825287318139010e-50},
+ {3.6846682995337232e-01, 1.0463640796159268e-17, 2.0554984738517304e-35, 1.0441861305618769e-51},
+ {3.7131719395183754e-01, 3.4749239648238266e-19, -7.5151053042866671e-37, -2.8153468438650851e-53},
+ {3.7416406297145799e-01, 8.0114103761962118e-18, 5.3429599813406052e-34, 1.0351378796539210e-50},
+ {3.7700741021641826e-01, -2.7255302041956930e-18, 6.3646586445018137e-35, 8.3048657176503559e-52},
+ {3.7984720892405116e-01, 9.9151305855172370e-18, 4.8761409697224886e-34, 1.4025084000776705e-50},
+ {3.8268343236508978e-01, -1.0050772696461588e-17, -2.0605316302806695e-34, -1.2717724698085205e-50},
+ {3.8551605384391885e-01, 1.5177665396472313e-17, 1.4198230518016535e-33, 5.8955167159904235e-50},
+ {3.8834504669882630e-01, -1.0053770598398717e-17, 7.5942999255057131e-34, -3.1967974046654219e-50},
+ {3.9117038430225387e-01, 1.7997787858243995e-17, -1.0613482402609856e-33, -5.4582148817791032e-50},
+ {3.9399204006104810e-01, 9.7649241641239336e-18, -2.1233599441284617e-34, -5.5529836795340819e-51},
+ {3.9680998741671031e-01, 2.0545063670840126e-17, 6.1347058801922842e-34, 1.0733788150636430e-50},
+ {3.9962419984564684e-01, -1.5065497476189372e-17, -9.9653258881867298e-34, -5.7524323712725355e-50},
+ {4.0243465085941843e-01, 1.0902619339328270e-17, 7.3998528125989765e-34, 2.2745784806823499e-50},
+ {4.0524131400498986e-01, 9.9111401942899884e-18, -2.5169070895434648e-34, 9.2772984818436573e-53},
+ {4.0804416286497869e-01, -7.0006015137351311e-18, -1.4108207334268228e-34, 1.5175546997577136e-52},
+ {4.1084317105790397e-01, -2.4219835190355499e-17, -1.1418902925313314e-33, -2.0996843165093468e-50},
+ {4.1363831223843456e-01, -1.0393984940597871e-17, -1.1481681174503880e-34, -2.0281052851028680e-51},
+ {4.1642956009763721e-01, -2.5475580413131732e-17, -3.4482678506112824e-34, 7.1788619351865480e-51},
+ {4.1921688836322396e-01, -4.2232463750110590e-18, -3.6053023045255790e-34, -2.2209673210025631e-50},
+ {4.2200027079979968e-01, 4.3543266994128527e-18, 3.1734310272251190e-34, -1.3573247980738668e-50},
+ {4.2477968120910881e-01, 2.7462312204277281e-17, -4.6552847802111948e-34, 6.5961781099193122e-51},
+ {4.2755509343028208e-01, 9.4111898162954726e-18, -1.7446682426598801e-34, -2.2054492626480169e-51},
+ {4.3032648134008261e-01, 2.2259686974092690e-17, 8.5972591314085075e-34, -2.9420897889003020e-50},
+ {4.3309381885315196e-01, 1.1224283329847517e-17, 5.3223748041075651e-35, 5.3926192627014212e-51},
+ {4.3585707992225547e-01, 1.6230515450644527e-17, -6.4371449063579431e-35, -6.9102436481386757e-51},
+ {4.3861623853852766e-01, -2.0883315831075090e-17, -1.4259583540891877e-34, 6.3864763590657077e-52},
+ {4.4137126873171667e-01, 2.2360783886964969e-17, 1.1864769603515770e-34, -3.8087003266189232e-51},
+ {4.4412214457042926e-01, -2.4218874422178315e-17, 2.2205230838703907e-34, 9.2133035911356258e-51},
+ {4.4686884016237421e-01, -1.9222136142309382e-17, -4.4425678589732049e-35, -1.3673609292149535e-51},
+ {4.4961132965460660e-01, 4.8831924232035243e-18, 2.7151084498191381e-34, -1.5653993171613154e-50},
+ {4.5234958723377089e-01, -1.4827977472196122e-17, -7.6947501088972324e-34, 1.7656856882031319e-50},
+ {4.5508358712634384e-01, -1.2379906758116472e-17, 5.5289688955542643e-34, -8.5382312840209386e-51},
+ {4.5781330359887723e-01, -8.4554254922295949e-18, -6.3770394246764263e-34, 3.1778253575564249e-50},
+ {4.6053871095824001e-01, 1.8488777492177872e-17, -1.0527732154209725e-33, 3.3235593490947102e-50},
+ {4.6325978355186020e-01, -7.3514924533231707e-18, 6.7175396881707035e-34, 3.9594127612123379e-50},
+ {4.6597649576796618e-01, -3.3023547778235135e-18, 3.4904677050476886e-35, 3.4483855263874246e-51},
+ {4.6868882203582796e-01, -2.2949251681845054e-17, -1.1364757641823658e-33, 6.8840522501918612e-50},
+ {4.7139673682599764e-01, 6.5166781360690130e-18, 2.9457546966235984e-34, -6.2159717738836630e-51},
+ {4.7410021465055002e-01, -8.1451601548978075e-18, -3.4789448555614422e-34, -1.1681943974658508e-50},
+ {4.7679923006332214e-01, -1.0293515338305794e-17, -3.6582045008369952e-34, 1.7424131479176475e-50},
+ {4.7949375766015301e-01, 1.8419999662684771e-17, -1.3040838621273312e-33, 1.0977131822246471e-50},
+ {4.8218377207912277e-01, -2.5861500925520442e-17, -6.2913197606500007e-36, 4.0802359808684726e-52},
+ {4.8486924800079112e-01, -1.8034004203262245e-17, -3.5244276906958044e-34, -1.7138318654749246e-50},
+ {4.8755016014843594e-01, 1.4231090931273653e-17, -1.8277733073262697e-34, -1.5208291790429557e-51},
+ {4.9022648328829116e-01, -5.1496145643440404e-18, -3.6903027405284104e-34, 1.5172940095151304e-50},
+ {4.9289819222978404e-01, -1.0257831676562186e-18, 6.9520817760885069e-35, -2.4260961214090389e-51},
+ {4.9556526182577254e-01, -9.4323241942365362e-18, 3.1212918657699143e-35, 4.2009072375242736e-52},
+ {4.9822766697278187e-01, -1.6126383830540798e-17, -1.5092897319298871e-33, 1.1049298890895917e-50},
+ {5.0088538261124083e-01, -3.9604015147074639e-17, -2.2208395201898007e-33, 1.3648202735839417e-49},
+ {5.0353838372571758e-01, -1.6731308204967497e-17, -1.0140233644074786e-33, 4.0953071937671477e-50},
+ {5.0618664534515534e-01, -4.8321592986493711e-17, 9.2858107226642252e-34, 4.2699802401037005e-50},
+ {5.0883014254310699e-01, 4.7836968268014130e-17, -1.0727022928806035e-33, 2.7309374513672757e-50},
+ {5.1146885043797041e-01, -1.3088001221007579e-17, 4.0929033363366899e-34, -3.7952190153477926e-50},
+ {5.1410274419322177e-01, -4.5712707523615624e-17, 1.5488279442238283e-33, -2.5853959305521130e-50},
+ {5.1673179901764987e-01, 8.3018617233836515e-18, 5.8251027467695202e-34, -2.2812397190535076e-50},
+ {5.1935599016558964e-01, -5.5331248144171145e-17, -3.1628375609769026e-35, -2.4091972051188571e-51},
+ {5.2197529293715439e-01, -4.6555795692088883e-17, 4.6378980936850430e-34, -3.3470542934689532e-51},
+ {5.2458968267846895e-01, -4.3068869040082345e-17, -4.2013155291932055e-34, -1.5096069926700274e-50},
+ {5.2719913478190139e-01, -4.2202983480560619e-17, 8.5585916184867295e-34, 7.9974339336732307e-50},
+ {5.2980362468629472e-01, -4.8067841706482342e-17, 5.8309721046630296e-34, -8.9740761521756660e-51},
+ {5.3240312787719801e-01, -4.1020306135800895e-17, -1.9239996374230821e-33, -1.5326987913812184e-49},
+ {5.3499761988709726e-01, -5.3683132708358134e-17, -1.3900569918838112e-33, 2.7154084726474092e-50},
+ {5.3758707629564551e-01, -2.2617365388403054e-17, -5.9787279033447075e-34, 3.1204419729043625e-51},
+ {5.4017147272989285e-01, 2.7072447965935839e-17, 1.1698799709213829e-33, -5.9094668515881500e-50},
+ {5.4275078486451589e-01, 1.7148261004757101e-17, -1.3525905925200870e-33, 4.9724411290727323e-50},
+ {5.4532498842204646e-01, -4.1517817538384258e-17, -1.5318930219385941e-33, 6.3629921101413974e-50},
+ {5.4789405917310019e-01, -2.4065878297113363e-17, -3.5639213669362606e-36, -2.6013270854271645e-52},
+ {5.5045797293660481e-01, -8.3319903015807663e-18, -2.3058454035767633e-34, -2.1611290432369010e-50},
+ {5.5301670558002758e-01, -4.7061536623798204e-17, -1.0617111545918056e-33, -1.6196316144407379e-50},
+ {5.5557023301960218e-01, 4.7094109405616768e-17, -2.0640520383682921e-33, 1.2290163188567138e-49},
+ {5.5811853122055610e-01, 1.3481176324765226e-17, -5.5016743873011438e-34, -2.3484822739335416e-50},
+ {5.6066157619733603e-01, -7.3956418153476152e-18, 3.9680620611731193e-34, 3.1995952200836223e-50},
+ {5.6319934401383409e-01, 2.3835775146854829e-17, 1.3511793173769814e-34, 9.3201311581248143e-51},
+ {5.6573181078361323e-01, -3.4096079596590466e-17, -1.7073289744303546e-33, 8.9147089975404507e-50},
+ {5.6825895267013160e-01, -5.0935673642769248e-17, -1.6274356351028249e-33, 9.8183151561702966e-51},
+ {5.7078074588696726e-01, 2.4568151455566208e-17, -1.2844481247560350e-33, -1.8037634376936261e-50},
+ {5.7329716669804220e-01, 8.5176611669306400e-18, -6.4443208788026766e-34, 2.2546105543273003e-50},
+ {5.7580819141784534e-01, -3.7909495458942734e-17, -2.7433738046854309e-33, 1.1130841524216795e-49},
+ {5.7831379641165559e-01, -2.6237691512372831e-17, 1.3679051680738167e-33, -3.1409808935335900e-50},
+ {5.8081395809576453e-01, 1.8585338586613408e-17, 2.7673843114549181e-34, 1.9605349619836937e-50},
+ {5.8330865293769829e-01, 3.4516601079044858e-18, 1.8065977478946306e-34, -6.3953958038544646e-51},
+ {5.8579785745643886e-01, -3.7485501964311294e-18, 2.7965403775536614e-34, -7.1816936024157202e-51},
+ {5.8828154822264533e-01, -2.9292166725006846e-17, -2.3744954603693934e-33, -1.1571631191512480e-50},
+ {5.9075970185887428e-01, -4.7013584170659542e-17, 2.4808417611768356e-33, 1.2598907673643198e-50},
+ {5.9323229503979980e-01, 1.2892320944189053e-17, 5.3058364776359583e-34, 4.1141674699390052e-50},
+ {5.9569930449243336e-01, -1.3438641936579467e-17, -6.7877687907721049e-35, -5.6046937531684890e-51},
+ {5.9816070699634227e-01, 3.8801885783000657e-17, -1.2084165858094663e-33, -4.0456610843430061e-50},
+ {6.0061647938386897e-01, -4.6398198229461932e-17, -1.6673493003710801e-33, 5.1982824378491445e-50},
+ {6.0306659854034816e-01, 3.7323357680559650e-17, 2.7771920866974305e-33, -1.6194229649742458e-49},
+ {6.0551104140432555e-01, -3.1202672493305677e-17, 1.2761267338680916e-33, -4.0859368598379647e-50},
+ {6.0794978496777363e-01, 3.5160832362096660e-17, -2.5546242776778394e-34, -1.4085313551220694e-50},
+ {6.1038280627630948e-01, -2.2563265648229169e-17, 1.3185575011226730e-33, 8.2316691420063460e-50},
+ {6.1281008242940971e-01, -4.2693476568409685e-18, 2.5839965886650320e-34, 1.6884412005622537e-50},
+ {6.1523159058062682e-01, 2.6231417767266950e-17, -1.4095366621106716e-33, 7.2058690491304558e-50},
+ {6.1764730793780398e-01, -4.7478594510902452e-17, -7.2986558263123996e-34, -3.0152327517439154e-50},
+ {6.2005721176328921e-01, -2.7983410837681118e-17, 1.1649951056138923e-33, -5.4539089117135207e-50},
+ {6.2246127937414997e-01, 5.2940728606573002e-18, -4.8486411215945827e-35, 1.2696527641980109e-52},
+ {6.2485948814238634e-01, 3.3671846037243900e-17, -2.7846053391012096e-33, 5.6102718120012104e-50},
+ {6.2725181549514408e-01, 3.0763585181253225e-17, 2.7068930273498138e-34, -1.1172240309286484e-50},
+ {6.2963823891492698e-01, 4.1115334049626806e-17, -1.9167473580230747e-33, 1.1118424028161730e-49},
+ {6.3201873593980906e-01, -4.0164942296463612e-17, -7.2208643641736723e-34, 3.7828920470544344e-50},
+ {6.3439328416364549e-01, 1.0420901929280035e-17, 4.1174558929280492e-34, -1.4464152986630705e-51},
+ {6.3676186123628420e-01, 3.1419048711901611e-17, -2.2693738415126449e-33, -1.6023584204297388e-49},
+ {6.3912444486377573e-01, 1.2416796312271043e-17, -6.2095419626356605e-34, 2.7762065999506603e-50},
+ {6.4148101280858316e-01, -9.9883430115943310e-18, 4.1969230376730128e-34, 5.6980543799257597e-51},
+ {6.4383154288979150e-01, -3.2084798795046886e-17, -1.2595311907053305e-33, -4.0205885230841536e-50},
+ {6.4617601298331639e-01, -2.9756137382280815e-17, -1.0275370077518259e-33, 8.0852478665893014e-51},
+ {6.4851440102211244e-01, 3.9870270313386831e-18, 1.9408388509540788e-34, -5.1798420636193190e-51},
+ {6.5084668499638088e-01, 3.9714670710500257e-17, 2.9178546787002963e-34, 3.8140635508293278e-51},
+ {6.5317284295377676e-01, 8.5695642060026238e-18, -6.9165322305070633e-34, 2.3873751224185395e-50},
+ {6.5549285299961535e-01, 3.5638734426385005e-17, 1.2695365790889811e-33, 4.3984952865412050e-50},
+ {6.5780669329707864e-01, 1.9580943058468545e-17, -1.1944272256627192e-33, 2.8556402616436858e-50},
+ {6.6011434206742048e-01, -1.3960054386823638e-19, 6.1515777931494047e-36, 5.3510498875622660e-52},
+ {6.6241577759017178e-01, -2.2615508885764591e-17, 5.0177050318126862e-34, 2.9162532399530762e-50},
+ {6.6471097820334490e-01, -3.6227793598034367e-17, -9.0607934765540427e-34, 3.0917036342380213e-50},
+ {6.6699992230363747e-01, 3.5284364997428166e-17, -1.0382057232458238e-33, 7.3812756550167626e-50},
+ {6.6928258834663612e-01, -5.4592652417447913e-17, -2.5181014709695152e-33, -1.6867875999437174e-49},
+ {6.7155895484701844e-01, -4.0489037749296692e-17, 3.1995835625355681e-34, -1.4044414655670960e-50},
+ {6.7382900037875604e-01, 2.3091901236161086e-17, 5.7428037192881319e-34, 1.1240668354625977e-50},
+ {6.7609270357531592e-01, 3.7256902248049466e-17, 1.7059417895764375e-33, 9.7326347795300652e-50},
+ {6.7835004312986147e-01, 1.8302093041863122e-17, 9.5241675746813072e-34, 5.0328101116133503e-50},
+ {6.8060099779545302e-01, 2.8473293354522047e-17, 4.1331805977270903e-34, 4.2579030510748576e-50},
+ {6.8284554638524808e-01, -1.2958058061524531e-17, 1.8292386959330698e-34, 3.4536209116044487e-51},
+ {6.8508366777270036e-01, 2.5948135194645137e-17, -8.5030743129500702e-34, -6.9572086141009930e-50},
+ {6.8731534089175916e-01, -5.5156158714917168e-17, 1.1896489854266829e-33, -7.8505896218220662e-51},
+ {6.8954054473706694e-01, -1.5889323294806790e-17, 9.1242356240205712e-34, 3.8315454152267638e-50},
+ {6.9175925836415775e-01, 2.7406078472410668e-17, 1.3286508943202092e-33, 1.0651869129580079e-51},
+ {6.9397146088965400e-01, 7.4345076956280137e-18, 7.5061528388197460e-34, -1.5928000240686583e-50},
+ {6.9617713149146299e-01, -4.1224081213582889e-17, -3.1838716762083291e-35, -3.9625587412119131e-51},
+ {6.9837624940897280e-01, 4.8988282435667768e-17, 1.9134010413244152e-33, 2.6161153243793989e-50},
+ {7.0056879394324834e-01, 3.1027960192992922e-17, 9.5638250509179997e-34, 4.5896916138107048e-51},
+ {7.0275474445722530e-01, 2.5278294383629822e-18, -8.6985561210674942e-35, -5.6899862307812990e-51},
+ {7.0493408037590488e-01, 2.7608725585748502e-17, 2.9816599471629137e-34, 1.1533044185111206e-50},
+ {7.0710678118654757e-01, -4.8336466567264567e-17, 2.0693376543497068e-33, 2.4677734957341755e-50}
+};
+
+static INLINE void
+sin_table_qd(double *s0, double *s1, double *s2, double *s3, double j)
+{
+ int int_j=(int)j;
+ s0[0]=s_table[int_j-1][0];
+ s1[0]=s_table[int_j-1][1];
+ s2[0]=s_table[int_j-1][2];
+ s3[0]=s_table[int_j-1][3];
+}
+
+static double c_table[265][4] = {
+ {9.9999529380957619e-01, -1.9668064285322189e-17, -6.3053955095883481e-34, 5.3266110855726731e-52},
+ {9.9998117528260111e-01, 3.3568103522895585e-17, -1.4740132559368063e-35, 9.8603097594755596e-52},
+ {9.9995764455196390e-01, -3.1527836866647287e-17, 2.6363251186638437e-33, 1.0007504815488399e-49},
+ {9.9992470183914450e-01, 3.7931082512668012e-17, -8.5099918660501484e-35, -4.9956973223295153e-51},
+ {9.9988234745421256e-01, -3.5477814872408538e-17, 1.7102001035303974e-33, -1.0725388519026542e-49},
+ {9.9983058179582340e-01, 1.8825140517551119e-17, -5.1383513457616937e-34, -3.8378827995403787e-50},
+ {9.9976940535121528e-01, 4.2681177032289012e-17, 1.9062302359737099e-33, -6.0221153262881160e-50},
+ {9.9969881869620425e-01, -2.9851486403799753e-17, -1.9084787370733737e-33, 5.5980260344029202e-51},
+ {9.9961882249517864e-01, -4.1181965521424734e-17, 2.0915365593699916e-33, 8.1403390920903734e-50},
+ {9.9952941750109314e-01, 2.0517917823755591e-17, -4.7673802585706520e-34, -2.9443604198656772e-50},
+ {9.9943060455546173e-01, 3.9644497752257798e-17, -2.3757223716722428e-34, -1.2856759011361726e-51},
+ {9.9932238458834954e-01, -4.2858538440845682e-17, 3.3235101605146565e-34, -8.3554272377057543e-51},
+ {9.9920475861836389e-01, 9.1796317110385693e-18, 5.5416208185868570e-34, 8.0267046717615311e-52},
+ {9.9907772775264536e-01, 2.1419007653587032e-17, -7.9048203318529618e-34, -5.3166296181112712e-50},
+ {9.9894129318685687e-01, -2.0610641910058638e-17, -1.2546525485913485e-33, -7.5175888806157064e-50},
+ {9.9879545620517241e-01, -1.2291693337075465e-17, 2.4468446786491271e-34, 1.0723891085210268e-50},
+ {9.9864021818026527e-01, -4.8690254312923302e-17, -2.9470881967909147e-34, -1.3000650761346907e-50},
+ {9.9847558057329477e-01, -2.2002931182778795e-17, -1.2371509454944992e-33, -2.4911225131232065e-50},
+ {9.9830154493389289e-01, -5.1869402702792278e-17, 1.0480195493633452e-33, -2.8995649143155511e-50},
+ {9.9811811290014918e-01, 2.7935487558113833e-17, 2.4423341255830345e-33, -6.7646699175334417e-50},
+ {9.9792528619859600e-01, 1.7143659778886362e-17, 5.7885840902887460e-34, -9.2601432603894597e-51},
+ {9.9772306664419164e-01, -2.6394475274898721e-17, -1.6176223087661783e-34, -9.9924942889362281e-51},
+ {9.9751145614030345e-01, 5.6007205919806937e-18, -5.9477673514685690e-35, -1.4166807162743627e-54},
+ {9.9729045667869021e-01, 9.1647695371101735e-18, 6.7824134309739296e-34, -8.6191392795543357e-52},
+ {9.9706007033948296e-01, 1.6734093546241963e-17, -1.3169951440780028e-33, 1.0311048767952477e-50},
+ {9.9682029929116567e-01, 4.7062820708615655e-17, 2.8412041076474937e-33, -8.0006155670263622e-50},
+ {9.9657114579055484e-01, 1.1707179088390986e-17, -7.5934413263024663e-34, 2.8474848436926008e-50},
+ {9.9631261218277800e-01, 1.1336497891624735e-17, 3.4002458674414360e-34, 7.7419075921544901e-52},
+ {9.9604470090125197e-01, 2.2870031707670695e-17, -3.9184839405013148e-34, -3.7081260416246375e-50},
+ {9.9576741446765982e-01, -2.3151908323094359e-17, -1.6306512931944591e-34, -1.5925420783863192e-51},
+ {9.9548075549192694e-01, 3.2084621412226554e-18, -4.9501292146013023e-36, -2.7811428850878516e-52},
+ {9.9518472667219693e-01, -4.2486913678304410e-17, 1.3315510772504614e-33, 6.7927987417051888e-50},
+ {9.9487933079480562e-01, 4.2130813284943662e-18, -4.2062597488288452e-35, 2.5157064556087620e-51},
+ {9.9456457073425542e-01, 3.6745069641528058e-17, -3.0603304105471010e-33, 1.0397872280487526e-49},
+ {9.9424044945318790e-01, 4.4129423472462673e-17, -3.0107231708238066e-33, 7.4201582906861892e-50},
+ {9.9390697000235606e-01, -1.8964849471123746e-17, -1.5980853777937752e-35, -8.5374807150597082e-52},
+ {9.9356413552059530e-01, 2.9752309927797428e-17, -4.5066707331134233e-34, -3.3548191633805036e-50},
+ {9.9321194923479450e-01, 3.3096906261272262e-17, 1.5592811973249567e-33, 1.4373977733253592e-50},
+ {9.9285041445986510e-01, -1.4094517733693302e-17, -1.1954558131616916e-33, 1.8761873742867983e-50},
+ {9.9247953459870997e-01, 3.1093055095428906e-17, -1.8379594757818019e-33, -3.9885758559381314e-51},
+ {9.9209931314219180e-01, -3.9431926149588778e-17, -6.2758062911047230e-34, -1.2960929559212390e-50},
+ {9.9170975366909953e-01, -2.3372891311883661e-18, 2.7073298824968591e-35, -1.2569459441802872e-51},
+ {9.9131085984611544e-01, -2.5192111583372105e-17, -1.2852471567380887e-33, 5.2385212584310483e-50},
+ {9.9090263542778001e-01, 1.5394565094566704e-17, -1.0799984133184567e-33, 2.7451115960133595e-51},
+ {9.9048508425645709e-01, -5.5411437553780867e-17, -1.4614017210753585e-33, -3.8339374397387620e-50},
+ {9.9005821026229712e-01, -1.7055485906233963e-17, 1.3454939685758777e-33, 7.3117589137300036e-50},
+ {9.8962201746320089e-01, -5.2398217968132530e-17, 1.3463189211456219e-33, 5.8021640554894872e-50},
+ {9.8917650996478101e-01, -4.0987309937047111e-17, -4.4857560552048437e-34, -3.9414504502871125e-50},
+ {9.8872169196032378e-01, -1.0976227206656125e-17, 3.2311342577653764e-34, 9.6367946583575041e-51},
+ {9.8825756773074946e-01, 2.7030607784372632e-17, 7.7514866488601377e-35, 2.1019644956864938e-51},
+ {9.8778414164457218e-01, -2.3600693397159021e-17, -1.2323283769707861e-33, 3.0130900716803339e-50},
+ {9.8730141815785843e-01, -5.2332261255715652e-17, -2.7937644333152473e-33, 1.2074160567958408e-49},
+ {9.8680940181418553e-01, -5.0287214351061075e-17, -2.2681526238144461e-33, 4.4003694320169133e-50},
+ {9.8630809724459867e-01, -2.1520877103013341e-17, 1.1866528054187716e-33, -7.8532199199813836e-50},
+ {9.8579750916756748e-01, -5.1439452979953012e-17, 2.6276439309996725e-33, 7.5423552783286347e-50},
+ {9.8527764238894122e-01, 2.3155637027900207e-17, -7.5275971545764833e-34, 1.0582231660456094e-50},
+ {9.8474850180190421e-01, 1.0548144061829957e-17, 2.8786145266267306e-34, -3.6782210081466112e-51},
+ {9.8421009238692903e-01, 4.7983922627050691e-17, 2.2597419645070588e-34, 1.7573875814863400e-50},
+ {9.8366241921173025e-01, 1.9864948201635255e-17, -1.0743046281211033e-35, 1.7975662796558100e-52},
+ {9.8310548743121629e-01, 4.2170007522888628e-17, 8.2396265656440904e-34, -8.0803700139096561e-50},
+ {9.8253930228744124e-01, 1.5149580813777224e-17, -4.1802771422186237e-34, -2.2150174326226160e-50},
+ {9.8196386910955524e-01, 2.1108443711513084e-17, -1.5253013442896054e-33, -6.8388082079337969e-50},
+ {9.8137919331375456e-01, 1.3428163260355633e-17, -6.5294290469962986e-34, 2.7965412287456268e-51},
+ {9.8078528040323043e-01, 1.8546939997825006e-17, -1.0696564445530757e-33, 6.6668174475264961e-50},
+ {9.8018213596811743e-01, -3.6801786963856159e-17, 6.3245171387992842e-34, 1.8600176137175971e-50},
+ {9.7956976568544052e-01, 1.5573991584990420e-17, -1.3401066029782990e-33, -1.7263702199862149e-50},
+ {9.7894817531906220e-01, -2.3817727961148053e-18, -1.0694750370381661e-34, -8.2293047196087462e-51},
+ {9.7831737071962765e-01, -2.1623082233344895e-17, 1.0970403012028032e-33, 7.7091923099369339e-50},
+ {9.7767735782450993e-01, 5.0514136167059628e-17, -1.3254751701428788e-33, 7.0161254312124538e-50},
+ {9.7702814265775439e-01, -4.3353875751555997e-17, 5.4948839831535478e-34, -9.2755263105377306e-51},
+ {9.7636973133002114e-01, 9.3093931526213780e-18, -4.1184949155685665e-34, -3.1913926031393690e-50},
+ {9.7570213003852857e-01, -2.5572556081259686e-17, -9.3174244508942223e-34, -8.3675863211646863e-51},
+ {9.7502534506699412e-01, 2.6642660651899135e-17, 1.7819392739353853e-34, -3.3159625385648947e-51},
+ {9.7433938278557586e-01, 2.3041221476151512e-18, 1.0758686005031430e-34, 5.1074116432809478e-51},
+ {9.7364424965081198e-01, -5.1729808691005871e-17, -1.5508473005989887e-33, -1.6505125917675401e-49},
+ {9.7293995220556018e-01, -3.1311211122281800e-17, -2.6874087789006141e-33, -2.1652434818822145e-51},
+ {9.7222649707893627e-01, 3.6461169785938221e-17, 3.0309636883883133e-33, -1.2702716907967306e-51},
+ {9.7150389098625178e-01, -7.9865421122289046e-18, -4.3628417211263380e-34, 3.4307517798759352e-51},
+ {9.7077214072895035e-01, -4.7992163325114922e-17, 3.0347528910975783e-33, 8.5989199506479701e-50},
+ {9.7003125319454397e-01, 1.8365300348428844e-17, -1.4311097571944918e-33, 8.5846781998740697e-51},
+ {9.6928123535654853e-01, -4.5663660261927896e-17, 9.6147526917239387e-34, 8.1267605207871330e-51},
+ {9.6852209427441727e-01, 4.9475074918244771e-17, 2.8558738351911241e-33, 6.2948422316507461e-50},
+ {9.6775383709347551e-01, -4.5512132825515820e-17, -1.4127617988719093e-33, -8.4620609089704578e-50},
+ {9.6697647104485207e-01, 3.8496228837337864e-17, -5.3881631542745647e-34, -3.5221863171458959e-50},
+ {9.6619000344541250e-01, 5.1298840401665493e-17, 1.4564075904769808e-34, 1.0095973971377432e-50},
+ {9.6539444169768940e-01, -2.3745389918392156e-17, 5.9221515590053862e-34, -3.8811192556231094e-50},
+ {9.6458979328981276e-01, -3.4189470735959786e-17, 2.2982074155463522e-33, -4.5128791045607634e-50},
+ {9.6377606579543984e-01, 2.6463950561220029e-17, -2.9073234590199323e-36, -1.2938328629395601e-52},
+ {9.6295326687368388e-01, 8.9341960404313634e-18, -3.9071244661020126e-34, 1.6212091116847394e-50},
+ {9.6212140426904158e-01, 1.5236770453846305e-17, -1.3050173525597142e-33, 7.9016122394092666e-50},
+ {9.6128048581132064e-01, 2.0933955216674039e-18, 1.0768607469015692e-34, -5.9453639304361774e-51},
+ {9.6043051941556579e-01, 2.4653904815317185e-17, -1.3792169410906322e-33, -4.7726598378506903e-51},
+ {9.5957151308198452e-01, 1.1000640085000957e-17, -4.2036030828223975e-34, 4.0023704842606573e-51},
+ {9.5870347489587160e-01, -4.3685014392372053e-17, 2.2001800662729131e-33, -1.0553721324358075e-49},
+ {9.5782641302753291e-01, -1.7696710075371263e-17, 1.9164034110382190e-34, 8.1489235071754813e-51},
+ {9.5694033573220882e-01, 4.0553869861875701e-17, -1.7147013364302149e-33, 2.5736745295329455e-50},
+ {9.5604525134999641e-01, 3.7705045279589067e-17, 1.9678699997347571e-33, 8.5093177731230180e-50},
+ {9.5514116830577067e-01, 5.0088652955014668e-17, -2.6983181838059211e-33, 1.0102323575596493e-49},
+ {9.5422809510910567e-01, -3.7545901690626874e-17, 1.4951619241257764e-33, -8.2717333151394973e-50},
+ {9.5330604035419386e-01, -2.5190738779919934e-17, -1.4272239821134379e-33, -4.6717286809283155e-50},
+ {9.5237501271976588e-01, -2.0269300462299272e-17, -1.0635956887246246e-33, -3.5514537666487619e-50},
+ {9.5143502096900834e-01, 3.1350584123266695e-17, -2.4824833452737813e-33, 9.5450335525380613e-51},
+ {9.5048607394948170e-01, 1.9410097562630436e-17, -8.1559393949816789e-34, -1.0501209720164562e-50},
+ {9.4952818059303667e-01, -7.5544151928043298e-18, -5.1260245024046686e-34, 1.8093643389040406e-50},
+ {9.4856134991573027e-01, 2.0668262262333232e-17, -5.9440730243667306e-34, 1.4268853111554300e-50},
+ {9.4758559101774109e-01, 4.3417993852125991e-17, -2.7728667889840373e-34, 5.5709160196519968e-51},
+ {9.4660091308328353e-01, 3.5056800210680730e-17, 9.8578536940318117e-34, 6.6035911064585197e-50},
+ {9.4560732538052128e-01, 4.6019102478523738e-17, -6.2534384769452059e-34, 1.5758941215779961e-50},
+ {9.4460483726148026e-01, 8.8100545476641165e-18, 5.2291695602757842e-34, -3.3487256018407123e-50},
+ {9.4359345816196039e-01, -2.4093127844404214e-17, 1.0283279856803939e-34, -2.3398232614531355e-51},
+ {9.4257319760144687e-01, 1.3235564806436886e-17, -5.7048262885386911e-35, 3.9947050442753744e-51},
+ {9.4154406518302081e-01, -2.7896379547698341e-17, 1.6273236356733898e-33, -5.3075944708471203e-51},
+ {9.4050607059326830e-01, 2.8610421567116268e-17, 2.9261501147538827e-33, -2.6849867690896925e-50},
+ {9.3945922360218992e-01, -7.0152867943098655e-18, -5.6395693818011210e-34, 3.5568142678987651e-50},
+ {9.3840353406310806e-01, 5.4242545044795490e-17, -1.9039966607859759e-33, -1.5627792988341215e-49},
+ {9.3733901191257496e-01, -3.6570926284362776e-17, -1.1902940071273247e-33, -1.1215082331583223e-50},
+ {9.3626566717027826e-01, -1.3013766145497654e-17, 5.2229870061990595e-34, -3.3972777075634108e-51},
+ {9.3518350993894761e-01, -3.2609395302485065e-17, -8.1813015218875245e-34, 5.5642140024928139e-50},
+ {9.3409255040425887e-01, 4.4662824360767511e-17, -2.5903243047396916e-33, 8.1505209004343043e-50},
+ {9.3299279883473885e-01, 4.2041415555384355e-17, 9.0285896495521276e-34, 5.3019984977661259e-50},
+ {9.3188426558166815e-01, -4.0785944377318095e-17, 1.7631450298754169e-33, 2.5776403305507453e-50},
+ {9.3076696107898371e-01, 1.9703775102838329e-17, 6.5657908718278205e-34, -1.9480347966259524e-51},
+ {9.2964089584318121e-01, 5.1282530016864107e-17, 2.3719739891916261e-34, -1.7230065426917127e-50},
+ {9.2850608047321559e-01, -2.3306639848485943e-17, -7.7799084333208503e-34, -5.8597558009300305e-50},
+ {9.2736252565040111e-01, -2.7677111692155437e-17, 2.2110293450199576e-34, 2.0349190819680613e-50},
+ {9.2621024213831138e-01, -3.7303754586099054e-17, 2.0464457809993405e-33, 1.3831799631231817e-49},
+ {9.2504924078267758e-01, 6.0529447412576159e-18, -8.8256517760278541e-35, 1.8285462122388328e-51},
+ {9.2387953251128674e-01, 1.7645047084336677e-17, -5.0442537321586818e-34, -4.0478677716823890e-50},
+ {9.2270112833387852e-01, 5.2963798918539814e-17, -5.7135699628876685e-34, 3.0163671797219087e-50},
+ {9.2151403934204190e-01, 4.1639843390684644e-17, 1.1891485604702356e-33, 2.0862437594380324e-50},
+ {9.2031827670911059e-01, -2.7806888779036837e-17, 2.7011013677071274e-33, 1.1998578792455499e-49},
+ {9.1911385169005777e-01, -2.6496484622344718e-17, 6.5403604763461920e-34, -2.8997180201186078e-50},
+ {9.1790077562139050e-01, -3.9074579680849515e-17, 2.3004636541490264e-33, 3.9851762744443107e-50},
+ {9.1667905992104270e-01, -4.1733978698287568e-17, 1.2094444804381172e-33, 4.9356916826097816e-50},
+ {9.1544871608826783e-01, -1.3591056692900894e-17, 5.9923027475594735e-34, 2.1403295925962879e-50},
+ {9.1420975570353069e-01, -3.6316182527814423e-17, -1.9438819777122554e-33, 2.8340679287728316e-50},
+ {9.1296219042839821e-01, -4.7932505228039469e-17, -1.7753551889428638e-33, 4.0607782903868160e-51},
+ {9.1170603200542988e-01, -2.6913273175034130e-17, -5.1928101916162528e-35, 1.1338175936090630e-51},
+ {9.1044129225806725e-01, -5.0433041673313820e-17, 1.0938746257404305e-33, 9.5378272084170731e-51},
+ {9.0916798309052238e-01, -3.6878564091359894e-18, 2.9951330310507693e-34, -1.2225666136919926e-50},
+ {9.0788611648766626e-01, -4.9459964301225840e-17, -1.6599682707075313e-33, -5.1925202712634716e-50},
+ {9.0659570451491533e-01, 3.0506718955442023e-17, -1.4478836557141204e-33, 1.8906373784448725e-50},
+ {9.0529675931811882e-01, -4.1153099826889901e-17, 2.9859368705184223e-33, 5.1145293917439211e-50},
+ {9.0398929312344334e-01, -6.6097544687484308e-18, 1.2728013034680357e-34, -4.3026097234014823e-51},
+ {9.0267331823725883e-01, -1.9250787033961483e-17, 1.3242128993244527e-33, -5.2971030688703665e-50},
+ {9.0134884704602203e-01, -1.3524789367698682e-17, 6.3605353115880091e-34, 3.6227400654573828e-50},
+ {9.0001589201616028e-01, -5.0639618050802273e-17, 1.0783525384031576e-33, 2.8130016326515111e-50},
+ {8.9867446569395382e-01, 2.6316906461033013e-17, 3.7003137047796840e-35, -2.3447719900465938e-51},
+ {8.9732458070541832e-01, -3.6396283314867290e-17, -2.3611649895474815e-33, 1.1837247047900082e-49},
+ {8.9596624975618511e-01, 4.9025099114811813e-17, -1.9440489814795326e-33, -1.7070486667767033e-49},
+ {8.9459948563138270e-01, -1.7516226396814919e-17, -1.3200670047246923e-33, -1.5953009884324695e-50},
+ {8.9322430119551532e-01, -4.1161239151908913e-18, 2.5380253805715999e-34, 4.2849455510516192e-51},
+ {8.9184070939234272e-01, 4.6690228137124547e-18, 1.6150254286841982e-34, -3.9617448820725012e-51},
+ {8.9044872324475788e-01, 1.1781931459051803e-17, -1.3346142209571930e-34, -9.4982373530733431e-51},
+ {8.8904835585466457e-01, -1.1164514966766675e-17, -3.4797636107798736e-34, -1.5605079997040631e-50},
+ {8.8763962040285393e-01, 1.2805091918587960e-17, 3.9948742059584459e-35, 3.8940716325338136e-51},
+ {8.8622253014888064e-01, -6.7307369600274315e-18, 1.2385593432917413e-34, 2.0364014759133320e-51},
+ {8.8479709843093779e-01, -9.4331469628972690e-18, -5.7106541478701439e-34, 1.8260134111907397e-50},
+ {8.8336333866573158e-01, 1.5822643380255127e-17, -7.8921320007588250e-34, -1.4782321016179836e-50},
+ {8.8192126434835505e-01, -1.9843248405890562e-17, -7.0412114007673834e-34, -1.0636770169389104e-50},
+ {8.8047088905216075e-01, 1.6311096602996350e-17, -5.7541360594724172e-34, -4.0128611862170021e-50},
+ {8.7901222642863353e-01, -4.7356837291118011e-17, 1.4388771297975192e-33, -2.9085554304479134e-50},
+ {8.7754529020726124e-01, 5.0113311846499550e-17, 2.8382769008739543e-34, 1.5550640393164140e-50},
+ {8.7607009419540660e-01, 5.8729024235147677e-18, 2.7941144391738458e-34, -1.8536073846509828e-50},
+ {8.7458665227817611e-01, -5.7216617730397065e-19, -2.9705811503689596e-35, 8.7389593969796752e-52},
+ {8.7309497841829009e-01, 7.8424672990129903e-18, -4.8685015839797165e-34, -2.2815570587477527e-50},
+ {8.7159508665595109e-01, -5.5272998038551050e-17, -2.2104090204984907e-33, -9.7749763187643172e-50},
+ {8.7008699110871146e-01, -4.1888510868549968e-17, 7.0900185861878415e-34, 3.7600251115157260e-50},
+ {8.6857070597134090e-01, 2.7192781689782903e-19, -1.6710140396932428e-35, -1.2625514734637969e-51},
+ {8.6704624551569265e-01, 3.0267859550930567e-18, -1.1559438782171572e-34, -5.3580556397808012e-52},
+ {8.6551362409056909e-01, -6.3723113549628899e-18, 2.3725520321746832e-34, 1.5911880348395175e-50},
+ {8.6397285612158670e-01, 4.1486355957361607e-17, 2.2709976932210266e-33, -8.1228385659479984e-50},
+ {8.6242395611104050e-01, 3.7008992527383130e-17, 5.2128411542701573e-34, 2.6945600081026861e-50},
+ {8.6086693863776731e-01, -3.0050048898573656e-17, -8.8706183090892111e-34, 1.5005320558097301e-50},
+ {8.5930181835700836e-01, 4.2435655816850687e-17, 7.6181814059912025e-34, -3.9592127850658708e-50},
+ {8.5772861000027212e-01, -4.8183447936336620e-17, -1.1044130517687532e-33, -8.7400233444645562e-50},
+ {8.5614732837519447e-01, 9.1806925616606261e-18, 5.6328649785951470e-34, 2.3326646113217378e-51},
+ {8.5455798836540053e-01, -1.2991124236396092e-17, 1.2893407722948080e-34, -3.6506925747583053e-52},
+ {8.5296060493036363e-01, 2.7152984251981370e-17, 7.4336483283120719e-34, 4.2162417622350668e-50},
+ {8.5135519310526520e-01, -5.3279874446016209e-17, 2.2281156380919942e-33, -4.0281886404138477e-50},
+ {8.4974176800085244e-01, 5.1812347659974015e-17, 3.0810626087331275e-33, -2.5931308201994965e-50},
+ {8.4812034480329723e-01, 1.8762563415239981e-17, 1.4048773307919617e-33, -2.4915221509958691e-50},
+ {8.4649093877405213e-01, -4.7969419958569345e-17, -2.7518267097886703e-33, -7.3518959727313350e-50},
+ {8.4485356524970712e-01, -4.3631360296879637e-17, -2.0307726853367547e-33, 4.3097229819851761e-50},
+ {8.4320823964184544e-01, 9.6536707005959077e-19, 2.8995142431556364e-36, 9.6715076811480284e-53},
+ {8.4155497743689844e-01, -3.4095465391321557e-17, -8.4130208607579595e-34, -4.9447283960568686e-50},
+ {8.3989379419599952e-01, -1.6673694881511411e-17, -1.4759184141750289e-33, -7.5795098161914058e-50},
+ {8.3822470555483808e-01, -3.5560085052855026e-17, 1.1689791577022643e-33, -5.8627347359723411e-50},
+ {8.3654772722351201e-01, -2.0899059027066533e-17, -9.8104097821002585e-35, -3.1609177868229853e-51},
+ {8.3486287498638001e-01, 4.6048430609159657e-17, -5.1827423265239912e-34, -7.0505343435504109e-51},
+ {8.3317016470191319e-01, 1.3275129507229764e-18, 4.8589164115370863e-35, 4.5422281300506859e-51},
+ {8.3146961230254524e-01, 1.4073856984728024e-18, 4.6951315383980830e-35, 5.1431906049905658e-51},
+ {8.2976123379452305e-01, -2.9349109376485597e-18, 1.1496917934149818e-34, 3.5186665544980233e-51},
+ {8.2804504525775580e-01, -4.4196593225871532e-17, 2.7967864855211251e-33, 1.0030777287393502e-49},
+ {8.2632106284566353e-01, -5.3957485453612902e-17, 6.8976896130138550e-34, 3.8106164274199196e-50},
+ {8.2458930278502529e-01, -2.6512360488868275e-17, 1.6916964350914386e-34, 6.7693974813562649e-51},
+ {8.2284978137582632e-01, 1.5193019034505495e-17, 9.6890547246521685e-34, 5.6994562923653264e-50},
+ {8.2110251499110465e-01, 3.0715131609697682e-17, -1.7037168325855879e-33, -1.1149862443283853e-49},
+ {8.1934752007679701e-01, -4.8200736995191133e-17, -1.5574489646672781e-35, -9.5647853614522216e-53},
+ {8.1758481315158371e-01, -1.4883149812426772e-17, -7.8273262771298917e-34, 4.1332149161031594e-50},
+ {8.1581441080673378e-01, 8.2652693782130871e-18, -2.3028778135179471e-34, 1.5102071387249843e-50},
+ {8.1403632970594841e-01, -5.2127351877042624e-17, -1.9047670611316360e-33, -1.6937269585941507e-49},
+ {8.1225058658520388e-01, 3.1054545609214803e-17, 2.2649541922707251e-34, -7.4221684154649405e-51},
+ {8.1045719825259477e-01, 2.3520367349840499e-17, -7.7530070904846341e-34, -7.2792616357197140e-50},
+ {8.0865618158817498e-01, 9.3251597879721674e-18, -7.1823301933068394e-34, 2.3925440846132106e-50},
+ {8.0684755354379922e-01, 4.9220603766095546e-17, 2.9796016899903487e-33, 1.5220754223615788e-49},
+ {8.0503133114296355e-01, 5.1368289568212149e-17, 6.3082807402256524e-34, 7.3277646085129827e-51},
+ {8.0320753148064494e-01, -3.3060609804814910e-17, -1.2242726252420433e-33, 2.8413673268630117e-50},
+ {8.0137617172314024e-01, -2.0958013413495834e-17, -4.3798162198006931e-34, 2.0235690497752515e-50},
+ {7.9953726910790501e-01, 2.0356723822005431e-17, -9.7448513696896360e-34, 5.3608109599696008e-52},
+ {7.9769084094339116e-01, -4.6730759884788944e-17, 2.3075897077191757e-33, 3.1605567774640253e-51},
+ {7.9583690460888357e-01, -3.0062724851910721e-17, -2.2496210832042235e-33, -6.5881774117183040e-50},
+ {7.9397547755433717e-01, -7.4194631759921416e-18, 2.4124341304631069e-34, -4.9956808616244972e-51},
+ {7.9210657730021239e-01, -3.7087850202326467e-17, -1.4874457267228264e-33, 2.9323097289153505e-50},
+ {7.9023022143731003e-01, 2.3056905954954492e-17, 1.4481080533260193e-33, -7.6725237057203488e-50},
+ {7.8834642762660623e-01, 3.4396993154059708e-17, 1.7710623746737170e-33, 1.7084159098417402e-49},
+ {7.8645521359908577e-01, -9.7841429939305265e-18, 3.3906063272445472e-34, 5.7269505320382577e-51},
+ {7.8455659715557524e-01, -8.5627965423173476e-18, -2.1106834459001849e-34, -1.6890322182469603e-50},
+ {7.8265059616657573e-01, 9.0745866975808825e-18, 6.7623847404278666e-34, -1.7173237731987271e-50},
+ {7.8073722857209449e-01, -9.9198782066678806e-18, -2.1265794012162715e-36, 3.0772165598957647e-54},
+ {7.7881651238147598e-01, -2.4891385579973807e-17, 6.7665497024807980e-35, -6.5218594281701332e-52},
+ {7.7688846567323244e-01, 7.7418602570672864e-18, -5.9986517872157897e-34, 3.0566548232958972e-50},
+ {7.7495310659487393e-01, -5.2209083189826433e-17, -9.6653593393686612e-34, 3.7027750076562569e-50},
+ {7.7301045336273699e-01, -3.2565907033649772e-17, 1.3860807251523929e-33, -3.9971329917586022e-50},
+ {7.7106052426181382e-01, -4.4558442347769265e-17, -2.9863565614083783e-33, -6.8795262083596236e-50},
+ {7.6910333764557959e-01, 5.1546455184564817e-17, 2.6142829553524292e-33, -1.6199023632773298e-49},
+ {7.6713891193582040e-01, -1.8885903683750782e-17, -1.3659359331495433e-33, -2.2538834962921934e-50},
+ {7.6516726562245896e-01, -3.2707225612534598e-17, 1.1177117747079528e-33, -3.7005182280175715e-50},
+ {7.6318841726338127e-01, 2.6314748416750748e-18, 1.4048039063095910e-34, 8.9601886626630321e-52},
+ {7.6120238548426178e-01, 3.5315510881690551e-17, 1.2833566381864357e-33, 8.6221435180890613e-50},
+ {7.5920918897838807e-01, -3.8558842175523123e-17, 2.9720241208332759e-34, -1.2521388928220163e-50},
+ {7.5720884650648457e-01, -1.9909098777335502e-17, 3.9409283266158482e-34, 2.0744254207802976e-50},
+ {7.5520137689653655e-01, -1.9402238001823017e-17, -3.7756206444727573e-34, -2.1212242308178287e-50},
+ {7.5318679904361252e-01, -3.7937789838736540e-17, -6.7009539920231559e-34, -6.7128562115050214e-51},
+ {7.5116513190968637e-01, 4.3499761158645868e-17, 2.5227718971102212e-33, -6.5969709212757102e-50},
+ {7.4913639452345937e-01, -4.4729078447011889e-17, -2.4206025249983768e-33, 1.1336681351116422e-49},
+ {7.4710060598018013e-01, 1.1874824875965430e-17, 2.1992523849833518e-34, 1.1025018564644483e-50},
+ {7.4505778544146595e-01, 1.5078686911877863e-17, 8.0898987212942471e-34, 8.2677958765323532e-50},
+ {7.4300795213512172e-01, -2.5144629669719265e-17, 7.1128989512526157e-34, 3.0181629077821220e-50},
+ {7.4095112535495911e-01, -1.4708616952297345e-17, -4.9550433827142032e-34, 3.1434132533735671e-50},
+ {7.3888732446061511e-01, 3.4324874808225091e-17, -1.3706639444717610e-33, -3.3520827530718938e-51},
+ {7.3681656887736990e-01, -2.8932468101656295e-17, -3.4649887126202378e-34, -1.8484474476291476e-50},
+ {7.3473887809596350e-01, -3.4507595976263941e-17, -2.3718000676666409e-33, -3.9696090387165402e-50},
+ {7.3265427167241282e-01, 1.8918673481573520e-17, -1.5123719544119886e-33, -9.7922152011625728e-51},
+ {7.3056276922782759e-01, -2.9689959904476928e-17, -1.1276871244239744e-33, -3.0531520961539007e-50},
+ {7.2846439044822520e-01, 1.1924642323370718e-19, 5.9001892316611011e-36, 1.2178089069502704e-52},
+ {7.2635915508434601e-01, -3.1917502443460542e-17, 7.7047912412039396e-34, 4.1455880160182123e-50},
+ {7.2424708295146689e-01, 2.9198471334403004e-17, 2.3027324968739464e-33, -1.2928820533892183e-51},
+ {7.2212819392921535e-01, -2.3871262053452047e-17, 1.0636125432862273e-33, -4.4598638837802517e-50},
+ {7.2000250796138165e-01, -2.5689658854462333e-17, -9.1492566948567925e-34, 4.4403780801267786e-50},
+ {7.1787004505573171e-01, 2.7006476062511453e-17, -2.2854956580215348e-34, 9.1726903890287867e-51},
+ {7.1573082528381871e-01, -5.1581018476410262e-17, -1.3736271349300259e-34, -1.2734611344111297e-50},
+ {7.1358486878079364e-01, -4.2342504403133584e-17, -4.2690366101617268e-34, -2.6352370883066522e-50},
+ {7.1143219574521643e-01, 7.9643298613856813e-18, 2.9488239510721469e-34, 1.6985236437666356e-50},
+ {7.0927282643886569e-01, -3.7597359110245730e-17, 1.0613125954645119e-34, 8.9465480185486032e-51},
+ {7.0710678118654757e-01, -4.8336466567264567e-17, 2.0693376543497068e-33, 2.4677734957341755e-50}
+};
+
+static INLINE void
+cos_table_qd(double *c0, double *c1, double *c2, double *c3, double j)
+{
+ int int_j=(int)j;
+ c0[0]=c_table[int_j-1][0];
+ c1[0]=c_table[int_j-1][1];
+ c2[0]=c_table[int_j-1][2];
+ c3[0]=c_table[int_j-1][3];
+ 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},
+ {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}
+};
+
+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
+ double thresh = 0.5*eps;
+ double y0,y1,y2,y3,r0,r1,r2,r3,t0,t1,t2,t3,p0,p1,p2,p3;
+ int i;
+
+ if(x0==0.0) {
+ c0[0]=1.0; c1[0]=0.0; c2[0]=0.0; c3[0]=0.0;
+ return;
+ }
+
+ i=1;
+ qd_mul_qd(&y0,&y1,&y2,&y3,x0,x1,x2,x3,x0,x1,x2,x3);
+ y0 = -y0; y1 = -y1; y2 = -y2; y3 = -y3;
+ r0=y0; r1=y1; r2=y2; r3=y3;
+ s_mul_qd(&p0,&p1,&p2,&p3,0.5,r0,r1,r2,r3);
+ qd_add_s(&c0[0],&c1[0],&c2[0],&c3[0],p0,p1,p2,p3,1.0);
+
+ 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(&c0[0],&c1[0],&c2[0],&c3[0],c0[0],c1[0],c2[0],c3[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,r2,inv_fact[i][0],inv_fact[i][1],inv_fact[i][2],inv_fact[i][3]);
+ qd_add_qd(&c0[0],&c1[0],&c2[0],&c3[0],c0[0],c1[0],c2[0],c3[0],t0,t1,t2,t3);
+ i=i+2;
+ }
+}
+
+static void
+sincos_taylor_qd(double *s0, double *s1, double *s2, double *s3, double *c0, double *c1, double *c2, double *c3, 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,p0,p1,p2,p3,q0,q1,q2,q3;
+ int i;
+
+ if(x0==0.0) {
+ s0[0]=0.0; s1[0]=0.0; s2[0]=0.0; s3[0]=0.0;
+ c0[0]=1.0; c1[0]=0.0; c2[0]=0.0; c3[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)||((int)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;
+ }
+ 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);
+}
+
+//--------------------------------------------------------------------------------------------
+//
+// quad-double sine
+//
+// args
+// a0, a1, a2, a3 : double numbers
+// a0 + a1 + a2 + a3 = qd number
+//
+// return (s0,s1,s2,s3) for qd number
+//
+//--------------------------------------------------------------------------------------------
+static void
+qdsin(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
+ double _pi2[4] = {1.570796326794896558e+00,6.123233995736766036e-17,-1.497384904859169833e-33,5.562271104316826408e-50}; // pi/2
+ double _pi1024[4] = {3.067961575771282340e-03,1.195944139792337116e-19,-2.924579892303066080e-36,1.086381075061880158e-52};
+ double u0,u1,u2,u3,v0,v1,v2,v3;
+ double sin0,sin1,sin2,sin3,cos0,cos1,cos2,cos3;
+ int int_j;
+
+ if(a0[0]==0) {
+ s0[0]=0.0; s1[0]=0.0; s2[0]=0.0; s3[0]=0.0;
+ return;
+ }
+
+ //approximately reduce modulo 2*pi
+ qd_div_qd(&p0,&p1,&p2,&p3,a0[0],a1[0],a2[0],a3[0],_2pi[0],_2pi[1],_2pi[2],_2pi[3]);
+ nint_qd(&z0,&z1,&z2,&z3,p0,p1,p2,p3);
+ qd_mul_qd(&q0,&q1,&q2,&q3,_2pi[0],_2pi[1],_2pi[2],_2pi[3],z0,z1,z2,z3);
+ qd_sub_qd(&r0,&r1,&r2,&r3,a0[0],a1[0],a2[0],a3[0],q0,q1,q2,q3);
+
+ //approximately reduce modulo pi/2 and then modulo pi/1024
+ j=floor(r0/_pi2[0]+0.5);
+ s_mul_qd(&p0, &p1, &p2, &p3, j, _pi2[0], _pi2[1], _pi2[2], _pi2[3]); // Modified,2012/01/16 Saito
+ qd_sub_qd(&t0,&t1,&t2,&t3,r0,r1,r2,r3,p0,p1,p2,p3);
+ k=floor(t0/_pi1024[0]+0.5);
+ s_mul_qd(&q0, &q1, &q2, &q3, k, _pi1024[0], _pi1024[1], _pi1024[2], _pi1024[3]); // Modified,2012/01/16 Saito
+ qd_sub_qd(&t0,&t1,&t2,&t3,t0,t1,t2,t3,q0,q1,q2,q3);
+ abs_k=(int)fabs(k);
+ int_j=(int)j;
+
+ //checking errors
+ if(j<-2 || j>2) {
+ s0[0]=0.0; s1[0]=0.0; s2[0]=0.0; s3[0]=1.0;
+ return;
+ }
+
+ if(abs_k >256) {
+ s0[0]=0.0; s1[0]=0.0; s2[0]=0.0; s3[0]=1.0;
+ return;
+ }
+
+ if(k==0) {
+ switch(int_j) {
+ case 0:
+ sin_taylor_qd(&s0[0],&s1[0],&s2[0],&s3[0],t0,t1,t2,t3);
+ return;
+ case 1:
+ cos_taylor_qd(&s0[0],&s1[0],&s2[0],&s3[0],t0,t1,t2,t3);
+ return;
+ case -1:
+ cos_taylor_qd(&s0[0],&s1[0],&s2[0],&s3[0],t0,t1,t2,t3);
+ s0[0]=-s0[0]; s1[0]=-s1[0]; s2[0]=-s2[0]; s3[0]=-s3[0];
+ return;
+ case 2:
+ case -2:
+ sin_taylor_qd(&s0[0],&s1[0],&s2[0],&s3[0],t0,t1,t2,t3);
+ s0[0]=-s0[0]; s1[0]=-s1[0]; s2[0]=-s2[0]; s3[0]=-s3[0];
+ return;
+ }
+ }
+
+ cos_table_qd(&u0,&u1,&u2,&u3,abs_k);
+ sin_table_qd(&v0,&v1,&v2,&v3,abs_k);
+ sincos_taylor_qd(&sin0,&sin1,&sin2,&sin3,&cos0,&cos1,&cos2,&cos3,t0,t1,t2,t3);
+
+ if(j==0) {
+ if(k>0) {
+ qd_mul_qd(&p0,&p1,&p2,&p3,u0,u1,u2,u3,sin0,sin1,sin2,sin3);
+ qd_mul_qd(&q0,&q1,&q2,&q3,v0,v1,v2,v3,cos0,cos1,cos2,cos3);
+ qd_add_qd(&s0[0],&s1[0],&s2[0],&s3[0],p0,p1,p2,p3,q0,q1,q2,q3);
+ }
+ else {
+ qd_mul_qd(&p0,&p1,&p2,&p3,u0,u1,u2,u3,sin0,sin1,sin2,sin3);
+ qd_mul_qd(&q0,&q1,&q2,&q3,v0,v1,v2,v3,cos0,cos1,cos2,cos3);
+ qd_sub_qd(&s0[0],&s1[0],&s2[0],&s3[0],p0,p1,p2,p3,q0,q1,q2,q3);
+ }
+ }
+ else if(j==1) {
+ if(k>0) {
+ qd_mul_qd(&p0,&p1,&p2,&p3,u0,u1,u2,u3,cos0,cos1,cos2,cos3);
+ qd_mul_qd(&q0,&q1,&q2,&q3,v0,v1,v2,v3,sin0,sin1,sin2,sin3);
+ qd_sub_qd(&s0[0],&s1[0],&s2[0],&s3[0],p0,p1,p2,p3,q0,q1,q2,q3);
+ }
+ else {
+ qd_mul_qd(&p0,&p1,&p2,&p3,u0,u1,u2,u3,cos0,cos1,cos2,cos3);
+ qd_mul_qd(&q0,&q1,&q2,&q3,v0,v1,v2,v3,sin0,sin1,sin2,sin3);
+ qd_add_qd(&s0[0],&s1[0],&s2[0],&s3[0],p0,p1,p2,p3,q0,q1,q2,q3);
+ }
+ }
+ else if(j==-1) {
+ if(k>0) {
+ qd_mul_qd(&p0,&p1,&p2,&p3,v0,v1,v2,v3,sin0,sin1,sin2,sin3);
+ qd_mul_qd(&q0,&q1,&q2,&q3,u0,u1,u2,u3,cos0,cos1,cos2,cos3);
+ qd_sub_qd(&s0[0],&s1[0],&s2[0],&s3[0],p0,p1,p2,p3,q0,q1,q2,q3);
+ }
+ else {
+ qd_mul_qd(&p0,&p1,&p2,&p3,u0,u1,u2,u3,cos0,cos1,cos2,cos3);
+ qd_mul_qd(&q0,&q1,&q2,&q3,v0,v1,v2,v3,sin0,sin1,sin2,sin3);
+ p0=-p0; p1=-p1; p2=-p2; p3=-p3;
+ qd_sub_qd(&s0[0],&s1[0],&s2[0],&s3[0],p0,p1,p2,p3,q0,q1,q2,q3);
+ }
+ }
+ else {
+ if(k>0) {
+ qd_mul_qd(&p0,&p1,&p2,&p3,u0,u1,u2,u3,sin0,sin1,sin2,sin3);
+ qd_mul_qd(&q0,&q1,&q2,&q3,v0,v1,v2,v3,cos0,cos1,cos2,cos3);
+ p0=-p0; p1=-p1; p2=-p2; p3=-p3;
+ qd_sub_qd(&s0[0],&s1[0],&s2[0],&s3[0],p0,p1,p2,p3,q0,q1,q2,q3);
+ }
+ else {
+ qd_mul_qd(&p0,&p1,&p2,&p3,v0,v1,v2,v3,cos0,cos1,cos2,cos3);
+ qd_mul_qd(&q0,&q1,&q2,&q3,u0,u1,u2,u3,sin0,sin1,sin2,sin3);
+ qd_sub_qd(&s0[0],&s1[0],&s2[0],&s3[0],p0,p1,p2,p3,q0,q1,q2,q3);
+ }
+ }
+}
+
+
+//--------------------------------------------------------------------------------------------
+//
+// quad-double cosine
+//
+// args
+// a0, a1, a2, a3 : double numbers
+// a0 + a1 + a2 + a3 = qd number
+//
+// return (c0,c1,c2,c3) for qd number
+//
+//--------------------------------------------------------------------------------------------
+static void
+qdcos(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
+ double _pi2[4] = {1.570796326794896558e+00,6.123233995736766036e-17,-1.497384904859169833e-33,5.562271104316826408e-50}; // pi/2
+ double _pi1024[4] = {3.067961575771282340e-03,1.195944139792337116e-19,-2.924579892303066080e-36,1.086381075061880158e-52};
+ double u0,u1,u2,u3,v0,v1,v2,v3;
+ double sin0,sin1,sin2,sin3,cos0,cos1,cos2,cos3;
+ int int_j;
+
+ if(a0[0]==0) {
+ c1[0]=1.0; c1[0]=0.0; c2[0]=0.0; c3[0]=0.0;
+ return;
+ }
+
+ //approximately reduce modulo 2*pi
+ qd_div_qd(&p0,&p1,&p2,&p3,a0[0],a1[0],a2[0],a3[0],_2pi[0],_2pi[1],_2pi[2],_2pi[3]);
+ nint_qd(&z0,&z1,&z2,&z3,p0,p1,p2,p3);
+ qd_mul_qd(&q0,&q1,&q2,&q3,_2pi[0],_2pi[1],_2pi[2],_2pi[3],z0,z1,z2,z3);
+ qd_sub_qd(&r0,&r1,&r2,&r3,a0[0],a1[0],a2[0],a3[0],q0,q1,q2,q3);
+
+ //approximately reduce modulo pi/2 and then modulo pi/1024
+ j=floor(r0/_pi2[0]+0.5);
+ s_mul_qd(&p0, &p1, &p2, &p3, j, _pi2[0], _pi2[1], _pi2[2], _pi2[3]); // Modified,2012/01/16 Saito
+ qd_sub_qd(&t0,&t1,&t2,&t3,r0,r1,r2,r3,p0,p1,p2,p3);
+ k=floor(t0/_pi1024[0]+0.5);
+ s_mul_qd(&q0, &q1, &q2, &q3, k, _pi1024[0], _pi1024[1], _pi1024[2], _pi1024[3]); // Modified,2012/01/16 Saito
+ qd_sub_qd(&t0,&t1,&t2,&t3,t0,t1,t2,t3,q0,q1,q2,q3);
+ abs_k=(int)fabs(k);
+ int_j=(int)j;
+
+ //checking errors
+ if(j<-2 || j>2) {
+ c0[0]=0.0; c1[0]=0.0; c2[0]=0.0; c3[0]=1.0;
+ return;
+ }
+
+ if(abs_k >256) {
+ c0[0]=0.0; c1[0]=0.0; c2[0]=0.0; c3[0]=1.0;
+ return;
+ }
+
+ if(k==0) {
+ switch(int_j) {
+ case 0:
+ cos_taylor_qd(&c0[0],&c1[0],&c2[0],&c3[0],t0,t1,t2,t3);
+ return;
+ case 1:
+ sin_taylor_qd(&c0[0],&c1[0],&c2[0],&c3[0],t0,t1,t2,t3);
+ c0[0]=-c0[0]; c1[0]=-c1[0]; c2[0]=-c2[0]; c3[0]=-c3[0];
+ return;
+ case -1:
+ sin_taylor_qd(&c0[0],&c1[0],&c2[0],&c3[0],t0,t1,t2,t3);
+ return;
+ case 2:
+ case -2:
+ cos_taylor_qd(&c0[0],&c1[0],&c2[2],&c3[0],t0,t1,t2,t3);
+ c0[0]=-c0[0]; c1[0]=-c1[0]; c2[0]=-c2[0]; c3[0]=-c3[0];
+ return;
+ }
+ }
+
+ cos_table_qd(&u0,&u1,&u2,&u3,abs_k);
+ sin_table_qd(&v0,&v1,&v2,&v3,abs_k);
+ sincos_taylor_qd(&sin0,&sin1,&sin2,&sin3,&cos0,&cos1,&cos2,&cos3,t0,t1,t2,t3);
+
+ if(j==0) {
+ if(k>0) {
+ //u * cos_t - v * sin_t;
+ qd_mul_qd(&p0,&p1,&p2,&p3,u0,u1,u2,u3,cos0,cos1,cos2,cos3);
+ qd_mul_qd(&q0,&q1,&q2,&q3,v0,v1,v2,v3,sin0,sin1,sin2,sin3);
+ qd_sub_qd(&c0[0],&c1[0],&c2[0],&c3[0],p0,p1,p2,p3,q0,q1,q2,q3);
+ }
+ else {
+ //u * cos_t + v * sin_t;
+ qd_mul_qd(&p0,&p1,&p2,&p3,u0,u1,u2,u3,cos0,cos1,cos2,cos3);
+ qd_mul_qd(&q0,&q1,&q2,&q3,v0,v1,v2,v3,sin0,sin1,sin2,sin3);
+ qd_add_qd(&c0[0],&c1[0],&c2[0],&c3[0],p0,p1,p2,p3,q0,q1,q2,q3);
+ }
+ }
+ else if(j==1) {
+ if(k>0) {
+ //-u * sin_t - v * cos_t;
+ qd_mul_qd(&p0,&p1,&p2,&p3,u0,u1,u2,u3,sin0,sin1,sin2,sin3);
+ qd_mul_qd(&q0,&q1,&q2,&q3,v0,v1,v2,v3,cos0,cos1,cos2,cos3);
+ p0=-p0; p1=-p1; p2=-p2; p3=-p3;
+ qd_sub_qd(&c0[0],&c1[0],&c2[0],&c3[0],p0,p1,p2,p3,q0,q1,q2,q3);
+ }
+ else {
+ //v * cos_t - u * sin_t;
+ qd_mul_qd(&p0,&p1,&p2,&p3,v0,v1,v2,v3,cos0,cos1,cos2,cos3);
+ qd_mul_qd(&q0,&q1,&q2,&q3,u0,u1,u2,u3,sin0,sin1,sin2,sin3);
+ qd_sub_qd(&c0[0],&c1[0],&c2[0],&c3[0],p0,p1,p2,p3,q0,q1,q2,q3);
+ }
+ }
+ else if(j==-1) {
+ if(k>0) {
+ //u * sin_t + v * cos_t;
+ qd_mul_qd(&p0,&p1,&p2,&p3,u0,u1,u2,u3,sin0,sin1,sin2,sin3);
+ qd_mul_qd(&q0,&q1,&q2,&q3,v0,v1,v2,v3,cos0,cos1,cos2,cos3);
+ qd_add_qd(&c0[0],&c1[0],&c2[0],&c3[0],p0,p1,p2,p3,q0,q1,q2,q3);
+ }
+ else {
+ //u * sin_t - v * cos_t;
+ qd_mul_qd(&p0,&p1,&p2,&p3,u0,u1,u2,u3,sin0,sin1,sin2,sin3);
+ qd_mul_qd(&q0,&q1,&q2,&q3,v0,v1,v2,v3,cos0,cos1,cos2,cos3);
+ qd_sub_qd(&c0[0],&c1[0],&c2[0],&c3[0],p0,p1,p2,p3,q0,q1,q2,q3);
+ }
+ }
+ else {
+ if(k>0) {
+ //v * sin_t - u * cos_t;
+ qd_mul_qd(&p0,&p1,&p2,&p3,v0,v1,v2,v3,sin0,sin1,sin2,sin3);
+ qd_mul_qd(&q0,&q1,&q2,&q3,u0,u1,u2,u3,cos0,cos1,cos2,cos3);
+ qd_sub_qd(&c0[0],&c1[0],&c2[0],&c3[0],p0,p1,p2,p3,q0,q1,q2,q3);
+ }
+ else {
+ //-u * cos_t - v * sin_t;
+ qd_mul_qd(&p0,&p1,&p2,&p3,u0,u1,u2,u3,cos0,cos1,cos2,cos3);
+ qd_mul_qd(&q0,&q1,&q2,&q3,v0,v1,v2,v3,sin0,sin1,sin2,sin3);
+ p0=-p0; p1=-p1; p2=-p2; p3=-p3;
+ qd_sub_qd(&c0[0],&c1[0],&c2[0],&c3[0],p0,p1,p2,p3,q0,q1,q2,q3);
+ }
+ }
+ return;
+}
+
+//--------------------------------------------------------------------------------------------
+//
+// quad-double sine and cosine
+//
+//--------------------------------------------------------------------------------------------
+static INLINE void
+sincos_qd(double *s0, double *s1, double *s2, double *s3, 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
+ double _pi2[4] = {1.570796326794896558e+00,6.123233995736766036e-17,-1.497384904859169833e-33,5.562271104316826408e-50}; // pi/2
+ double _pi1024[4] = {3.067961575771282340e-03,1.195944139792337116e-19,-2.924579892303066080e-36,1.086381075061880158e-52};
+ double u0,u1,u2,u3,v0,v1,v2,v3;
+ double sin0,sin1,sin2,sin3,cos0,cos1,cos2,cos3;
+ int int_j;
+
+ if(a0==0.0) {
+ s0[0]=0.0; s1[0]=0.0; s2[0]=0.0; s3[0]=0.0;
+ c0[0]=1.0; c1[0]=0.0; c2[0]=0.0; c3[0]=0.0;
+ return;
+ }
+
+ //approximately reduce modulo 2*pi
+ qd_div_qd(&p0,&p1,&p2,&p3,a0,a1,a2,a3,_2pi[0],_2pi[1],_2pi[2],_2pi[3]);
+ nint_qd(&z0,&z1,&z2,&z3,p0,p1,p2,p3);
+ qd_mul_qd(&q0,&q1,&q2,&q3,_2pi[0],_2pi[1],_2pi[2],_2pi[3],z0,z1,z2,z3);
+ qd_sub_qd(&r0,&r1,&r2,&r3,a0,a1,a2,a3,q0,q1,q2,q3);
+
+ //approximately reduce modulo pi/2 and then modulo pi/1024
+ j=floor(r0/_pi2[0]+0.5);
+ s_mul_qd(&p0, &p1, &p2, &p3, j, _pi2[0], _pi2[1], _pi2[2], _pi2[3]); // Modified,2012/01/16 Saito
+ qd_sub_qd(&t0,&t1,&t2,&t3,r0,r1,r2,r3,p0,p1,p2,p3);
+ k=floor(t0/_pi1024[0]+0.5);
+ s_mul_qd(&q0, &q1, &q2, &q3, k, _pi1024[0], _pi1024[1], _pi1024[2], _pi1024[3]); // Modified,2012/01/16 Saito
+ qd_sub_qd(&t0,&t1,&t2,&t3,t0,t1,t2,t3,q0,q1,q2,q3);
+ abs_k=(int)fabs(k);
+ int_j=(int)j;
+
+ //checking errors
+ if(j<-2 || j>2) {
+ s0[0]=0.0; s1[0]=0.0; s2[0]=0.0; s3[0]=1.0;
+ c0[0]=0.0; c1[0]=0.0; c2[0]=0.0; c3[0]=1.0;
+ return;
+ }
+
+ if(abs_k >256) {
+ s0[0]=0.0; s1[0]=0.0; s2[0]=0.0; s3[0]=1.0;
+ c0[0]=0.0; c1[0]=0.0; c2[0]=0.0; c3[0]=1.0;
+ return;
+ }
+
+ sincos_taylor_qd(&sin0,&sin1,&sin2,&sin3,&cos0,&cos1,&cos2,&cos3,t0,t1,t2,t3);
+
+ if(k==0) {
+ if(j==0) {
+ s0[0]=sin0; s1[0]=sin1; s2[0]=sin2; s3[0]=sin3;
+ c0[0]=cos0; c1[0]=cos1; c2[0]=cos2; c3[0]=cos3;
+ return;
+ }
+ else if(j==1) {
+ s0[0]=cos0; s1[0]=cos1; s2[0]=cos2; s3[0]=cos3;
+ c0[0]=-sin0; c1[0]=-sin1; c2[0]=-sin2; c3[0]=-sin3;
+ return;
+ }
+ else if(j==-1) {
+ s0[0]=-cos0; s1[0]=-cos1; s2[0]=-cos2; s3[0]=-cos3;
+ c0[0]=sin0; c1[0]=sin1; c2[0]=sin2; c3[0]=sin3;
+ return;
+ }
+ else {
+ s0[0]=-sin0; s1[0]=-sin1; s2[0]=-sin2; s3[0]=-sin3;
+ c0[0]=-cos0; c1[0]=-cos1; c2[0]=-cos2; c3[0]=-cos3;
+ return;
+ }
+ }
+
+ cos_table_qd(&u0,&u1,&u2,&u3,abs_k);
+ sin_table_qd(&v0,&v1,&v2,&v3,abs_k);
+
+ if(j==0) {
+ if(k>0) {
+ qd_mul_qd(&p0,&p1,&p2,&p3,u0,u1,u2,u3,sin0,sin1,sin2,sin3);
+ qd_mul_qd(&q0,&q1,&q2,&q3,v0,v1,v2,v3,cos0,cos1,cos2,cos3);
+ qd_add_qd(&s0[0],&s1[0],&s2[0],&s3[0],p0,p1,p2,p3,q0,q1,q2,q3);
+
+ qd_mul_qd(&p0,&p1,&p2,&p3,u0,u1,u2,u3,cos0,cos1,cos2,cos3);
+ qd_mul_qd(&q0,&q1,&q2,&q3,v0,v1,v2,v3,sin0,sin1,sin2,sin3);
+ qd_sub_qd(&c0[0],&c1[0],&c2[0],&c3[0],p0,p1,p2,p3,q0,q1,q2,q3);
+ }
+ else {
+ qd_mul_qd(&p0,&p1,&p2,&p3,u0,u1,u2,u3,sin0,sin1,sin2,sin3);
+ qd_mul_qd(&q0,&q1,&q2,&q3,v0,v1,v2,v3,cos0,cos1,cos2,cos3);
+ qd_sub_qd(&s0[0],&s1[0],&s2[0],&s3[0],p0,p1,p2,p3,q0,q1,q2,q3);
+
+ qd_mul_qd(&p0,&p1,&p2,&p3,u0,u1,u2,u3,cos0,cos1,cos2,cos3);
+ qd_mul_qd(&q0,&q1,&q2,&q3,v0,v1,v2,v3,sin0,sin1,sin2,sin3);
+ qd_add_qd(&c0[0],&c1[0],&c2[0],&c3[0],p0,p1,p2,p3,q0,q1,q2,q3);
+ }
+ }
+ else if(j==1) {
+ if(k>0) {
+ qd_mul_qd(&p0,&p1,&p2,&p3,u0,u1,u2,u3,cos0,cos1,cos2,cos3);
+ qd_mul_qd(&q0,&q1,&q2,&q3,v0,v1,v2,v3,sin0,sin1,sin2,sin3);
+ qd_sub_qd(&s0[0],&s1[0],&s2[0],&s3[0],p0,p1,p2,p3,q0,q1,q2,q3);
+
+ qd_mul_qd(&p0,&p1,&p2,&p3,u0,u1,u2,u3,sin0,sin1,sin2,sin3);
+ qd_mul_qd(&q0,&q1,&q2,&q3,v0,v1,v2,v3,cos0,cos1,cos2,cos3);
+ p0=-p0; p1=-p1; p2=-p2; p3=-p3;
+ qd_sub_qd(&c0[0],&c1[0],&c2[0],&c3[0],p0,p1,p2,p3,q0,q1,q2,q3);
+ }
+ else {
+ qd_mul_qd(&p0,&p1,&p2,&p3,u0,u1,u2,u3,cos0,cos1,cos2,cos3);
+ qd_mul_qd(&q0,&q1,&q2,&q3,v0,v1,v2,v3,sin0,sin1,sin2,sin3);
+ qd_add_qd(&s0[0],&s1[0],&s2[0],&s3[0],p0,p1,p2,p3,q0,q1,q2,q3);
+
+ qd_mul_qd(&p0,&p1,&p2,&p3,v0,v1,v2,v3,cos0,cos1,cos2,cos3);
+ qd_mul_qd(&q0,&q1,&q2,&q3,u0,u1,u2,u3,sin0,sin1,sin2,sin3);
+ qd_sub_qd(&c0[0],&c1[0],&c2[0],&c3[0],p0,p1,p2,p3,q0,q1,q2,q3);
+ }
+ }
+ else if(j==-1) {
+ if(k>0) {
+ qd_mul_qd(&p0,&p1,&p2,&p3,v0,v1,v2,v3,sin0,sin1,sin2,sin3);
+ qd_mul_qd(&q0,&q1,&q2,&q3,u0,u1,u2,u3,cos0,cos1,cos2,cos3);
+ qd_sub_qd(&s0[0],&s1[0],&s2[0],&s3[0],p0,p1,p2,p3,q0,q1,q2,q3);
+
+ qd_mul_qd(&p0,&p1,&p2,&p3,u0,u1,u2,u3,sin0,sin1,sin2,sin3);
+ qd_mul_qd(&q0,&q1,&q2,&q3,v0,v1,v2,v3,cos0,cos1,cos2,cos3);
+ qd_add_qd(&c0[0],&c1[0],&c2[0],&c3[0],p0,p1,p2,p3,q0,q1,q2,q3);
+ }
+ else {
+ qd_mul_qd(&p0,&p1,&p2,&p3,u0,u1,u2,u3,cos0,cos1,cos2,cos3);
+ qd_mul_qd(&q0,&q1,&q2,&q3,v0,v1,v2,v3,sin0,sin1,sin2,sin3);
+ p0=-p0; p1=-p1; p2=-p2; p3=-p3;
+ qd_sub_qd(&s0[0],&s1[0],&s2[0],&s3[0],p0,p1,p2,p3,q0,q1,q2,q3);
+
+ qd_mul_qd(&p0,&p1,&p2,&p3,u0,u1,u2,u3,sin0,sin1,sin2,sin3);
+ qd_mul_qd(&q0,&q1,&q2,&q3,v0,v1,v2,v3,cos0,cos1,cos2,cos3);
+ qd_sub_qd(&c0[0],&c1[0],&c2[0],&c3[0],p0,p1,p2,p3,q0,q1,q2,q3);
+ }
+ }
+
+ else {
+ if(k>0) {
+ qd_mul_qd(&p0,&p1,&p2,&p3,u0,u1,u2,u3,sin0,sin1,sin2,sin3);
+ qd_mul_qd(&q0,&q1,&q2,&q3,v0,v1,v2,v3,cos0,cos1,cos2,cos3);
+ p0=-p0; p1=-p1; p2=-p2; p3=-p3;
+ qd_sub_qd(&s0[0],&s1[0],&s2[0],&s3[0],p0,p1,p2,p3,q0,q1,q2,q3);
+
+ qd_mul_qd(&p0,&p1,&p2,&p3,v0,v1,v2,v3,sin0,sin1,sin2,sin3);
+ qd_mul_qd(&q0,&q1,&q2,&q3,u0,u1,u2,u3,cos0,cos1,cos2,cos3);
+ qd_sub_qd(&c0[0],&c1[0],&c2[0],&c3[0],p0,p1,p2,p3,q0,q1,q2,q3);
+ }
+ else {
+ qd_mul_qd(&p0,&p1,&p2,&p3,v0,v1,v2,v3,cos0,cos1,cos2,cos3);
+ qd_mul_qd(&q0,&q1,&q2,&q3,u0,u1,u2,u3,sin0,sin1,sin2,sin3);
+ qd_sub_qd(&s0[0],&s1[0],&s2[0],&s3[0],p0,p1,p2,p3,q0,q1,q2,q3);
+
+ qd_mul_qd(&p0,&p1,&p2,&p3,u0,u1,u2,u3,cos0,cos1,cos2,cos3);
+ qd_mul_qd(&q0,&q1,&q2,&q3,v0,v1,v2,v3,sin0,sin1,sin2,sin3);
+ p0=-p0; p1=-p1; p2=-p2; p3=-p3;
+ qd_sub_qd(&c0[0],&c1[0],&c2[0],&c3[0],p0,p1,p2,p3,q0,q1,q2,q3);
+ }
+ }
+}
+
+//--------------------------------------------------------------------------------------------
+//
+// quad-double tangent
+//
+// args
+// a0, a1, a2, a3 : double numbers
+// a0 + a1 + a2 + a3 = qd number
+//
+// return (t0,t1,t2,t3) for qd number
+//
+//--------------------------------------------------------------------------------------------
+static void
+qdtan(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]);
+ qd_div_qd(&t0[0],&t1[0],&t2[0],&t3[0],sin0,sin1,sin2,sin3,cos0,cos1,cos2,cos3);
+}
+
+//--------------------------------------------------------------------------------------------
+//
+// quad-double exponent
+//
+// args
+// x0, x1, x2, x3 : double numbers
+// x0 + x1 + x2 + x3 = qd number
+//
+// return (e0, e1, e2, e3) for qd number
+//
+//--------------------------------------------------------------------------------------------
+static void
+qdexp(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;
+ double log_2[4] = {6.931471805599452862e-01,2.319046813846299558e-17,5.707708438416212066e-34,-3.582432210601811423e-50};
+ double m, p0,p1,p2,p3, q0,q1,q2,q3, r0,r1,r2,r3, s0,s1,s2,s3, t0,t1,t2,t3, v0,v1,v2,v3;
+ int i;
+ double t=1.0;
+ double eps = 1.21543267145725e-63; // = 2^-209
+ double thresh = inv_k * eps;
+
+ if((x0[0]==0.0) && (x1[0]==0.0) && (x2[0]==0.0) && (x3[0]==0.0)) {
+ e0[0] = 1.0; e1[0] = 0.0; e2[0] = 0.0; e3[0] = 0.0;
+ return;
+ }
+
+ if((x0[0]==1.0) && (x1[0]==0.0) && (x2[0]==0.0) && (x3[0]==0.0)) {
+ e0[0] = 2.7182818284590451;
+ e1[0] = 1.4456468917292502e-16;
+ e2[0] = -2.127717108038176765e-33;
+ e3[0] = 1.515630159841218954e-49;
+ return;
+ }
+
+ if(x0[0]<=-709) { //underflow
+ e0[0] = 0.0; e1[0] = 0.0; e2[0] = 0.0; e3[0] = 0.0;
+ return;
+ }
+
+ if(x0[0]>=709) { //overflow
+ e0[0] = 0.0; e1[0] = 1.0; e2[0] = 2.0; e3[0] = 3.0;
+ return;
+ }
+
+ m = floor(x0[0] / log_2[0] + 0.5);
+
+ s_mul_qd(&p0, &p1, &p2, &p3, m, log_2[0], log_2[1], log_2[2], log_2[3]);
+ qd_sub_qd(&q0, &q1, &q2, &q3, x0[0], x1[0], x2[0], x3[0], p0, p1, p2, p3);
+ r0 = q0 * inv_k;
+ r1 = q1 * inv_k;
+ r2 = q2 * inv_k;
+ r3 = q3 * inv_k;
+
+ qd_sqr(&p0, &p1, &p2, &p3, r0, r1, r2, r3);
+ qd_add_qd(&s0, &s1, &s2, &s3, r0, r1, r2, r3, 0.5*p0, 0.5*p1, 0.5*p2, 0.5*p3);
+ i = 0;
+ do {
+ qd_mul_qd(&p0, &p1, &p2, &p3, p0, p1, p2, p3, r0, r1, r2, r3);
+ qd_mul_qd(&t0, &t1, &t2, &t3, p0, p1, p2, p3, inv_fact[i][0], inv_fact[i][1], inv_fact[i][2], inv_fact[i][3]);
+ i = i+1;
+ qd_add_qd(&s0, &s1, &s2, &s3, s0, s1, s2, s3, t0, t1, t2, t3);
+ } while ((i<=17) && (fabs(t0)>thresh));
+
+
+ qd_sqr(&v0, &v1, &v2, &v3, s0, s1, s2, s3);
+ qd_add_qd(&s0, &s1, &s2, &s3, 2.0*s0, 2.0*s1, 2.0*s2, 2.0*s3, v0, v1, v2, v3); //1
+ qd_sqr(&v0, &v1, &v2, &v3, s0, s1, s2, s3);
+ qd_add_qd(&s0, &s1, &s2, &s3, 2.0*s0, 2.0*s1, 2.0*s2, 2.0*s3, v0, v1, v2, v3); //2
+ qd_sqr(&v0, &v1, &v2, &v3, s0, s1, s2, s3);
+ qd_add_qd(&s0, &s1, &s2, &s3, 2.0*s0, 2.0*s1, 2.0*s2, 2.0*s3, v0, v1, v2, v3); //3
+ qd_sqr(&v0, &v1, &v2, &v3, s0, s1, s2, s3);
+ qd_add_qd(&s0, &s1, &s2, &s3, 2.0*s0, 2.0*s1, 2.0*s2, 2.0*s3, v0, v1, v2, v3); //4
+ qd_sqr(&v0, &v1, &v2, &v3, s0, s1, s2, s3);
+ qd_add_qd(&s0, &s1, &s2, &s3, 2.0*s0, 2.0*s1, 2.0*s2, 2.0*s3, v0, v1, v2, v3); //5
+ qd_sqr(&v0, &v1, &v2, &v3, s0, s1, s2, s3);
+ qd_add_qd(&s0, &s1, &s2, &s3, 2.0*s0, 2.0*s1, 2.0*s2, 2.0*s3, v0, v1, v2, v3); //6
+ qd_sqr(&v0, &v1, &v2, &v3, s0, s1, s2, s3);
+ qd_add_qd(&s0, &s1, &s2, &s3, 2.0*s0, 2.0*s1, 2.0*s2, 2.0*s3, v0, v1, v2, v3); //7
+ qd_sqr(&v0, &v1, &v2, &v3, s0, s1, s2, s3);
+ qd_add_qd(&s0, &s1, &s2, &s3, 2.0*s0, 2.0*s1, 2.0*s2, 2.0*s3, v0, v1, v2, v3); //8
+ qd_sqr(&v0, &v1, &v2, &v3, s0, s1, s2, s3);
+ qd_add_qd(&s0, &s1, &s2, &s3, 2.0*s0, 2.0*s1, 2.0*s2, 2.0*s3, v0, v1, v2, v3); //9
+ qd_sqr(&v0, &v1, &v2, &v3, s0, s1, s2, s3);
+ qd_add_qd(&s0, &s1, &s2, &s3, 2.0*s0, 2.0*s1, 2.0*s2, 2.0*s3, v0, v1, v2, v3); //10
+ qd_sqr(&v0, &v1, &v2, &v3, s0, s1, s2, s3);
+ qd_add_qd(&s0, &s1, &s2, &s3, 2.0*s0, 2.0*s1, 2.0*s2, 2.0*s3, v0, v1, v2, v3); //11
+ qd_sqr(&v0, &v1, &v2, &v3, s0, s1, s2, s3);
+ qd_add_qd(&s0, &s1, &s2, &s3, 2.0*s0, 2.0*s1, 2.0*s2, 2.0*s3, v0, v1, v2, v3); //12
+ qd_sqr(&v0, &v1, &v2, &v3, s0, s1, s2, s3);
+ qd_add_qd(&s0, &s1, &s2, &s3, 2.0*s0, 2.0*s1, 2.0*s2, 2.0*s3, v0, v1, v2, v3); //13
+ qd_sqr(&v0, &v1, &v2, &v3, s0, s1, s2, s3);
+ qd_add_qd(&s0, &s1, &s2, &s3, 2.0*s0, 2.0*s1, 2.0*s2, 2.0*s3, v0, v1, v2, v3); //14
+ qd_sqr(&v0, &v1, &v2, &v3, s0, s1, s2, s3);
+ qd_add_qd(&s0, &s1, &s2, &s3, 2.0*s0, 2.0*s1, 2.0*s2, 2.0*s3, v0, v1, v2, v3); //15
+ qd_sqr(&v0, &v1, &v2, &v3, s0, s1, s2, s3);
+ qd_add_qd(&s0, &s1, &s2, &s3, 2.0*s0, 2.0*s1, 2.0*s2, 2.0*s3, v0, v1, v2, v3); //16
+ qd_add_s(&s0, &s1, &s2, &s3, s0, s1, s2, s3, 1.0);
+
+ for(i=0; i<m; i++) {
+ t=t*2;
+ }
+
+ e0[0] = s0*t; e1[0] = s1*t; e2[0] = s2*t; e3[0] = s3*t;
+ return;
+}
+
#if 0
/*********** Basic Functions ************/
/* Computes fl(a+b) and err(a+b). Assumes |a| >= |b|. */
-inline double
+INLINE double
quick_two_sum(double a, double b, double *errPtr) {
double s = a + b;
*errPtr = b - (s - a);
@@ -390,7 +2372,7 @@
}
/* Computes fl(a-b) and err(a-b). Assumes |a| >= |b| */
-inline double
+INLINE double
quick_two_diff(double a, double b, double *errPtr) {
double s = a - b;
*errPtr = (a - s) - b;
@@ -398,7 +2380,7 @@
}
/* Computes fl(a+b) and err(a+b). */
-inline double
+INLINE double
two_sum(double a, double b, double *errPtr) {
double s = a + b;
double bb = s - a;
@@ -407,7 +2389,7 @@
}
/* Computes fl(a-b) and err(a-b). */
-inline double
+INLINE double
two_diff(double a, double b, double *errPtr) {
double s = a - b;
double bb = s - a;
@@ -417,7 +2399,7 @@
#ifndef QD_FMS
/* Computes high word and lo word of a */
-inline void
+INLINE void
split(double a, double *hiPtr, double *loPtr) {
double temp;
if (a > _QD_SPLIT_THRESH || a < -_QD_SPLIT_THRESH) {
@@ -436,7 +2418,7 @@
#endif
/* Computes fl(a*b) and err(a*b). */
-inline double
+INLINE double
two_prod(double a, double b, double *errPtr) {
#ifdef QD_FMS
double p = a * b;
@@ -453,7 +2435,7 @@
}
/* Computes fl(a*a) and err(a*a). Faster than the above method. */
-inline double
+INLINE double
two_sqr(double a, double *errPtr) {
#ifdef QD_FMS
double p = a * a;
@@ -469,7 +2451,7 @@
}
/* Computes the nearest integer to d. */
-inline double
+INLINE double
nint(double d) {
if (d == floor(d))
return d;
@@ -477,12 +2459,12 @@
}
/* Computes the truncated integer. */
-inline double
+INLINE double
aint(double d) {
return (d >= 0.0) ? floor(d) : ceil(d);
}
-inline void
+INLINE void
renorm4(double *c0Ptr, double *c1Ptr, double *c2Ptr, double *c3Ptr) {
double s0, s1, s2 = 0.0, s3 = 0.0;
double c0 = *c0Ptr;
@@ -515,7 +2497,7 @@
*c3Ptr = s3;
}
-inline void
+INLINE void
renorm5(double *c0Ptr, double *c1Ptr, double *c2Ptr, double *c3Ptr, double *c4Ptr) {
double s0, s1, s2 = 0.0, s3 = 0.0;
@@ -568,7 +2550,7 @@
*c3Ptr = s3;
}
-inline void
+INLINE void
three_sum(double *aPtr, double *bPtr, double *cPtr) {
double t1, t2, t3;
t1 = two_sum(*aPtr, *bPtr, &t2);
@@ -576,7 +2558,7 @@
*bPtr = two_sum(t2, t3, cPtr);
}
-inline void three_sum2(double *aPtr, double *bPtr, double *cPtr) {
+INLINE void three_sum2(double *aPtr, double *bPtr, double *cPtr) {
double t1, t2, t3;
t1 = two_sum(*aPtr, *bPtr, &t2);
*aPtr = two_sum(*cPtr, t1, &t3);
@@ -587,23 +2569,23 @@
#if 0
/* These are provided to give consistent
interface for double with double-double and quad-double. */
-inline void
+INLINE void
sincosh(double t, double &sinh_t, double &cosh_t) {
sinh_t = sinh(t);
cosh_t = cosh(t);
}
-inline double
+INLINE double
sqr(double t) {
return t * t;
}
-inline double
+INLINE double
to_double(double a) {
return a;
}
-inline int
+INLINE int
to_int(double a) {
return static_cast<int>(a);
}
@@ -639,47 +2621,47 @@
QDoubles give you roughly 200 bit or approx. 60 decimal digits of precision.
Representation:
- QDoubles use 4 IEEE doubles, each keeping 53 bits of precision.
- A qDouble's value is the sum of those 4 doubles,
- and a qDouble keeps this unevaluated sum as its state.
- (due to overlap and rounding, the final precision is less than 53*4)
- The exponent range is still the double exponent range,
- but the number of mantissa bits is rougly multiplied by 4.
+ QDoubles use 4 IEEE doubles, each keeping 53 bits of precision.
+ A qDouble's value is the sum of those 4 doubles,
+ and a qDouble keeps this unevaluated sum as its state.
+ (due to overlap and rounding, the final precision is less than 53*4)
+ The exponent range is still the double exponent range,
+ but the number of mantissa bits is rougly multiplied by 4.
Range and Precision of Storage Formats: see LimitedPrecisionReal >> documentation
The number of decmal digits:
- QDouble decimalPrecision -> 61
- LongFloat decimalPrecision -> 19
- Float decimalPrecision -> 16
- ShortFloat decimalPrecision -> 7
+ QDouble decimalPrecision -> 61
+ LongFloat decimalPrecision -> 19
+ Float decimalPrecision -> 16
+ ShortFloat decimalPrecision -> 7
The number of bits:
- QDouble precision -> 204
- LongFloat precision -> 64
- Float precision -> 53
- ShortFloat precision -> 24
+ QDouble precision -> 204
+ LongFloat precision -> 64
+ Float precision -> 53
+ ShortFloat precision -> 24
Notice:
- when assigning a converted double precision number as in:
- qd := 1.0 asQDouble.
- you still get only a regular double precision approximation to 0.1
- because the error is already inherit in the double.
-
- For a full precision constant, you (currently) need to convert from a string
- (because the compilers do not know about them, yet):
- qd := QDouble readFrom:'0.1'.
-
- To see the error of the double precision version, compute:
- (0.1 asQDouble) - (QDouble readFrom:'0.1')
+ when assigning a converted double precision number as in:
+ qd := 1.0 asQDouble.
+ you still get only a regular double precision approximation to 0.1
+ because the error is already inherit in the double.
+
+ For a full precision constant, you (currently) need to convert from a string
+ (because the compilers do not know about them, yet):
+ qd := QDouble readFrom:'0.1'.
+
+ To see the error of the double precision version, compute:
+ (0.1 asQDouble) - (QDouble readFrom:'0.1')
[author:]
- Claus Gittinger
+ Claus Gittinger
[see also:]
- Number
- Float ShortFloat LongFloat
- Fraction FixedPoint Integer Complex
- FloatArray DoubleArray
+ Number
+ Float ShortFloat LongFloat
+ Fraction FixedPoint Integer Complex
+ FloatArray DoubleArray
"
!
@@ -735,17 +2717,12 @@
%{ /* NOCONTEXT */
#ifdef __SCHTEAM__
- ERROR("trying to instantiate a quad double");
+ ERROR("trying to instantiate a qDouble");
#else
- OBJ newFloat;
-
- __qNew(newFloat, sizeof(struct __quadDoubleStruct));
- __stx_setClass(newFloat, QDouble);
- __QuadDoubleInstPtr(newFloat)->d_quadDoubleValue[0] = 0.0;
- __QuadDoubleInstPtr(newFloat)->d_quadDoubleValue[1] = 0.0;
- __QuadDoubleInstPtr(newFloat)->d_quadDoubleValue[2] = 0.0;
- __QuadDoubleInstPtr(newFloat)->d_quadDoubleValue[3] = 0.0;
- RETURN (newFloat);
+ OBJ newQD;
+
+ __qNew_qdReal(newQD, 0.0, 0.0, 0.0, 0.0);
+ RETURN (newQD);
#endif /* not SCHTEAM */
%}
@@ -761,7 +2738,7 @@
%{ /* NOCONTEXT */
#ifdef __SCHTEAM__
- ERROR("trying to instantiate a quad double");
+ ERROR("trying to instantiate a qDouble");
#else
OBJ newQD;
@@ -792,16 +2769,13 @@
%{ /* NOCONTEXT */
#ifdef __SCHTEAM__
- ERROR("trying to instantiate a quad double");
+ ERROR("trying to instantiate a qDouble");
#else
OBJ newQD;
if (__isDoubleArray(aDoubleArray)) {
- __qNew_qdReal(newQD,
- __DoubleArrayInstPtr(aDoubleArray)->d_element[0],
- __DoubleArrayInstPtr(aDoubleArray)->d_element[1],
- __DoubleArrayInstPtr(aDoubleArray)->d_element[2],
- __DoubleArrayInstPtr(aDoubleArray)->d_element[3]);
+ double* __d__ = __DoubleArrayInstPtr(aDoubleArray)->d_element;
+ __qNew_qdReal(newQD, __d__[0], __d__[1], __d__[2], __d__[3]);
RETURN (newQD);
}
#endif
@@ -824,10 +2798,10 @@
%{ /* NOCONTEXT */
#ifdef __SCHTEAM__
- ERROR("trying to instantiate a quad double");
+ ERROR("trying to instantiate a qDouble");
#else
double dVal;
- OBJ newFloat;
+ OBJ newQD;
if (__isFloatLike(aFloat)) {
dVal = __floatVal(aFloat);
@@ -837,19 +2811,14 @@
goto badArg;
}
- __qNew(newFloat, sizeof(struct __quadDoubleStruct));
- __stx_setClass(newFloat, QDouble);
- __QuadDoubleInstPtr(newFloat)->d_quadDoubleValue[0] = dVal;
- __QuadDoubleInstPtr(newFloat)->d_quadDoubleValue[1] = 0.0;
- __QuadDoubleInstPtr(newFloat)->d_quadDoubleValue[2] = 0.0;
- __QuadDoubleInstPtr(newFloat)->d_quadDoubleValue[3] = 0.0;
- RETURN (newFloat);
+ __qNew_qdReal(newQD, dVal, 0.0, 0.0, 0.0);
+ RETURN (newQD);
badArg: ;
#endif
%}.
- self error:'invalid argument'
+ self argumentError:'invalid (non-float) argument'
"
self fromFloat:1.0
@@ -862,32 +2831,34 @@
"return a new quad-precision double from coercing anInteger"
%{ /* NOCONTEXT */
-#ifndef __SCHTEAM__
- OBJ newFloat;
+#ifdef __SCHTEAM__
+ ERROR("trying to instantiate a qDouble");
+#else
+ OBJ newQD;
if (__isSmallInteger(anInteger)) {
INT iVal = __smallIntegerVal(anInteger);
- double *a;
-
- __qNew(newFloat, sizeof(struct __quadDoubleStruct));
- __stx_setClass(newFloat, QDouble);
-
- a = __QuadDoubleInstPtr(newFloat)->d_quadDoubleValue;
- a[1] = 0.0;
- a[2] = 0.0;
- a[3] = 0.0;
+ 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;
// need more than 52bits?
if ((sizeof(INT) > 52)
- && ((iVal > 0xFFFFFFFFFFFFF)
- || (iVal < -0xFFFFFFFFFFFFF))) {
- a[0] = (double)(iVal & ~0xFFFFFFFF);
- a[1] = (double)(iVal & 0xFFFFFFFF);
- // m_renorm4(a[0], a[1], a[2], a[3]);
+ && ((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);
+ // renorm4(&(a[0]), &(a[1]), &(a[2]), &(a[3]));
} else {
- a[0] = (double)iVal;
+ __d__[0] = (double)iVal;
}
- RETURN (newFloat);
+ RETURN (newQD);
}
#endif
%}.
@@ -939,18 +2910,18 @@
(returns approx. 200 bits of precision)"
E isNil ifTrue:[
- E := self d0: 2.718281828459045091e+00
- d1: 1.445646891729250158e-16
- d2: -2.127717108038176765e-33
- d3: 1.515630159841218954e-49
+ E := self d0: 2.718281828459045091e+00
+ d1: 1.445646891729250158e-16
+ d2: -2.127717108038176765e-33
+ d3: 1.515630159841218954e-49
].
^ E
"
- self e printfPrintString:'%.61f'
+ self e printfPrintString:'%.61f'
-> '2.7182818284590452353602874713526624977572470936999595749669676'
Wolfram says:
- 2.71828182845904523536028747135266249775724709369995957496696762772407663035354759457138217852516642742746
+ 2.71828182845904523536028747135266249775724709369995957496696762772407663035354759457138217852516642742746
"
"Created: / 12-06-2017 / 18:29:36 / cg"
@@ -961,10 +2932,10 @@
(returns approx. 200 bits of precision)"
FMax isNil ifTrue:[
- FMax := self d0: 1.797693134862314E+308
- d1: 9.97920154767359795037e+291
- d2: 5.53956966280111259858e+275
- d3: 3.07507889307840487279e+259
+ FMax := self d0: 1.797693134862314E+308
+ d1: 9.97920154767359795037e+291
+ d2: 5.53956966280111259858e+275
+ d3: 3.07507889307840487279e+259
].
^ FMax
@@ -980,7 +2951,7 @@
"return the smallest representable instance of this class"
FMin isNil ifTrue:[
- FMin := 1.6259745436952323e-260 asQDouble
+ FMin := Float fmin asQDouble. "/ 1.6259745436952323e-260 asQDouble
].
^ FMin
@@ -998,54 +2969,6 @@
"Created: / 18-06-2017 / 23:41:07 / cg"
!
-invFact
- "table returning 1/n!!
- (for taylor series)"
-
- InvFact isNil ifTrue:[
- InvFact := Array new:15.
- InvFact at:1 put:(self d0:1.66666666666666657e-01 d1:9.25185853854297066e-18 d2:5.13581318503262866e-34 d3:2.85094902409834186e-50).
- InvFact at:2 put:(self d0:4.16666666666666644e-02 d1:2.31296463463574266e-18 d2:1.28395329625815716e-34 d3:7.12737256024585466e-51).
- InvFact at:3 put:(self d0:8.33333333333333322e-03 d1:1.15648231731787138e-19 d2:1.60494162032269652e-36 d3:2.22730392507682967e-53).
- InvFact at:4 put:(self d0:1.38888888888888894e-03 d1:-5.30054395437357706e-20 d2:-1.73868675534958776e-36 d3:-1.63335621172300840e-52).
- InvFact at:5 put:(self d0:1.98412698412698413e-04 d1:1.72095582934207053e-22 d2:1.49269123913941271e-40 d3:1.29470326746002471e-58).
- InvFact at:6 put:(self d0:2.48015873015873016e-05 d1:2.15119478667758816e-23 d2:1.86586404892426588e-41 d3:1.61837908432503088e-59).
- InvFact at:7 put:(self d0:2.75573192239858925e-06 d1:-1.85839327404647208e-22 d2:8.49175460488199287e-39 d3:-5.72661640789429621e-55).
- InvFact at:8 put:(self d0:2.75573192239858883e-07 d1:2.37677146222502973e-23 d2:-3.26318890334088294e-40 d3:1.61435111860404415e-56).
- InvFact at:9 put:(self d0:2.50521083854417202e-08 d1:-1.44881407093591197e-24 d2:2.04267351467144546e-41 d3:-8.49632672007163175e-58).
- InvFact at:10 put:(self d0:2.08767569878681002e-09 d1:-1.20734505911325997e-25 d2:1.70222792889287100e-42 d3:1.41609532150396700e-58).
- InvFact at:11 put:(self d0:1.60590438368216133e-10 d1:1.25852945887520981e-26 d2:-5.31334602762985031e-43 d3:3.54021472597605528e-59).
- InvFact at:12 put:(self d0:1.14707455977297245e-11 d1:2.06555127528307454e-28 d2:6.88907923246664603e-45 d3:5.72920002655109095e-61).
- InvFact at:13 put:(self d0:7.64716373181981641e-13 d1:7.03872877733453001e-30 d2:-7.82753927716258345e-48 d3:1.92138649443790242e-64).
- InvFact at:14 put:(self d0:4.77947733238738525e-14 d1:4.39920548583408126e-31 d2:-4.89221204822661465e-49 d3:1.20086655902368901e-65).
- InvFact at:15 put:(self d0:2.81145725434552060e-15 d1:1.65088427308614326e-31 d2:-2.87777179307447918e-50 d3:4.27110689256293549e-67).
- ].
- ^ InvFact
-
- "
- 1.0 / (3 factorial) 0.166666666666667
-
- 1 asQDouble / (3 factorial) - (self invFact at:1)
- 1 asQDouble / (4 factorial) - (self invFact at:2)
- 1 asQDouble / (5 factorial) - (self invFact at:3)
- 1 asQDouble / (6 factorial) - (self invFact at:4)
- 1 asQDouble / (7 factorial) - (self invFact at:5)
- 1 asQDouble / (8 factorial) - (self invFact at:6)
- 1 asQDouble / (9 factorial) - (self invFact at:7)
- 1 asQDouble / (10 factorial) - (self invFact at:8)
- 1 asQDouble / (11 factorial) - (self invFact at:9)
- 1 asQDouble / (12 factorial) - (self invFact at:10)
- 1 asQDouble / (13 factorial) - (self invFact at:11)
- 1 asQDouble / (14 factorial) - (self invFact at:12)
- 1 asQDouble / (15 factorial) - (self invFact at:13)
- 1 asQDouble / (16 factorial) - (self invFact at:14)
- 1 asQDouble / (17 factorial) - (self invFact at:15)
- "
-
- "Created: / 19-06-2017 / 02:22:23 / cg"
- "Modified (comment): / 20-06-2017 / 13:04:54 / cg"
-!
-
ln10
"return the constant e as quad precision double.
(returns approx. 200 bits of precision)"
@@ -1073,18 +2996,18 @@
(returns approx. 200 bits of precision)"
Ln2 isNil ifTrue:[
- Ln2 := self d0: 6.931471805599452862e-01
- d1: 2.319046813846299558e-17
- d2: 5.707708438416212066e-34
- d3: -3.582432210601811423e-50
+ Ln2 := self d0: 6.931471805599452862e-01
+ d1: 2.319046813846299558e-17
+ d2: 5.707708438416212066e-34
+ d3: -3.582432210601811423e-50
].
^ Ln2
"
- self ln2 printfPrintString:'%.61f'
- -> '0.6931471805599452709398341558750792990469129794959648865081141'
+ self ln2 printfPrintString:'%.61f'
+ -> '0.6931471805599452709398341558750792990469129794959648865081141'
Wolfram says:
- 0.69314718055994530941723212145817656807550013436025525412068000949339362196969471560586332699641868754200148102057...
+ 0.69314718055994530941723212145817656807550013436025525412068000949339362196969471560586332699641868754200148102057...
"
"Created: / 12-06-2017 / 18:31:34 / cg"
@@ -1101,21 +3024,21 @@
(returns approx. 200 bits of precision)"
Pi isNil ifTrue:[
- Pi := self d0: 3.141592653589793116e+00
- d1: 1.224646799147353207e-16
- d2: -2.994769809718339666e-33
- d3: 1.112454220863365282e-49
- ].
+ Pi := self d0: 3.141592653589793116e+00
+ d1: 1.224646799147353207e-16
+ d2: -2.994769809718339666e-33
+ d3: 1.112454220863365282e-49
+ ].
^ Pi
"
- self pi printfPrintString:'%.60f'
- '3.141592653589793238462643383279502884197169399375105820974945'
+ self pi printfPrintString:'%.60f'
+ '3.141592653589793238462643383279502884197169399375105820974945'
Wolfram says:
- 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068
+ 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068
(QDouble readFrom:'3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253')
- printfPrintString:'%.60f'
+ printfPrintString:'%.60f'
"
@@ -1160,12 +3083,12 @@
^ 30
"
- ShortFloat defaultPrintPrecision
- Float defaultPrintPrecision
- LongFloat defaultPrintPrecision
- QDouble defaultPrintPrecision
- QuadFloat defaultPrintPrecision
- OctaFloat defaultPrintPrecision
+ ShortFloat defaultPrintPrecision
+ Float defaultPrintPrecision
+ LongFloat defaultPrintPrecision
+ QDouble defaultPrintPrecision
+ QuadFloat defaultPrintPrecision
+ OctaFloat defaultPrintPrecision
"
"Created: / 17-06-2017 / 02:58:51 / cg"
@@ -1179,12 +3102,12 @@
"/ ^ 1.2154326714572500565324311366323150942261000827598106963711353e-63
Epsilon isNil ifTrue:[
- Epsilon := self computeEpsilon.
+ Epsilon := self computeEpsilon.
].
^ Epsilon
"
- Float epsilon -> 2.22044604925031E-16
+ Float epsilon -> 2.22044604925031E-16
ShortFloat epsilon -> 1.19209289550781E-07
LongFloat epsilon -> 1.0842021724855E-19
QDouble epsilon -> 7.77876909732643E-62 / (1.215432671457250056532e-63 read comment in precision)
@@ -1236,7 +3159,7 @@
"answer the number of bits in the mantissa"
"/ subtract some due to overlap in the component numbers
- "/ actual precision seems to be more like:
+ "/ actual precision seems to be more like:
"/ ^ (Float precision) * 4 - 3 + 1.
"/ but I am a bit conservative here:
^ (Float precision - 2) * 4
@@ -1245,7 +3168,7 @@
ShortFloat precision -> 24
Float precision -> 53
LongFloat precision -> 64
- QDouble precision -> 204
+ QDouble precision -> 204
QuadFloat precision -> 113
OctaFloat precision -> 237
@@ -1278,22 +3201,84 @@
* aNumber
"return the product of the receiver and the argument, aNumber"
+%{
+ if (__isFloatLike(aNumber)) {
+ double *a = __QDoubleInstPtr(self)->d_qDoubleValue;
+ double b = __floatVal(aNumber);
+ double c0, c1, c2, c3;
+ OBJ newQD;
+ int savedCV;
+
+ fpu_fix_start(&savedCV);
+ s_mul_qd(&c0, &c1, &c2, &c3, b, a[0], a[1], a[2], a[3]);
+ fpu_fix_end(&savedCV);
+ __qNew_qdReal(newQD, c0, c1, c2, c3);
+ RETURN( newQD );
+ }
+ if (__isQDouble(aNumber)) {
+ double *a = __QDoubleInstPtr(self)->d_qDoubleValue;
+ double *b = __QDoubleInstPtr(aNumber)->d_qDoubleValue;
+ double c0, c1, c2, c3;
+ OBJ newQD;
+ int savedCV;
+
+ fpu_fix_start(&savedCV);
+ qd_mul_qd(&c0, &c1, &c2, &c3, a[0], a[1], a[2], a[3], b[0], b[1], b[2], b[3]);
+ fpu_fix_end(&savedCV);
+ __qNew_qdReal(newQD, c0, c1, c2, c3);
+ RETURN( newQD );
+ }
+%}.
^ aNumber productFromQDouble:self
"
- (((QDouble fromFloat:1e20) + (QDouble fromFloat:1.0)) * (QDouble fromFloat:2.0)) asDoubleArray
+ (QDouble fromFloat:1e20) * 2.0
+ (QDouble fromFloat:1e20) * 1e20
+ (QDouble fromFloat:1e20) * (QDouble fromFloat:1e20)
+ ((QDouble fromFloat:1e20) * (QDouble fromFloat:2.0)) asDoubleArray
+ ((QDouble fromFloat:1e-20) * (QDouble fromFloat:2.0)) asDoubleArray
+ ((QDouble fromFloat:2.0) * (QDouble fromFloat:2.0)) asDoubleArray
"
- "Created: / 13-06-2017 / 01:00:47 / cg"
- "Modified (comment): / 14-06-2017 / 12:08:50 / cg"
+ "Created: / 12-06-2017 / 23:41:39 / cg"
+ "Modified (comment): / 15-06-2017 / 00:34:41 / cg"
!
+ aNumber
"return the sum of the receiver and the argument, aNumber"
+%{
+ if (__isFloatLike(aNumber)) {
+ double *a = __QDoubleInstPtr(self)->d_qDoubleValue;
+ double b = __floatVal(aNumber);
+ double c0, c1, c2, c3;
+ OBJ newQD;
+ int savedCV;
+
+ fpu_fix_start(&savedCV);
+ qd_add_s(&c0, &c1, &c2, &c3, a[0], a[1], a[2], a[3], b);
+ fpu_fix_end(&savedCV);
+ __qNew_qdReal(newQD, c0, c1, c2, c3);
+ RETURN( newQD );
+ }
+ if (__isQDouble(aNumber)) {
+ double *a = __QDoubleInstPtr(self)->d_qDoubleValue;
+ double *b = __QDoubleInstPtr(aNumber)->d_qDoubleValue;
+ double c0, c1, c2, c3;
+ OBJ newQD;
+ int savedCV;
+
+ fpu_fix_start(&savedCV);
+ qd_add_qd(&c0, &c1, &c2, &c3, a[0], a[1], a[2], a[3], b[0], b[1], b[2], b[3]);
+ fpu_fix_end(&savedCV);
+ __qNew_qdReal(newQD, c0, c1, c2, c3);
+ RETURN( newQD );
+ }
+%}.
^ aNumber sumFromQDouble:self
"
+ ((QDouble fromFloat:1e20) + 1.0) asDoubleArray
((QDouble fromFloat:1e20) + (QDouble fromFloat:1.0)) asDoubleArray
"
@@ -1304,10 +3289,38 @@
- aNumber
"return the sum of the receiver and the argument, aNumber"
+%{
+ if (__isFloatLike(aNumber)) {
+ double *a = __QDoubleInstPtr(self)->d_qDoubleValue;
+ double b = __floatVal(aNumber);
+ double c0, c1, c2, c3;
+ OBJ newQD;
+ int savedCV;
+
+ fpu_fix_start(&savedCV);
+ qd_add_s(&c0, &c1, &c2, &c3, a[0], a[1], a[2], a[3], -b);
+ fpu_fix_end(&savedCV);
+ __qNew_qdReal(newQD, c0, c1, c2, c3);
+ RETURN( newQD );
+ }
+ if (__isQDouble(aNumber)) {
+ double *a = __QDoubleInstPtr(self)->d_qDoubleValue;
+ double *b = __QDoubleInstPtr(aNumber)->d_qDoubleValue;
+ double c0, c1, c2, c3;
+ OBJ newQD;
+ int savedCV;
+
+ fpu_fix_start(&savedCV);
+ qd_sub_qd(&c0, &c1, &c2, &c3, a[0], a[1], a[2], a[3], b[0], b[1], b[2], b[3]);
+ fpu_fix_end(&savedCV);
+ __qNew_qdReal(newQD, c0, c1, c2, c3);
+ RETURN( newQD );
+ }
+%}.
^ self + (aNumber negated)
"
- (QDouble fromFloat:1e20) asDoubleArray
+ (QDouble fromFloat:1e20) - 1.0
((QDouble fromFloat:1e20) - (QDouble fromFloat:1.0)) asDoubleArray
(QDouble fromFloat:1e-20) asDoubleArray
((QDouble fromFloat:1e-20) - (QDouble fromFloat:1.0)) asDoubleArray
@@ -1324,6 +3337,34 @@
/ aNumber
"return the quotient of the receiver and the argument, aNumber"
+%{
+ if (__isFloatLike(aNumber)) {
+ double *a = __QDoubleInstPtr(self)->d_qDoubleValue;
+ double b = __floatVal(aNumber);
+ double c0, c1, c2, c3;
+ OBJ newQD;
+ int savedCV;
+
+ fpu_fix_start(&savedCV);
+ qd_div_s(&c0, &c1, &c2, &c3, a[0], a[1], a[2], a[3], b);
+ fpu_fix_end(&savedCV);
+ __qNew_qdReal(newQD, c0, c1, c2, c3);
+ RETURN( newQD );
+ }
+ if (__isQDouble(aNumber)) {
+ double *a = __QDoubleInstPtr(self)->d_qDoubleValue;
+ double *b = __QDoubleInstPtr(aNumber)->d_qDoubleValue;
+ double c0, c1, c2, c3;
+ OBJ newQD;
+ int savedCV;
+
+ fpu_fix_start(&savedCV);
+ qd_div_qd(&c0, &c1, &c2, &c3, a[0], a[1], a[2], a[3], b[0], b[1], b[2], b[3]);
+ fpu_fix_end(&savedCV);
+ __qNew_qdReal(newQD, c0, c1, c2, c3);
+ RETURN( newQD );
+ }
+%}.
^ aNumber quotientFromQDouble:self
"
@@ -1342,14 +3383,14 @@
asDoubleArray
^ DoubleArray
- with:self d0
- with:self d1
- with:self d2
- with:self d3.
+ with:self d0
+ with:self d1
+ with:self d2
+ with:self d3.
"
- (QDouble fromFloat:1.0) asDoubleArray
- (1.0 asQDouble + 1e-40) asDoubleArray
+ (QDouble fromFloat:1.0) asDoubleArray
+ (1.0 asQDouble + 1e-40) asDoubleArray
(QDouble fromFloat:2.0) asDoubleArray
"
@@ -1421,12 +3462,9 @@
"return a QDouble with same value as myself."
^ self
-
-
!
asTrueFraction
-self halt.
^ self d0 asTrueFraction
+ self d1 asTrueFraction
+ self d2 asTrueFraction
@@ -1530,41 +3568,44 @@
"return true, if the argument, aNumber has the same value as than the receiver"
%{
- if (__Class(aNumber) == QDouble) {
- double *a = __QuadDoubleInstPtr(self)->d_quadDoubleValue;
- double *b = __QuadDoubleInstPtr(aNumber)->d_quadDoubleValue;
-
- RETURN ((a[0] == b[0]
- && a[1] == b[1]
- && a[2] == b[2]
- && a[3] == b[3]) ? true : false);
+ if (__isSmallInteger(aNumber)) {
+ double *a = __QDoubleInstPtr(self)->d_qDoubleValue;
+ double b = (double)(__intVal(aNumber));
+
+ RETURN ((a[0] == b
+ && a[1] == 0.0
+ && a[2] == 0.0
+ && a[3] == 0.0) ? true : false);
+ }
+ if (aNumber == nil) {
+ RETURN(false);
}
- if (__Class(aNumber) == Float) {
- double *a = __QuadDoubleInstPtr(self)->d_quadDoubleValue;
- double b = __floatVal(aNumber);
-
- RETURN ((a[0] == b
- && a[1] == 0.0
- && a[2] == 0.0
- && a[3] == 0.0) ? true : false);
+ if (__qClass(aNumber) == QDouble) {
+ double *a = __QDoubleInstPtr(self)->d_qDoubleValue;
+ double *b = __QDoubleInstPtr(aNumber)->d_qDoubleValue;
+
+ RETURN ((a[0] == b[0]
+ && a[1] == b[1]
+ && a[2] == b[2]
+ && a[3] == b[3]) ? true : false);
}
- if (__isSmallInteger(aNumber)) {
- double *a = __QuadDoubleInstPtr(self)->d_quadDoubleValue;
- double b = (double)(__intVal(aNumber));
-
- RETURN ((a[0] == b
- && a[1] == 0.0
- && a[2] == 0.0
- && a[3] == 0.0) ? true : false);
+ if (__qClass(aNumber) == Float) {
+ double *a = __QDoubleInstPtr(self)->d_qDoubleValue;
+ double b = __floatVal(aNumber);
+
+ RETURN ((a[0] == b
+ && a[1] == 0.0
+ && a[2] == 0.0
+ && a[3] == 0.0) ? true : false);
}
%}.
^ aNumber equalFromQDouble:self
"
- 1.0 asQDouble = 1.0 asQDouble
- 1.0 asQDouble = 1.0
- 1.0 asQDouble = 1
- 1.0 asQDouble = 2
+ 1.0 asQDouble = 1.0 asQDouble
+ 1.0 asQDouble = 1.0
+ 1.0 asQDouble = 1
+ 1.0 asQDouble = 2
"
"Created: / 13-06-2017 / 17:12:09 / cg"
@@ -1573,9 +3614,23 @@
!QDouble methodsFor:'double dispatching'!
differenceFromFloat:aFloat
- "aFloat - self"
-
- ^ aFloat + (self negated)
+%{
+ if (__isFloatLike(aFloat)) {
+ double *a = __QDoubleInstPtr(self)->d_qDoubleValue;
+ double b = __floatVal(aFloat);
+ double c0, c1, c2, c3;
+ double e;
+ OBJ newQD;
+ int savedCV;
+
+ fpu_fix_start(&savedCV);
+ s_sub_qd(&c0, &c1, &c2, &c3, b, a[0], a[1], a[2], a[3]);
+ fpu_fix_end(&savedCV);
+ __qNew_qdReal(newQD, c0, c1, c2, c3);
+ RETURN( newQD );
+ }
+%}.
+ ^ super differenceFromFloat:aFloat.
"
1.0 - (QDouble fromFloat:1.0)
@@ -1588,33 +3643,45 @@
(1e20 - (QDouble fromFloat:1.0) + 1e-20) asDoubleArray
"
- "Created: / 12-06-2017 / 23:38:05 / cg"
+
!
differenceFromQDouble:aQDouble
- "aQDouble - self"
-
- ^ aQDouble + (self negated)
+%{
+ if (__isQDouble(aQDouble)) {
+ double *a = __QDoubleInstPtr(aQDouble)->d_qDoubleValue;
+ double *b = __QDoubleInstPtr(self)->d_qDoubleValue;
+ double c0, c1, c2, c3;
+ OBJ newQD;
+ int savedCV;
+
+ fpu_fix_start(&savedCV);
+ qd_sub_qd(&c0, &c1, &c2, &c3, a[0], a[1], a[2], a[3], b[0], b[1], b[2], b[3]);
+ fpu_fix_end(&savedCV);
+ __qNew_qdReal(newQD, c0, c1, c2, c3);
+ RETURN( newQD );
+ }
+%}.
+ ^ super differenceFromQDouble:aQDouble
"
+ (QDouble fromFloat:1.0) - (QDouble fromFloat:1.0)
+ (QDouble fromFloat:1.0) - 1.0
1.0 - (QDouble fromFloat:1.0)
- 1e20 - (QDouble fromFloat:1.0)
- (1.0 - (QDouble fromFloat:1.0)) asFloat
- (1e20 - (QDouble fromFloat:1.0)) asFloat
-
+
+ ((QDouble fromFloat:1.0) - (QDouble fromFloat:1.0)) asDoubleArray
+ ((QDouble fromFloat:1.0) - 1.0) asDoubleArray
(1.0 - (QDouble fromFloat:1.0)) asDoubleArray
+ (1e-20 - (QDouble fromFloat:1.0)) asDoubleArray
(1e20 - (QDouble fromFloat:1.0)) asDoubleArray
- (1e20 - (QDouble fromFloat:1.0) + 1e-20) asDoubleArray
- "
-
- "Created: / 12-06-2017 / 23:38:19 / cg"
+ "
!
equalFromQDouble:aQDouble
%{
if (__Class(aQDouble) == QDouble) {
- double *a = __QuadDoubleInstPtr(self)->d_quadDoubleValue;
- double *b = __QuadDoubleInstPtr(aQDouble)->d_quadDoubleValue;
+ double *a = __QDoubleInstPtr(self)->d_qDoubleValue;
+ double *b = __QDoubleInstPtr(aQDouble)->d_qDoubleValue;
RETURN ((a[0] == b[0]
&& a[1] == b[1]
@@ -1643,8 +3710,8 @@
%{
if (__Class(aQDouble) == QDouble) {
- double *a = __QuadDoubleInstPtr(aQDouble)->d_quadDoubleValue;
- double *b = __QuadDoubleInstPtr(self)->d_quadDoubleValue;
+ double *a = __QDoubleInstPtr(aQDouble)->d_qDoubleValue;
+ double *b = __QDoubleInstPtr(self)->d_qDoubleValue;
// now compare if a < b!
RETURN
@@ -1673,37 +3740,17 @@
productFromFloat:aFloat
%{
if (__isFloatLike(aFloat)) {
- double *a = __QuadDoubleInstPtr(self)->d_quadDoubleValue;
- double b = __floatVal(aFloat);
- double p0, p1, p2, p3;
- double q0, q1, q2;
- double s0, s1, s2, s3, s4;
- OBJ newQD;
- int savedCV;
-
- fpu_fix_start(&savedCV);
-
- m_two_prod(p0, a[0], b, q0);
- m_two_prod(p1, a[1], b, q1);
- m_two_prod(p2, a[2], b, q2);
- p3 = a[3] * b;
-
- s0 = p0;
-
- m_two_sum(s1, q0, p1, s2);
-
- m_three_sum(s2, q1, p2);
-
- m_three_sum2(q1, q2, p3);
- s3 = q1;
-
- s4 = q2 + p2;
-
- m_renorm5(s0, s1, s2, s3, s4);
- fpu_fix_end(&savedCV);
-
- __qNew_qdReal(newQD, s0, s1, s2, s3);
- RETURN( newQD );
+ double a = __floatVal(aFloat);
+ double *b = __QDoubleInstPtr(self)->d_qDoubleValue;
+ double c0, c1, c2, c3;
+ OBJ newQD;
+ int savedCV;
+
+ fpu_fix_start(&savedCV);
+ s_mul_qd(&c0, &c1, &c2, &c3, a, b[0], b[1], b[2], b[3]);
+ fpu_fix_end(&savedCV);
+ __qNew_qdReal(newQD, c0, c1, c2, c3);
+ RETURN( newQD );
}
%}.
^ super productFromFloat:aFloat.
@@ -1711,17 +3758,19 @@
"
loosing bits here:
- (1e20+1.0)*2.0
- (1e20+1.0)*1e20
- (1e40+1.0)*2.0
+ (1e20+1.0)*2.0 - 2E20 -> 0.0
+ (1e20+1.0)*100.0 - 1E+22 -> 0.0
+ (1e20+1.0)*1000.0 - 1E+23 -> 0.0
+ (1e20+1.0)*1e20 - 1E+40 -> 0.0
+ (1e40+1.0)*2.0 - 2E+40 -> 0.0
but not here:
- (1.0 asQDouble) * 2.0
- ((1e20 asQDouble) + (1.0)) * 2.0
- ((1e20 asQDouble) + (1.0)) * 100.0 10000000000000000000100.0
- ((1e20 asQDouble) + (1.0)) * 1000.0
- ((1e40 asQDouble) + (1.0)) * 2.0
+ ((1e20 asQDouble) + (1.0)) * 2.0 - 2E20 -> 2.0
+ ((1e20 asQDouble) + (1.0)) * 100.0 - 1E+22 -> 100.0
+ ((1e20 asQDouble) + (1.0)) * 1000.0 - 1E+23 -> 8389608.0 WRONG
+ ((1e20 asQDouble) + (1.0)) * 1e20 - 1E+40 ->
+ ((1e40 asQDouble) + (1.0)) * 2.0 - 2E+40 ->
2.0 * (QDouble fromFloat:1.0)
2.0 * (QDouble fromFloat:3.0)
@@ -1751,137 +3800,18 @@
productFromQDouble:aQDouble
%{
- if (__Class(aQDouble) == QDouble) {
- double *a = __QuadDoubleInstPtr(self)->d_quadDoubleValue;
- double *b = __QuadDoubleInstPtr(aQDouble)->d_quadDoubleValue;
- OBJ newQD;
- int savedCV;
-
-#define QD_IEEE_MUL
-
-#ifndef QD_IEEE_MUL
- // sloppy
- double p0, p1, p2, p3, p4, p5;
- double q0, q1, q2, q3, q4, q5;
- double t0, t1;
- double s0, s1, s2;
-
- fpu_fix_start(&savedCV);
-
- m_two_prod(p0, a[0], b[0], q0);
- // fprintf(stderr, "%f * %f -> %f, %f\n", a[0], b[0], p0, q0);
-
- m_two_prod(p1, a[0], b[1], q1);
- m_two_prod(p2, a[1], b[0], q2);
- // fprintf(stderr, "%f * %f -> %f, %f\n", a[0], b[1], p1, q1);
- // fprintf(stderr, "%f * %f -> %f, %f\n", a[1], b[0], p2, q2);
-
- m_two_prod(p3, a[0], b[2], q3);
- m_two_prod(p4, a[1], b[1], q4);
- m_two_prod(p5, a[2], b[0], q5);
-
- // fprintf(stderr, "%f * %f -> %f, %f\n", a[0], b[2], p3, q3);
- // fprintf(stderr, "%f * %f -> %f, %f\n", a[1], b[1], p4, q4);
- // fprintf(stderr, "%f * %f -> %f, %f\n", a[2], b[0], p5, q5);
-
- /* Start Accumulation */
- m_three_sum(p1, p2, q0);
- // fprintf(stderr, "mul7: after three_sum: p1:%e p2:%e q0:%e\n", p1, p2, q0);
-
- /* Six-Three Sum of p2, q1, q2, p3, p4, p5. */
- m_three_sum(p2, q1, q2);
- // fprintf(stderr, "mul8: after three_sum: p2:%e q1:%e q2:%e\n", p2, q1, q2);
- m_three_sum(p3, p4, p5);
- // fprintf(stderr, "mul9: after three_sum: p3:%e p4:%e p5:%e\n", p3, p4, p5);
-
- /* compute (s0, s1, s2) = (p2, q1, q2) + (p3, p4, p5). */
- m_two_sum(s0, p2, p3, t0);
- // fprintf(stderr, "mul10: after two_sum: s0:%e p2:%e p3:%e t0:%e\n", s0, p2, p3, t0);
- m_two_sum(s1, q1, p4, t1);
- // fprintf(stderr, "mul11: after two_sum: s1:%e q1:%e p4:%e t1:%e\n", s1, q1, p4, t1);
- s2 = q2 + p5;
- m_two_sum(s1, s1, t0, t0);
- // fprintf(stderr, "mul12: after two_sum: s1:%e s1:%e t0:%e t0:%e\n", s1, s1, t0, t0);
- s2 += (t0 + t1);
-
- /* O(eps^3) order terms */
- s1 += a[0]*b[3] + a[1]*b[2] + a[2]*b[1] + a[3]*b[0] + q0 + q3 + q4 + q5;
-
- // fprintf(stderr, "before renorm5: p0:%e p2:%e s0:%e s1:%e s2:%e\n", p0, p1, s0, s1, s2);
- m_renorm5(p0, p1, s0, s1, s2);
- // fprintf(stderr, "after renorm5: p0:%e p2:%e s0:%e s1:%e s2:%e\n", p0, p1, s0, s1, s2);
- fpu_fix_end(&savedCV);
-
- __qNew_qdReal(newQD, p0, p1, s0, s1);
-#else
- // accurate
- double p0, p1, p2, p3, p4, p5;
- double q0, q1, q2, q3, q4, q5;
- double p6, p7, p8, p9;
- double q6, q7, q8, q9;
- double r0, r1;
- double t0, t1;
- double s0, s1, s2;
-
- fpu_fix_start(&savedCV);
-
- m_two_prod(p0, a[0], b[0], q0);
-
- m_two_prod(p1, a[0], b[1], q1);
- m_two_prod(p2, a[1], b[0], q2);
-
- m_two_prod(p3, a[0], b[2], q3);
- m_two_prod(p4, a[1], b[1], q4);
- m_two_prod(p5, a[2], b[0], q5);
-
- /* Start Accumulation */
- m_three_sum(p1, p2, q0);
-
- /* Six-Three Sum of p2, q1, q2, p3, p4, p5. */
- m_three_sum(p2, q1, q2);
- m_three_sum(p3, p4, p5);
- /* compute (s0, s1, s2) = (p2, q1, q2) + (p3, p4, p5). */
- m_two_sum(s0, p2, p3, t0);
- m_two_sum(s1, q1, p4, t1);
- s2 = q2 + p5;
- m_two_sum(s1, s1, t0, t0);
- s2 += (t0 + t1);
-
- /* O(eps^3) order terms */
- m_two_prod(p6, a[0], b[3], q6);
- m_two_prod(p7, a[1], b[2], q7);
- m_two_prod(p8, a[2], b[1], q8);
- m_two_prod(p9, a[3], b[0], q9);
-
- /* Nine-Two-Sum of q0, s1, q3, q4, q5, p6, p7, p8, p9. */
- m_two_sum(q0, q0, q3, q3);
- m_two_sum(q4, q4, q5, q5);
- m_two_sum(p6, p6, p7, p7);
- m_two_sum(p8, p8, p9, p9);
- /* Compute (t0, t1) = (q0, q3) + (q4, q5). */
- m_two_sum(t0, q0, q4, t1);
- t1 += (q3 + q5);
- /* Compute (r0, r1) = (p6, p7) + (p8, p9). */
- m_two_sum(r0, p6, p8, r1);
- r1 += (p7 + p9);
- /* Compute (q3, q4) = (t0, t1) + (r0, r1). */
- m_two_sum(q3, t0, r0, q4);
- q4 += (t1 + r1);
- /* Compute (t0, t1) = (q3, q4) + s1. */
- m_two_sum(t0, q3, s1, t1);
- t1 += q4;
-
- /* O(eps^4) terms -- Nine-One-Sum */
- t1 += a[1] * b[3] + a[2] * b[2] + a[3] * b[1] + q6 + q7 + q8 + q9 + s2;
-
- // fprintf(stderr, "before accur-renorm5: p0:%e p1:%e s0:%e t0:%e t1:%e\n", p0, p1, s0, t0, t1);
- m_renorm5(p0, p1, s0, t0, t1);
- // fprintf(stderr, "after accur-renorm5: p0:%e p1:%e s0:%e t0:%e t1:%e\n", p0, p1, s0, t0, t1);
- fpu_fix_end(&savedCV);
-
- __qNew_qdReal(newQD, p0, p1, s0, t0);
-#endif
- RETURN( newQD );
+ if (__isQDouble(aQDouble)) {
+ double *a = __QDoubleInstPtr(aQDouble)->d_qDoubleValue;
+ double *b = __QDoubleInstPtr(self)->d_qDoubleValue;
+ double c0, c1, c2, c3;
+ OBJ newQD;
+ int savedCV;
+
+ fpu_fix_start(&savedCV);
+ qd_mul_qd(&c0, &c1, &c2, &c3, a[0], a[1], a[2], a[3], b[0], b[1], b[2], b[3]);
+ fpu_fix_end(&savedCV);
+ __qNew_qdReal(newQD, c0, c1, c2, c3);
+ RETURN( newQD );
}
%}.
^ super productFromQDouble:aQDouble.
@@ -1891,9 +3821,9 @@
2.0 * (QDouble fromFloat:1.0)
(QDouble fromFloat:1.0) * (QDouble fromFloat:2.0)
- 1e20 * (QDouble fromFloat:1.0)
- (2.0 * (QDouble fromFloat:1.0)) asFloat
- (1e20 * (QDouble fromFloat:1.0)) asFloat
+ 1e20 * (QDouble fromFloat:2.0)
+ 2.0 * (QDouble fromFloat:1e20)
+ (QDouble fromFloat:1e20) * (QDouble fromFloat:1e20)
(1e20 * (QDouble fromFloat:1.0) * 1e-20) asDoubleArray
@@ -1904,180 +3834,29 @@
"Modified: / 05-07-2017 / 11:07:16 / cg"
!
-quotientFromQDouble:aQDouble
- "sloppy"
-
- |q0 q1 q2 q3 r|
-
- q0 := aQDouble d0 / self d0.
- "/ Stdout showCR:('q0: %1 (a[0]=%2; b[0]=%3)\n' bindWith:q0 with:self d0 with:aQDouble d0).
- r := aQDouble - (self * q0).
- "/ Stdout showCR:('r: %1\n' bindWith:r asDoubleArray).
-
- q1 := r d0 / self d0.
- "/ Stdout showCR:('q1: %1 (r[0]=%2; b[0]=%3)\n' bindWith:q1 with:r d0 with:aQDouble d0).
- r := r - (self * q1).
- "/ Stdout showCR:('r: %1\n' bindWith:r asDoubleArray).
-
- q2 := r d0 / self d0.
- "/ Stdout showCR:('q2: %1 (r[0]=%2; b[0]=%3)\n' bindWith:q2 with:r d0 with:aQDouble d0).
- r := r - (self * q2).
- "/ Stdout showCR:('r: %1\n' bindWith:r asDoubleArray).
-
- q3 := r d0 / self d0.
- "/ Stdout showCR:('q3: %1 (r[0]=%2; b[0]=%3)\n' bindWith:q3 with:r d0 with:aQDouble d0).
-
- r := self class d0:q0 d1:q1 d2:q2 d3:q3.
- "/ Stdout showCR:('before renorm: %1\n' bindWith:r asDoubleArray).
- r renorm.
- "/ Stdout showCR:('after renorm: %1\n' bindWith:r asDoubleArray).
- ^ r
+quotientFromFloat:aFloat
+%{
+ if (__isFloatLike(aFloat)) {
+ double a = __floatVal(aFloat);
+ double *b = __QDoubleInstPtr(self)->d_qDoubleValue;
+ double c0, c1, c2, c3;
+ OBJ newQD;
+ int savedCV;
+
+ fpu_fix_start(&savedCV);
+ s_div_qd(&c0, &c1, &c2, &c3, a, b[0], b[1], b[2], b[3]);
+ fpu_fix_end(&savedCV);
+ __qNew_qdReal(newQD, c0, c1, c2, c3);
+ RETURN( newQD );
+ }
+%}.
+ ^ super quotientFromFloat:aFloat.
"
2.0 / (QDouble fromFloat:2.0)
2.0 / (QDouble fromFloat:1.0)
1e20 / (QDouble fromFloat:1.0)
- (2.0 / (QDouble fromFloat:1.0)) asFloat
- (1e20 / (QDouble fromFloat:1.0)) asFloat
-
- (QDouble fromFloat:2.0) / 2.0
- (QDouble fromFloat:1e20) / 2.0
- ((QDouble fromFloat:1.0) / 2.0) asFloat
- ((QDouble fromFloat:1e20 / 2.0)) asFloat
-
- ((1e20 + (QDouble fromFloat:1.0) + 1e-20) / 2.0) asDoubleArray
-
- ((QDouble fromFloat:10.0) quotientFromQDouble: (QDouble fromFloat:1.234)) asDoubleArray
- ((QDouble fromFloat:1.234) / (QDouble fromFloat:10.0)) asDoubleArray
-
-q0: 1.234000e-01 (a[0]=1.234000e+00; b[0]=1.000000e+01)
-a: 1.234000e+00/0.000000e+00/0.000000e+00/0.000000e+00
-b: 1.000000e+01/0.000000e+00/0.000000e+00/0.000000e+00
-(b * q0): 1.234000e+00/-2.775558e-17/0.000000e+00/0.000000e+00
-r: 2.775558e-17/0.000000e+00/0.000000e+00/0.000000e+00
-
-q1: 2.775558e-18 (r[0]=2.775558e-17; b[0]=1.000000e+01)
-r: -1.540744e-33/0.000000e+00/0.000000e+00/0.000000e+00
-
-q2: -1.540744e-34 (r[0]=-1.540744e-33; b[0]=1.000000e+01)
-r: 8.552847e-50/0.000000e+00/0.000000e+00/0.000000e+00
-
-q3: 8.552847e-51 (r[0]=8.552847e-50; b[0]=1.000000e+01)
-
-before renorm: 1.234000e-01/2.775558e-18/-1.540744e-34/8.552847e-51
-after renorm: 1.234000e-01/2.775558e-18/-1.540744e-34/8.552847e-51
-1.234/10.0 is: 0.123400 / 0.000000 / -0.000000 / 0.000000
-
-
- "
-
- "Created: / 13-06-2017 / 17:50:35 / cg"
- "Modified (comment): / 15-06-2017 / 01:02:05 / cg"
-!
-
-quotientFromQDouble_accurate:aQDouble
- |q0 q1 q2 q3 q4 r|
-
-%{
-/* Computes fl(a+b) and err(a+b). Assumes |a| >= |b|. */
-#define quick_two_sum(s, a, b, err) \
- { \
- s = a + b; \
- err = b - (s - a); \
- }
-
-#define renorm5(c0, c1, c2, c3, c4) \
- { \
- double s0, s1, s2 = 0.0, s3 = 0.0; \
- \
- quick_two_sum(s0, c3, c4, c4); \
- quick_two_sum(s0, c2, s0, c3); \
- quick_two_sum(s0, c1, s0, c2); \
- quick_two_sum(c0, c0, s0, c1); \
- \
- s0 = c0; \
- s1 = c1; \
- \
- quick_two_sum(s0, c0, c1, s1); \
- if (s1 != 0.0) { \
- quick_two_sum(s1, s1, c2, s2); \
- if (s2 != 0.0) { \
- quick_two_sum(s2, s2, c3, s3); \
- if (s3 != 0.0) { \
- s3 += c4; \
- } else { \
- s2 += c4; \
- } \
- } else { \
- quick_two_sum(s1, s1, c3, s2); \
- if (s2 != 0.0) { \
- quick_two_sum(s2, s2, c4, s3); \
- } else { \
- quick_two_sum(s1, s1, c4, s2); \
- } \
- } \
- } else { \
- quick_two_sum(s0, s0, c2, s1); \
- if (s1 != 0.0) { \
- quick_two_sum(s1, s1, c3, s2); \
- if (s2 != 0.0) { \
- quick_two_sum(s2, s2, c4, s3); \
- } else { \
- quick_two_sum(s1, s1, c4, s2); \
- } \
- } else { \
- quick_two_sum(s0, s0, c3, s1); \
- if (s1 != 0.0) { \
- quick_two_sum(s1, s1, c4, s2); \
- } else { \
- quick_two_sum(s0, s0, c4, s1); \
- } \
- } \
- } \
- \
- c0 = s0; \
- c1 = s1; \
- c2 = s2; \
- c3 = s3; \
- }
-%}.
- q0 := aQDouble d0 / self d0.
- r := aQDouble - (self * q0).
-
- q1 := r d0 / self d0.
- r := r - (self * q1).
-
- q2 := r d0 / self d0.
- r := r - (self * q2).
-
- q3 := r d0 / self d0.
- r := r - (self * q3).
-
- q4 := r d0 / self d0.
- q0 isFinite ifTrue:[
-%{
- double cq0, cq1, cq2, cq3, cq4;
- cq0 = __floatVal(q0);
- cq1 = __floatVal(q1);
- cq2 = __floatVal(q2);
- cq3 = __floatVal(q3);
- cq4 = __floatVal(q4);
- renorm5(cq0, cq1, cq2, cq3, cq4)
- {
- OBJ newQD;
- __qNew_qdReal(newQD, cq0, cq1, cq2, cq3);
- RETURN (newQD);
- }
-#undef quick_two_sum
-#undef renorm5
-%}.
- ].
- ^ self class d0:q0 d1:0.0 d2:0.0 d3:0.0
-
- "
- 2.0 / (QDouble fromFloat:2.0)
- 2.0 / (QDouble fromFloat:1.0)
- 1e20 / (QDouble fromFloat:1.0)
+ 1e20 / (QDouble fromFloat:2.0)
(2.0 / (QDouble fromFloat:1.0)) asFloat
(1e20 / (QDouble fromFloat:1.0)) asFloat
@@ -2091,41 +3870,34 @@
((QDouble fromFloat:10.0) quotientFromQDouble: (QDouble fromFloat:1.234)) asDoubleArray
((QDouble fromFloat:1.234) / (QDouble fromFloat:10.0)) asDoubleArray
"
+
+ "Created: / 13-06-2017 / 17:50:35 / cg"
+ "Modified (comment): / 15-06-2017 / 01:02:05 / cg"
!
-quotientFromQDouble_sloppy:aQDouble
- "sloppy"
-
- |q0 q1 q2 q3 r|
-
- q0 := aQDouble d0 / self d0.
- "/ Stdout showCR:('q0: %1 (a[0]=%2; b[0]=%3)\n' bindWith:q0 with:self d0 with:aQDouble d0).
- r := aQDouble - (self * q0).
- "/ Stdout showCR:('r: %1\n' bindWith:r asDoubleArray).
-
- q1 := r d0 / self d0.
- "/ Stdout showCR:('q1: %1 (r[0]=%2; b[0]=%3)\n' bindWith:q1 with:r d0 with:aQDouble d0).
- r := r - (self * q1).
- "/ Stdout showCR:('r: %1\n' bindWith:r asDoubleArray).
-
- q2 := r d0 / self d0.
- "/ Stdout showCR:('q2: %1 (r[0]=%2; b[0]=%3)\n' bindWith:q2 with:r d0 with:aQDouble d0).
- r := r - (self * q2).
- "/ Stdout showCR:('r: %1\n' bindWith:r asDoubleArray).
-
- q3 := r d0 / self d0.
- "/ Stdout showCR:('q3: %1 (r[0]=%2; b[0]=%3)\n' bindWith:q3 with:r d0 with:aQDouble d0).
-
- r := self class d0:q0 d1:q1 d2:q2 d3:q3.
- "/ Stdout showCR:('before renorm: %1\n' bindWith:r asDoubleArray).
- r renorm.
- "/ Stdout showCR:('after renorm: %1\n' bindWith:r asDoubleArray).
- ^ r
+quotientFromQDouble:aQDouble
+%{
+ if (__isQDouble(aQDouble)) {
+ double *a = __QDoubleInstPtr(aQDouble)->d_qDoubleValue;
+ double *b = __QDoubleInstPtr(self)->d_qDoubleValue;
+ double c0, c1, c2, c3;
+ OBJ newQD;
+ int savedCV;
+
+ fpu_fix_start(&savedCV);
+ qd_div_qd(&c0, &c1, &c2, &c3, a[0], a[1], a[2], a[3], b[0], b[1], b[2], b[3]);
+ fpu_fix_end(&savedCV);
+ __qNew_qdReal(newQD, c0, c1, c2, c3);
+ RETURN( newQD );
+ }
+%}.
+ ^ super quotientFromQDouble:aQDouble.
"
2.0 / (QDouble fromFloat:2.0)
2.0 / (QDouble fromFloat:1.0)
1e20 / (QDouble fromFloat:1.0)
+ 1e20 / (QDouble fromFloat:2.0)
(2.0 / (QDouble fromFloat:1.0)) asFloat
(1e20 / (QDouble fromFloat:1.0)) asFloat
@@ -2138,26 +3910,6 @@
((QDouble fromFloat:10.0) quotientFromQDouble: (QDouble fromFloat:1.234)) asDoubleArray
((QDouble fromFloat:1.234) / (QDouble fromFloat:10.0)) asDoubleArray
-
-q0: 1.234000e-01 (a[0]=1.234000e+00; b[0]=1.000000e+01)
-a: 1.234000e+00/0.000000e+00/0.000000e+00/0.000000e+00
-b: 1.000000e+01/0.000000e+00/0.000000e+00/0.000000e+00
-(b * q0): 1.234000e+00/-2.775558e-17/0.000000e+00/0.000000e+00
-r: 2.775558e-17/0.000000e+00/0.000000e+00/0.000000e+00
-
-q1: 2.775558e-18 (r[0]=2.775558e-17; b[0]=1.000000e+01)
-r: -1.540744e-33/0.000000e+00/0.000000e+00/0.000000e+00
-
-q2: -1.540744e-34 (r[0]=-1.540744e-33; b[0]=1.000000e+01)
-r: 8.552847e-50/0.000000e+00/0.000000e+00/0.000000e+00
-
-q3: 8.552847e-51 (r[0]=8.552847e-50; b[0]=1.000000e+01)
-
-before renorm: 1.234000e-01/2.775558e-18/-1.540744e-34/8.552847e-51
-after renorm: 1.234000e-01/2.775558e-18/-1.540744e-34/8.552847e-51
-1.234/10.0 is: 0.123400 / 0.000000 / -0.000000 / 0.000000
-
-
"
"Created: / 13-06-2017 / 17:50:35 / cg"
@@ -2167,7 +3919,7 @@
sumFromFloat:aFloat
%{
if (__isFloatLike(aFloat)) {
- double *a = __QuadDoubleInstPtr(self)->d_quadDoubleValue;
+ double *a = __QDoubleInstPtr(self)->d_qDoubleValue;
double b = __floatVal(aFloat);
double c0, c1, c2, c3;
double e;
@@ -2175,13 +3927,7 @@
int savedCV;
fpu_fix_start(&savedCV);
-
- m_two_sum(c0, a[0], b, e);
- m_two_sum(c1 ,a[1], e, e);
- m_two_sum(c2, a[2], e, e);
- m_two_sum(c3, a[3], e, e);
-
- m_renorm5(c0, c1, c2, c3, e);
+ qd_add_s(&c0, &c1, &c2, &c3, a[0], a[1], a[2], a[3], b);
fpu_fix_end(&savedCV);
__qNew_qdReal(newQD, c0, c1, c2, c3);
RETURN( newQD );
@@ -2219,150 +3965,18 @@
sumFromQDouble:aQDouble
%{
-
-/* s = quick_three_accum(a, b, c) adds c to the dd-pair (a, b).
- * If the result does not fit in two doubles, then the sum is
- * output into s and (a,b) contains the remainder. Otherwise
- * s is zero and (a,b) contains the sum. */
-# define m_quick_three_accum(q3_outRef, q3_a, q3_b, q3_c) \
-{ \
- double s0, s;\
- int za, zb;\
-\
- m_two_sum(s0, q3_b, q3_c, q3_b); \
- m_two_sum(s, q3_a, s0, q3_a); \
-\
- za = (q3_a != 0.0);\
- zb = (q3_b != 0.0);\
-\
- if (za && zb) {\
- q3_outRef = s;\
- } else {\
- if (!zb) {\
- q3_b = q3_a;\
- q3_a = s;\
- } else {\
- q3_a = s;\
- }\
- q3_outRef = 0.0;\
- }\
-}
-
- if (__Class(aQDouble) == QDouble) {
- double *a = __QuadDoubleInstPtr(self)->d_quadDoubleValue;
- double *b = __QuadDoubleInstPtr(aQDouble)->d_quadDoubleValue;
- OBJ newQD;
- int savedCV;
-
-#define SLOPPY_ADD
-#ifdef SLOPPY_ADD
- /* Addition re-organized to minimize
- data dependency ... unfortunately some compilers are
- not very smart to do this automatically */
- double s0, s1, s2, s3;
- double t0, t1, t2, t3;
- double v0, v1, v2, v3;
- double u0, u1, u2, u3;
- double w0, w1, w2, w3;
-
- fpu_fix_start(&savedCV);
- s0 = a[0] + b[0];
- s1 = a[1] + b[1];
- s2 = a[2] + b[2];
- s3 = a[3] + b[3];
-
- v0 = s0 - a[0];
- v1 = s1 - a[1];
- v2 = s2 - a[2];
- v3 = s3 - a[3];
-
- u0 = s0 - v0;
- u1 = s1 - v1;
- u2 = s2 - v2;
- u3 = s3 - v3;
-
- w0 = a[0] - u0;
- w1 = a[1] - u1;
- w2 = a[2] - u2;
- w3 = a[3] - u3;
-
- u0 = b[0] - v0;
- u1 = b[1] - v1;
- u2 = b[2] - v2;
- u3 = b[3] - v3;
-
- t0 = w0 + u0;
- t1 = w1 + u1;
- t2 = w2 + u2;
- t3 = w3 + u3;
-
- m_two_sum(s1, s1, t0, t0);
- m_three_sum(s2, t0, t1);
- m_three_sum2(s3, t0, t2);
- t0 = t0 + t1 + t3;
-
- /* renormalize */
- m_renorm5(s0, s1, s2, s3, t0);
- fpu_fix_end(&savedCV);
- __qNew_qdReal(newQD, s0, s1, s2, s3);
-#else
- // ieee_add...
- int i, j, k;
- double s, t;
- double u, v; /* double-length accumulator */
- double x[4] = {0.0, 0.0, 0.0, 0.0};
-
- fpu_fix_start(&savedCV);
-
- i = j = k = 0;
- if (abs(a[i]) > abs(b[j]))
- u = a[i++];
- else
- u = b[j++];
- if (abs(a[i]) > abs(b[j]))
- v = a[i++];
- else
- v = b[j++];
-
- m_quick_two_sum(u, u, v, v);
-
- while (k < 4) {
- if (i >= 4 && j >= 4) {
- x[k] = u;
- if (k < 3) {
- x[++k] = v;
- }
- break;
- }
-
- if (i >= 4) {
- t = b[j++];
- } else if (j >= 4) {
- t = a[i++];
- } else if (abs(a[i]) > abs(b[j])) {
- t = a[i++];
- } else {
- t = b[j++];
- }
- m_quick_three_accum(s, u, v, t);
-
- if (s != 0.0) {
- x[k++] = s;
- }
- }
-
- /* add the rest. */
- for (k = i; k < 4; k++) {
- x[3] += a[k];
- }
- for (k = j; k < 4; k++) {
- x[3] += b[k];
- }
- m_renorm4(x[0], x[1], x[2], x[3]);
- fpu_fix_end(&savedCV);
- __qNew_qdReal(newQD, x[0], x[1], x[2], x[3]);
-#endif
- RETURN(newQD);
+ if (__isQDouble(aQDouble)) {
+ double *a = __QDoubleInstPtr(self)->d_qDoubleValue;
+ double *b = __QDoubleInstPtr(aQDouble)->d_qDoubleValue;
+ double c0, c1, c2, c3;
+ OBJ newQD;
+ int savedCV;
+
+ fpu_fix_start(&savedCV);
+ qd_add_qd(&c0, &c1, &c2, &c3, a[0], a[1], a[2], a[3], b[0], b[1], b[2], b[3]);
+ fpu_fix_end(&savedCV);
+ __qNew_qdReal(newQD, c0, c1, c2, c3);
+ RETURN( newQD );
}
%}.
^ super sumFromQDouble:aQDouble
@@ -2398,181 +4012,46 @@
!QDouble methodsFor:'mathematical functions'!
-exp
- "/ the exp_sloppy code is the algorithm from the original C++ qd package;
- "/ however, it is inexact in the 37th digit
- "/ Therefore, use the inherited code, which is slow, but more precise.
-
- ^ super exp
-
- "
- 2.0 exp -> 7.38905609893065
- 2.0 asQDouble exp -> 7.38905609893065022723
-
- 2 asQDouble exp printfPrintString:'%70.68f'
- -> 7.389056098930650227230427460575007813180315570551847324087127822432669
-
- actual result:
- -> 7.3890560989306502272304274605750078131803155705518473240871278225225737960790577633843124850791217947737531612654...
+cos
+
+%{
+ double *a = __QDoubleInstPtr(self)->d_qDoubleValue;
+ double q0, q1, q2, q3;
+ OBJ newQD;
+ int savedCV;
+
+ fpu_fix_start(&savedCV);
+ qdcos(&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 );
+%}.
"
-
- "Created: / 19-06-2017 / 01:49:32 / cg"
- "Modified (comment): / 22-06-2017 / 14:32:53 / cg"
+ 1.0 cos
+ (QDouble fromFloat:1.0) cos
+ "
!
-exp_sloppy
- "/ this is the algorithm from the qd C++ package;
- "/ however, it is inexact in the 37th digit
- "/
- "/ the inherited exp code is much slower, but more precise.
-
- "/ Strategy: We first reduce the size of x by noting that
- "/
- "/ exp(kr + m * log(2)) = 2^m * exp(r)^k
- "/
- "/ where m and k are integers. By choosing m appropriately
- "/ we can make |kr| <= log(2) / 2 = 0.347. Then exp(r) is
- "/ evaluated using the familiar Taylor series. Reducing the
- "/ argument substantially speeds up the convergence. */
-
- |k inv_k d0 m mul_pwr2 mul2 r
- s p t thresh eps i|
-
- eps := 1.21543267145725e-63. "/ = 2^-209.
-
- d0 := self d0.
- (d0 <= -709.0) ifTrue:[
- ^ 0.0 asQDouble.
- ].
- (d0 >= 709.0) ifTrue:[
- ^ Infinity positive
- ].
-
- (d0 = 0.0) ifTrue:[
- ^ 1.0 asQDouble.
- ].
- (d0 = 1.0) ifTrue:[
- self isOne ifTrue:[
- ^ self class e
- ].
- ].
-
- k := 65536.0. "/ 1.0 ldexp:16.
- inv_k := 1.0 / k.
-
- m := (d0 / Float ln2 + 0.5) floor.
-
- mul_pwr2 := [:a :b | QDouble d0:(a d0 * b) d1:(a d1 * b) d2:(a d2 * b) d3:(a d3 * b) ].
- mul2 := [:a | QDouble d0:(a d0 * 2.0) d1:(a d1 * 2.0) d2:(a d2 * 2.0) d3:(a d3 * 2.0) ].
-
- r := mul_pwr2 value:(self - (self class ln2 * m)) value:inv_k.
- thresh := inv_k * eps.
-
- p := r squared.
- s := r + (mul_pwr2 value:p value:0.5).
- i := 1.
- [
- p := p * r.
- t := p * (self class invFact at:i).
- i := i+1.
- s := s + t.
- ] doWhile:[ (t asFloat abs > thresh) and:[i < 10] ].
-
- s := (mul2 value:s) + s squared.
- s := (mul2 value:s) + s squared.
- s := (mul2 value:s) + s squared.
- s := (mul2 value:s) + s squared.
- s := (mul2 value:s) + s squared.
- s := (mul2 value:s) + s squared.
- s := (mul2 value:s) + s squared.
- s := (mul2 value:s) + s squared.
- s := (mul2 value:s) + s squared.
- s := (mul2 value:s) + s squared.
- s := (mul2 value:s) + s squared.
- s := (mul2 value:s) + s squared.
- s := (mul2 value:s) + s squared.
- s := (mul2 value:s) + s squared.
- s := (mul2 value:s) + s squared.
- s := (mul2 value:s) + s squared.
-
- s := s + 1.0.
- ^ s ldexp:m asInteger.
+exp
+
+%{
+ double *a = __QDoubleInstPtr(self)->d_qDoubleValue;
+ double q0, q1, q2, q3;
+ OBJ newQD;
+ int savedCV;
+
+ fpu_fix_start(&savedCV);
+ qdexp(&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 );
+%}.
"
- 1.0 exp -> 2.71828182845905
- 1.0 asQDouble exp -> 2.718281828459045235360287471353
-
- 10.0 exp -> 22026.4657948067
- 10.0 asQDouble exp -> 22026.46579480671651695790064528
-
- Wolfram:
- 22026.46579480671651695790064528424436635351261855678107423...
-
- 1000.0 exp -> INF
- 1000.0 asQDouble exp -> 22026.46579480671651695790064528
+ 1.0 exp
+ (QDouble fromFloat:1.0) exp
"
-
- "Created: / 21-06-2017 / 13:48:27 / cg"
-!
-
-exp_withAccuracy:epsilon
- "compute e^x of the receiver using a taylor series approximation.
- This method is only invoked for limitedPrecisionReal classes, which do not compute
- exp themself (i.e. QDouble)"
-
- "/ uses taylor series:
- "/ x x^2 x^3
- "/ e^x = 1 + --- + --- + --- ...
- "/ 1!! 2!! 3!!
-
- |x2 n num den approx delta|
-
- x2 := self asLimitedPrecisionReal squared.
-
- num := x2.
- den := self coerce:2.
- n := 3.
- approx := self + 1 + (num / den).
-
- [
- n := n + 1.
- den := den * n.
- num := num * self.
-
- delta := num / den.
- "/ Transcript showCR:delta.
- delta isNaN ifTrue:[self halt:'nan when dividing for delta'. num / den].
- delta = 0 ifTrue:[self halt:'zero delta'].
- approx := approx + delta.
- ] doUntil:[delta abs <= epsilon].
-
- ^ approx
-
- "
- wolfram:
- 7.389056098930650227230427460575007813180315570551847324087
-
- 2 asQDouble exp_withAccuracy:QDouble epsilon 7.38905609893065022723
- 2 asQDouble exp_withAccuracy:LongFloat epsilon 7.38905609893065022723
- 2 asQDouble exp_withAccuracy:Float epsilon 7.38905609893065022489
-
- 2.718281828459045235360287471352662497757247093699959574966...
-
- 1 asQDouble exp_withAccuracy:QDouble epsilon 2.71828182845904523536
- 1 asQDouble exp_withAccuracy:LongFloat epsilon 2.71828182845904523536
- 1 asQDouble exp_withAccuracy:Float epsilon 2.71828182845904522671
-
- 1.7392749415205010473946813036112352261479840577250084... × 10^18
-
- 42 asQDouble exp_withAccuracy:QDouble epsilon NAN
- 42 asQDouble exp_withAccuracy:LongFloat epsilon 1.73927494152050104739e18
- 42 asQDouble exp_withAccuracy:Float epsilon 1.73927494152050104739e18
-
- "
-
- "Created: / 04-07-2017 / 11:58:26 / cg"
- "Modified: / 10-10-2017 / 16:03:34 / cg"
!
ldexp:exp
@@ -2582,13 +4061,13 @@
mantissa and exponent: (f mantissa ldexp:f exponent) = f"
^ self class
- d0:(self d0 ldexp:exp)
- d1:(self d1 ldexp:exp)
- d2:(self d2 ldexp:exp)
- d3:(self d3 ldexp:exp)
+ d0:(self d0 ldexp:exp)
+ d1:(self d1 ldexp:exp)
+ d2:(self d2 ldexp:exp)
+ d3:(self d3 ldexp:exp)
"
- |f| f := 1 asQDouble. (f mantissa ldexp:f exponent) -> 1.0
- |f| f := (1e40 asQDouble + 1e-40). (f mantissa ldexp:f exponent) -> (1e40 asQDouble + 1e-40)
+ |f| f := 1 asQDouble. (f mantissa ldexp:f exponent) -> 1.0
+ |f| f := (1e40 asQDouble + 1e-40). (f mantissa ldexp:f exponent) -> (1e40 asQDouble + 1e-40)
1.0 ldexp:16 -> 65536.0
1.0 asQDouble ldexp:16 -> 65536.0
@@ -2612,17 +4091,17 @@
d0 := self d0.
d0 = 1.0 ifTrue:[
- self isOne ifTrue:[ ^ self class zero ].
+ self isOne ifTrue:[ ^ self class zero ].
].
d0 > 0.0 ifTrue:[
- "/ initial approx.
- x := d0 ln asQDouble.
- "/ three more iterations of newton...
- x := x + (self / (x exp)) - 1.0.
- x := x + (self / (x exp)) - 1.0.
- x := x + (self / (x exp)) - 1.0.
-
- ^ x
+ "/ initial approx.
+ x := d0 ln asQDouble.
+ "/ three more iterations of newton...
+ x := x + (self / (x exp)) - 1.0.
+ x := x + (self / (x exp)) - 1.0.
+ x := x + (self / (x exp)) - 1.0.
+
+ ^ x
].
"/ now done via trapInfinity; was:
@@ -2632,11 +4111,11 @@
"/ if you need -INF for a zero receiver, try Number trapInfinity:[...]
^ self class
- raise:(self = 0 ifTrue:[#infiniteResultSignal] ifFalse:[#domainErrorSignal])
- receiver:self
- selector:#ln
- arguments:#()
- errorString:'bad receiver in ln (not strictly positive)'
+ raise:(self = 0 ifTrue:[#infiniteResultSignal] ifFalse:[#domainErrorSignal])
+ receiver:self
+ selector:#ln
+ arguments:#()
+ errorString:'bad receiver in ln (not strictly positive)'
"
-1 ln
@@ -2646,22 +4125,22 @@
1.0 asQDouble ln
3.0 ln printfPrintString:'%60.58lf'
- -> 1.0986122886681097821082175869378261268138885498046875000000'
- ^
+ -> 1.0986122886681097821082175869378261268138885498046875000000'
+ ^
3.0 asQDouble ln printfPrintString:'%60.58f'
- -> 1.0986122886681096913952452369225257046474905578227494517347
+ -> 1.0986122886681096913952452369225257046474905578227494517347
3.0 asQDouble ln printfPrintString:'%70.68f'
- -> 1.09861228866810969139524523692252570464749055782274945173469433364779
+ -> 1.09861228866810969139524523692252570464749055782274945173469433364779
(3.0 asQDouble ln_withAccuracy:1e-64) printfPrintString:'%70.68f'
- 1.09861228866810969139524523692252570464749055782274945173469433364475
+ 1.09861228866810969139524523692252570464749055782274945173469433364475
(3.0 asQDouble ln_withAccuracy:1e-100) printfPrintString:'%70.68f'
- '1.098612288668109691395245236922525704647490557822749451734694333656909'
+ '1.098612288668109691395245236922525704647490557822749451734694333656909'
actual result:
- -> 1.0986122886681096913952452369225257046474905578227494517346943336374942932186089668736157548137320887879700290659...
+ -> 1.0986122886681096913952452369225257046474905578227494517346943336374942932186089668736157548137320887879700290659...
"
"Created: / 18-06-2017 / 23:32:54 / cg"
@@ -2687,93 +4166,89 @@
"Modified (comment): / 12-06-2017 / 23:46:57 / cg"
!
+raisedToInteger:n
+
+%{
+ if (__isSmallInteger(n)) {
+ double *a = __QDoubleInstPtr(self)->d_qDoubleValue;
+ double q0, q1, q2, q3;
+ OBJ newQD;
+ int savedCV;
+
+ fpu_fix_start(&savedCV);
+ qdpow(&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 );
+ }
+%}.
+ ^ super raisedToInteger:n.
+
+ "
+ (QDouble fromFloat:4.0) raisedToInteger:4
+ (QDouble fromFloat:10.0) raisedToInteger:10
+ (QDouble fromFloat:10.0000000000001) raisedToInteger:10
+ 10.0000000000001 raisedToInteger:10
+ "
+!
+
+sin
+
+%{
+ double *a = __QDoubleInstPtr(self)->d_qDoubleValue;
+ double q0, q1, q2, q3;
+ OBJ newQD;
+ int savedCV;
+
+ fpu_fix_start(&savedCV);
+ qdsin(&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 );
+%}.
+
+ "
+ 1.0 sin
+ (QDouble fromFloat:1.0) sin
+ "
+!
+
sqrt
"Return the square root of the receiver"
- "this computes a roughly 65 digits precision result,
- using sqrt from the double as an initial guess"
-
- |guess|
-
- guess := self d0 sqrt asQDouble.
- ^ self sqrt_withAccuracy:(self epsilon) fromInitialGuess:guess
-
- "FloatD is only correct in roughly 17 digits
-
- 2 sqrt printfPrintString:'%50.48f'
- -> 1.414213562373095145474621858738828450441360473633
- ^
- QDouble gives you much more:
-
- 2 asQDouble sqrt printfPrintString:'%50.48f'
- -> 1.414213562373095048801688724209698078569671875377
-
- 2 asQDouble sqrt printfPrintString:'%60.58f'
- -> 1.4142135623730950488016887242096980785696718753769480731767'
-
- 2 asQDouble sqrt printfPrintString:'%70.68f'
- -> 1.41421356237309504880168872420969807856967187537694807317667973799602
-
- actual digits (wolfram):
- -> 1.4142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276415727350138462309...
- "
-
- "Created: / 22-06-2017 / 13:53:47 / cg"
- "Modified: / 22-06-2017 / 17:08:06 / mawalch"
- "Modified (format): / 03-07-2017 / 12:09:04 / cg"
+%{
+ double *a = __QDoubleInstPtr(self)->d_qDoubleValue;
+ double q0, q1, q2, q3;
+ OBJ newQD;
+ int savedCV;
+
+ fpu_fix_start(&savedCV);
+ qdsqrt(&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 );
+%}.
+
+ "
+ (QDouble fromFloat:4.0) sqrt
+ (QDouble fromFloat:2.0) sqrt
+ (QDouble fromFloat:1e20) sqrt
+ "
!
squared
"return receiver * receiver"
%{
- double *a = __QuadDoubleInstPtr(self)->d_quadDoubleValue;
- double p0, p1, p2, p3, p4, p5;
+ double *a = __QDoubleInstPtr(self)->d_qDoubleValue;
double q0, q1, q2, q3;
- double s0, s1;
- double t0, t1;
- double t;
OBJ newQD;
-
- // code makes use of some symetries to avoid a few multiplications...
-
- m_two_sqr(p0, a[0], q0);
- t = 2.0 * a[0];
-
- m_two_prod(p1, t, a[1], q1);
- m_two_prod(p2, t, a[2], q2);
- m_two_sqr(p3, a[1], q3);
-
- m_two_sum(p1, q0, p1, q0);
-
- m_two_sum(q0, q0, q1, q1);
- m_two_sum(p2, p2, p3, p3);
-
- m_two_sum(s0, q0, p2, t0);
- m_two_sum(s1, q1, p3, t1);
-
- m_two_sum(s1, s1, t0, t0);
- t0 += t1;
-
- m_quick_two_sum(s1, s1, t0, t0);
- m_quick_two_sum(p2, s0, s1, t1);
- m_quick_two_sum(p3, t1, t0, q0);
-
- p4 = 2.0 * a[0] * a[3];
- p5 = 2.0 * a[1] * a[2];
-
- m_two_sum(p4, p4, p5, p5);
- m_two_sum(q2, q2, q3, q3);
-
- m_two_sum(t0, p4, q2, t1);
- t1 = t1 + p5 + q3;
-
- m_two_sum(p3, p3, t0, p4);
- p4 = p4 + q0 + t1;
-
- m_renorm5(p0, p1, p2, p3, p4);
-
- __qNew_qdReal(newQD, p0, p1, s0, s1);
+ int savedCV;
+
+ fpu_fix_start(&savedCV);
+ qd_sqr(&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 );
%}.
@@ -2784,6 +4259,27 @@
"Created: / 13-06-2017 / 01:27:58 / cg"
"Modified: / 22-06-2017 / 14:08:31 / cg"
+!
+
+tan
+
+%{
+ double *a = __QDoubleInstPtr(self)->d_qDoubleValue;
+ double q0, q1, q2, q3;
+ OBJ newQD;
+ int savedCV;
+
+ fpu_fix_start(&savedCV);
+ qdtan(&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 );
+%}.
+
+ "
+ 1.0 tan
+ (QDouble fromFloat:1.0) tan
+ "
! !
!QDouble methodsFor:'printing & storing'!
@@ -2924,20 +4420,26 @@
"return a printed representation of the receiver.
Notice:
- this code was adapted from an ugly piece of c++ code,
- which was obviously hacked.
- It does need a rework.
- As an alternative, use the printf functions, which should also deal wth QDoubles"
+ this code was adapted from an ugly piece of c++ code,
+ which was obviously hacked.
+ It does need a rework.
+ As an alternative, use the printf functions, which should also deal wth QDoubles"
self d1 = 0.0 ifTrue:[
- self d0 printOn:aStream.
- ^ self
+ self d0 printOn:aStream.
+ ^ self
].
thisContext isRecursive ifTrue:[
- aStream nextPutAll:'aQDouble (error while printing)'.
- ^ self.
+ aStream nextPutAll:'aQDouble (error while printing)'.
+ ^ self.
].
+ self d0 printOn:aStream. aStream nextPutAll:'/'.
+ self d1 printOn:aStream. aStream nextPutAll:'/'.
+ self d2 printOn:aStream. aStream nextPutAll:'/'.
+ self d3 printOn:aStream.
+ ^ self.
+
PrintfScanf printf:'%g' on:aStream argument:self.
"/ self
@@ -2945,16 +4447,16 @@
"/ fixed:true showPositive:false uppercase:false fillChar:(Character space)
"
- (1.2345 asQDouble) printString
- (2 asQDouble squared) printString
-
- (1.2345 asQDouble) printString.
- (1.2345 asFloat) printString.
- (1.2345 asLongFloat) printString.
+ (1.2345 asQDouble) printString
+ (2 asQDouble squared) printString
+
+ (1.2345 asQDouble) printString.
+ (1.2345 asFloat) printString.
+ (1.2345 asLongFloat) printString.
(1.2345 asShortFloat) printString.
- ((QDouble fromFloat:1.2345) / 10.0) printString
- ((QDouble fromFloat:1.2345) / 10000.0) printString
+ ((QDouble fromFloat:1.2345) / 10.0) printString
+ ((QDouble fromFloat:1.2345) / 10000.0) printString
((QDouble fromFloat:1.2345) / 1000000000.0) printString -> '0.0000123449999999999987156270014193593714e-4'
(1.2345 / 1000000000.0) printString -> '1.2345E-09'
"
@@ -3222,9 +4724,9 @@
%{
/* Computes the nearest integer to d. */
-#define nint(d) (((d) == floor(d)) ? (d) : floor((d) + 0.5))
-
- double *a = __QuadDoubleInstPtr(self)->d_quadDoubleValue;
+#define nint(d) (((d) == floor(d)) ? (d) : floor((d) + 0.5))
+
+ double *a = __QDoubleInstPtr(self)->d_qDoubleValue;
OBJ newQD;
double x0, x1, x2, x3;
@@ -3232,31 +4734,31 @@
x1 = x2 = x3 = 0.0;
if (x0 == a[0]) {
- /* First double is already an integer. */
- x1 = nint(a[1]);
-
- if (x1 == a[1]) {
- /* Second double is already an integer. */
- x2 = nint(a[2]);
-
- if (x2 == a[2]) {
- /* Third double is already an integer. */
- x3 = nint(a[3]);
- } else {
- if (abs(x2 - a[2]) == 0.5 && a[3] < 0.0) {
- x2 -= 1.0;
- }
- }
- } else {
- if (abs(x1 - a[1]) == 0.5 && a[2] < 0.0) {
- x1 -= 1.0;
- }
- }
+ /* First double is already an integer. */
+ x1 = nint(a[1]);
+
+ if (x1 == a[1]) {
+ /* Second double is already an integer. */
+ x2 = nint(a[2]);
+
+ if (x2 == a[2]) {
+ /* Third double is already an integer. */
+ x3 = nint(a[3]);
+ } else {
+ if (abs(x2 - a[2]) == 0.5 && a[3] < 0.0) {
+ x2 -= 1.0;
+ }
+ }
+ } else {
+ if (abs(x1 - a[1]) == 0.5 && a[2] < 0.0) {
+ x1 -= 1.0;
+ }
+ }
} else {
- /* First double is not an integer. */
- if (abs(x0 - a[0]) == 0.5 && a[1] < 0.0) {
- x0 -= 1.0;
- }
+ /* First double is not an integer. */
+ if (abs(x0 - a[0]) == 0.5 && a[1] < 0.0) {
+ x0 -= 1.0;
+ }
}
m_renorm4(x0, x1, x2, x3);
__qNew_qdReal(newQD, x0, x1, x2, x3);
@@ -3264,28 +4766,28 @@
%}.
"
- (QDouble fromFloat:4.0) roundedAsFloat
- (QDouble fromFloat:4.6) roundedAsFloat
- (QDouble fromFloat:4.50000001) roundedAsFloat
- (QDouble fromFloat:4.5) roundedAsFloat
- (QDouble fromFloat:4.49999999) roundedAsFloat
- (QDouble fromFloat:4.4) roundedAsFloat
- (QDouble fromFloat:4.1) roundedAsFloat
- (QDouble fromFloat:0.1) roundedAsFloat
- (QDouble fromFloat:0.5) roundedAsFloat
- (QDouble fromFloat:0.49999) roundedAsFloat
- (QDouble fromFloat:0.4) roundedAsFloat
-
- (QDouble fromFloat:-4.0) roundedAsFloat
- (QDouble fromFloat:-4.6) roundedAsFloat
- (QDouble fromFloat:-4.4) roundedAsFloat
- (QDouble fromFloat:-4.499999999) roundedAsFloat
- (QDouble fromFloat:-4.5) roundedAsFloat
- (QDouble fromFloat:-4.5000000001) roundedAsFloat
- (QDouble fromFloat:-4.1) roundedAsFloat
- (QDouble fromFloat:-0.1) roundedAsFloat
- (QDouble fromFloat:-0.5) roundedAsFloat
- (QDouble fromFloat:-0.4) roundedAsFloat
+ (QDouble fromFloat:4.0) roundedAsFloat
+ (QDouble fromFloat:4.6) roundedAsFloat
+ (QDouble fromFloat:4.50000001) roundedAsFloat
+ (QDouble fromFloat:4.5) roundedAsFloat
+ (QDouble fromFloat:4.49999999) roundedAsFloat
+ (QDouble fromFloat:4.4) roundedAsFloat
+ (QDouble fromFloat:4.1) roundedAsFloat
+ (QDouble fromFloat:0.1) roundedAsFloat
+ (QDouble fromFloat:0.5) roundedAsFloat
+ (QDouble fromFloat:0.49999) roundedAsFloat
+ (QDouble fromFloat:0.4) roundedAsFloat
+
+ (QDouble fromFloat:-4.0) roundedAsFloat
+ (QDouble fromFloat:-4.6) roundedAsFloat
+ (QDouble fromFloat:-4.4) roundedAsFloat
+ (QDouble fromFloat:-4.499999999) roundedAsFloat
+ (QDouble fromFloat:-4.5) roundedAsFloat
+ (QDouble fromFloat:-4.5000000001) roundedAsFloat
+ (QDouble fromFloat:-4.1) roundedAsFloat
+ (QDouble fromFloat:-0.1) roundedAsFloat
+ (QDouble fromFloat:-0.5) roundedAsFloat
+ (QDouble fromFloat:-0.4) roundedAsFloat
"
"Created: / 13-06-2017 / 01:52:44 / cg"
@@ -3295,7 +4797,7 @@
renorm
"destructive renormalization"
%{
- double *a = __QuadDoubleInstPtr(self)->d_quadDoubleValue;
+ double *a = __QDoubleInstPtr(self)->d_qDoubleValue;
double c0, c1, c2, c3;
c0 = a[0];
@@ -3324,7 +4826,7 @@
d0
"the most significant (and highest valued) 53 bits of precision"
%{
- RETURN ( __MKFLOAT(__QuadDoubleInstPtr(self)->d_quadDoubleValue[0]) );
+ RETURN ( __MKFLOAT(__QDoubleInstPtr(self)->d_qDoubleValue[0]) );
%}
"Created: / 12-06-2017 / 20:15:12 / cg"
@@ -3334,7 +4836,7 @@
d1
"the next most significant (and next highest valued) 53 bits of precision"
%{
- RETURN ( __MKFLOAT(__QuadDoubleInstPtr(self)->d_quadDoubleValue[1]) );
+ RETURN ( __MKFLOAT(__QDoubleInstPtr(self)->d_qDoubleValue[1]) );
%}
"Created: / 12-06-2017 / 20:15:12 / cg"
@@ -3343,7 +4845,7 @@
d2
%{
- RETURN ( __MKFLOAT(__QuadDoubleInstPtr(self)->d_quadDoubleValue[2]) );
+ RETURN ( __MKFLOAT(__QDoubleInstPtr(self)->d_qDoubleValue[2]) );
%}
"Created: / 12-06-2017 / 20:15:29 / cg"
@@ -3352,7 +4854,7 @@
d3
"the least significant (and smallest valued) 53 bits of precision"
%{
- RETURN ( __MKFLOAT(__QuadDoubleInstPtr(self)->d_quadDoubleValue[3]) );
+ RETURN ( __MKFLOAT(__QDoubleInstPtr(self)->d_qDoubleValue[3]) );
%}
"Created: / 12-06-2017 / 20:15:32 / cg"
@@ -3421,9 +4923,9 @@
"
(QDouble fromFloat:1.0) positive
(QDouble fromFloat:-1.0) positive
- (1.0 asQDouble + 1e-100) positive
- (0.0 asQDouble + 1e-100) positive
- (0.0 asQDouble - 1e-100) positive
+ (1.0 asQDouble + 1e-100) positive
+ (0.0 asQDouble + 1e-100) positive
+ (0.0 asQDouble - 1e-100) positive
"
"Created: / 13-06-2017 / 01:56:53 / cg"
@@ -3437,23 +4939,23 @@
^ self d0 sign
"
- Float nan isNaN
- Float nan sign
- Float infinity sign
- Float infinity negated sign
-
- ShortFloat nan isNaN
- ShortFloat nan sign
- ShortFloat infinity sign
- ShortFloat infinity negated sign
-
- QDouble nan isNaN
- QDouble nan sign
- QDouble infinity sign
- QDouble infinity negated sign
- 0 asQDouble sign
- 1 asQDouble sign
- -1 asQDouble sign
+ Float nan isNaN
+ Float nan sign
+ Float infinity sign
+ Float infinity negated sign
+
+ ShortFloat nan isNaN
+ ShortFloat nan sign
+ ShortFloat infinity sign
+ ShortFloat infinity negated sign
+
+ QDouble nan isNaN
+ QDouble nan sign
+ QDouble infinity sign
+ QDouble infinity negated sign
+ 0 asQDouble sign
+ 1 asQDouble sign
+ -1 asQDouble sign
"
! !
@@ -3468,12 +4970,12 @@
^ f d0 asInteger + f d1 asInteger + f d2 asInteger + f d3 asInteger
"
- (QDouble fromFloat:4.0) ceiling
- (QDouble fromFloat:4.1) ceiling
- (QDouble fromFloat:0.1) ceiling
- (0.1 + (QDouble fromFloat:1.0)) ceiling
+ (QDouble fromFloat:4.0) ceiling
+ (QDouble fromFloat:4.1) ceiling
+ (QDouble fromFloat:0.1) ceiling
+ (0.1 + (QDouble fromFloat:1.0)) ceiling
(1e20 + (QDouble fromFloat:1.0)) ceiling
- (1e20 + (QDouble fromFloat:1.1)) ceiling
+ (1e20 + (QDouble fromFloat:1.1)) ceiling
(QDouble fromFloat:1.5) ceiling
(QDouble fromFloat:0.5) ceiling
@@ -3489,8 +4991,9 @@
It may be useful, if the result is to be further used in another float-operation."
%{
- double *a = __QuadDoubleInstPtr(self)->d_quadDoubleValue;
+ double *a = __QDoubleInstPtr(self)->d_qDoubleValue;
OBJ newQD;
+ int savedCV;
double x0, x1, x2, x3;
x1 = x2 = x3 = 0.0;
@@ -3504,8 +5007,9 @@
x3 = ceil(a[3]);
}
}
-
+ fpu_fix_start(&savedCV);
m_renorm4(x0, x1, x2, x3);
+ fpu_fix_end(&savedCV);
}
__qNew_qdReal(newQD, x0, x1, x2, x3);
RETURN( newQD );
@@ -3554,8 +5058,9 @@
"return the receiver truncated towards negative infinity"
%{
- double *a = __QuadDoubleInstPtr(self)->d_quadDoubleValue;
+ double *a = __QDoubleInstPtr(self)->d_qDoubleValue;
OBJ newQD;
+ int savedCV;
double x0, x1, x2, x3;
x1 = x2 = x3 = 0.0;
@@ -3569,8 +5074,9 @@
x3 = floor(a[3]);
}
}
-
+ fpu_fix_start(&savedCV);
m_renorm4(x0, x1, x2, x3);
+ fpu_fix_end(&savedCV);
}
__qNew_qdReal(newQD, x0, x1, x2, x3);
RETURN( newQD );
@@ -3603,28 +5109,28 @@
^ f d0 asInteger + f d1 asInteger + f d2 asInteger + f d3 asInteger
"
- (QDouble fromFloat:4.0) rounded
- (QDouble fromFloat:4.6) rounded
- (QDouble fromFloat:4.50000001) rounded
- (QDouble fromFloat:4.5) rounded
- (QDouble fromFloat:4.49999999) rounded
- (QDouble fromFloat:4.4) rounded
- (QDouble fromFloat:4.1) rounded
- (QDouble fromFloat:0.1) rounded
- (QDouble fromFloat:0.5) rounded
- (QDouble fromFloat:0.49999) rounded
- (QDouble fromFloat:0.4) rounded
-
- (QDouble fromFloat:-4.0) rounded
- (QDouble fromFloat:-4.6) rounded
- (QDouble fromFloat:-4.4) rounded
- (QDouble fromFloat:-4.499999999) rounded
- (QDouble fromFloat:-4.5) rounded
- (QDouble fromFloat:-4.5000000001) rounded
- (QDouble fromFloat:-4.1) rounded
- (QDouble fromFloat:-0.1) rounded
- (QDouble fromFloat:-0.5) rounded
- (QDouble fromFloat:-0.4) rounded
+ (QDouble fromFloat:4.0) rounded
+ (QDouble fromFloat:4.6) rounded
+ (QDouble fromFloat:4.50000001) rounded
+ (QDouble fromFloat:4.5) rounded
+ (QDouble fromFloat:4.49999999) rounded
+ (QDouble fromFloat:4.4) rounded
+ (QDouble fromFloat:4.1) rounded
+ (QDouble fromFloat:0.1) rounded
+ (QDouble fromFloat:0.5) rounded
+ (QDouble fromFloat:0.49999) rounded
+ (QDouble fromFloat:0.4) rounded
+
+ (QDouble fromFloat:-4.0) rounded
+ (QDouble fromFloat:-4.6) rounded
+ (QDouble fromFloat:-4.4) rounded
+ (QDouble fromFloat:-4.499999999) rounded
+ (QDouble fromFloat:-4.5) rounded
+ (QDouble fromFloat:-4.5000000001) rounded
+ (QDouble fromFloat:-4.1) rounded
+ (QDouble fromFloat:-0.1) rounded
+ (QDouble fromFloat:-0.5) rounded
+ (QDouble fromFloat:-0.4) rounded
"
!
@@ -3632,33 +5138,33 @@
"return the receiver truncated towards negative infinity"
self positive ifTrue:[
- ^ self nintAsFloat
+ ^ self nintAsFloat
].
^ self negated nintAsFloat negated
"
- (QDouble fromFloat:4.0) roundedAsFloat
- (QDouble fromFloat:4.6) roundedAsFloat
- (QDouble fromFloat:4.50000001) roundedAsFloat
- (QDouble fromFloat:4.5) roundedAsFloat
- (QDouble fromFloat:4.49999999) roundedAsFloat
- (QDouble fromFloat:4.4) roundedAsFloat
- (QDouble fromFloat:4.1) roundedAsFloat
- (QDouble fromFloat:0.1) roundedAsFloat
- (QDouble fromFloat:0.5) roundedAsFloat
- (QDouble fromFloat:0.49999) roundedAsFloat
- (QDouble fromFloat:0.4) roundedAsFloat
-
- (QDouble fromFloat:-4.0) roundedAsFloat
- (QDouble fromFloat:-4.6) roundedAsFloat
- (QDouble fromFloat:-4.4) roundedAsFloat
- (QDouble fromFloat:-4.499999999) roundedAsFloat
- (QDouble fromFloat:-4.5) roundedAsFloat
- (QDouble fromFloat:-4.5000000001) roundedAsFloat
- (QDouble fromFloat:-4.1) roundedAsFloat
- (QDouble fromFloat:-0.1) roundedAsFloat
- (QDouble fromFloat:-0.5) roundedAsFloat
- (QDouble fromFloat:-0.4) roundedAsFloat
+ (QDouble fromFloat:4.0) roundedAsFloat
+ (QDouble fromFloat:4.6) roundedAsFloat
+ (QDouble fromFloat:4.50000001) roundedAsFloat
+ (QDouble fromFloat:4.5) roundedAsFloat
+ (QDouble fromFloat:4.49999999) roundedAsFloat
+ (QDouble fromFloat:4.4) roundedAsFloat
+ (QDouble fromFloat:4.1) roundedAsFloat
+ (QDouble fromFloat:0.1) roundedAsFloat
+ (QDouble fromFloat:0.5) roundedAsFloat
+ (QDouble fromFloat:0.49999) roundedAsFloat
+ (QDouble fromFloat:0.4) roundedAsFloat
+
+ (QDouble fromFloat:-4.0) roundedAsFloat
+ (QDouble fromFloat:-4.6) roundedAsFloat
+ (QDouble fromFloat:-4.4) roundedAsFloat
+ (QDouble fromFloat:-4.499999999) roundedAsFloat
+ (QDouble fromFloat:-4.5) roundedAsFloat
+ (QDouble fromFloat:-4.5000000001) roundedAsFloat
+ (QDouble fromFloat:-4.1) roundedAsFloat
+ (QDouble fromFloat:-0.1) roundedAsFloat
+ (QDouble fromFloat:-0.5) roundedAsFloat
+ (QDouble fromFloat:-0.4) roundedAsFloat
"
"Created: / 13-06-2017 / 01:52:44 / cg"