"{ Encoding: utf8 }"
"
COPYRIGHT (c) 2017 by eXept Software AG
All Rights Reserved
This software is furnished under a license and may be used
only in accordance with the terms of that license and with the
inclusion of the above copyright notice. This software may not
be provided or otherwise made available to, or used by, any
other person. No title to or ownership of the software is
hereby transferred.
"
"{ Package: 'stx:libbasic2' }"
"{ NameSpace: Smalltalk }"
LimitedPrecisionReal variableByteSubclass:#QDouble
instanceVariableNames:''
classVariableNames:'DefaultPrintFormat E Epsilon FMax FMin InvFact Ln10 Ln2 NaN Pi
QDoubleOne QDoubleZero'
poolDictionaries:''
category:'Magnitude-Numbers'
!
!QDouble primitiveDefinitions!
%{
#include <stdio.h>
#include <errno.h>
#define __USE_ISOC9X 1
#define __USE_ISOC99 1
#include <math.h>
/*
* on some systems, errno is a macro ... check for it here
*/
#ifndef errno
extern errno;
#endif
#if !defined (__win32__)
# include <locale.h>
#endif
#if defined (__aix__)
# include <float.h>
#endif
#if defined(__irix__) || defined(__solaris__) || defined(__sunos__)
# include <nan.h>
#endif
#if defined(__linux__)
# ifndef NAN
# include <bits/nan.h>
# endif
#endif
#ifdef __win32__
# ifndef isinf
# define isinf(x) \
((((unsigned int *)(&x))[0] == 0x00000000) && \
((((unsigned int *)(&x))[1] & 0x7FF00000) == 0x7FF00000))
# endif
#endif
#if defined(__x86__) || defined(__x86_64__)
# ifndef _FPU_EXTENDED
# define _FPU_EXTENDED 0x0300
# endif
# ifndef _FPU_DOUBLE
# define _FPU_DOUBLE 0x0200
# endif
# if defined( __win32__ ) && (defined( __BORLANDC__ ) || defined( __VISUALC__ ))
# define fpu_fix_start(old_cw_ptr)\
{\
*old_cw_ptr = _control87(0, 0); \
_control87(_FPU_DOUBLE, _FPU_EXTENDED);\
}
# define fpu_fix_end(old_cw_ptr)\
{\
_control87(*old_cw_ptr, _FPU_EXTENDED);\
}
# else // assume MINGW, GCC or CLANG
# ifndef _FPU_GETCW
# define _FPU_GETCW(x) asm volatile ("fnstcw %0":"=m" (x));
# endif
# ifndef _FPU_SETCW
# define _FPU_SETCW(x) asm volatile ("fldcw %0": :"m" (x));
# endif
# define fpu_fix_start(old_cw_ptr)\
{\
volatile unsigned short cw, new_cw;\
_FPU_GETCW(cw);\
new_cw = (cw & ~_FPU_EXTENDED) | _FPU_DOUBLE;\
_FPU_SETCW(new_cw);\
*old_cw_ptr = cw;\
}
# define fpu_fix_end(old_cw_ptr)\
{\
volatile unsigned short cw = *old_cw_ptr;\
_FPU_SETCW(cw);\
}
# endif
#endif
#define __qNew_qdReal(newQD, d0,d1,d2,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)
#define _QD_SPLITTER 134217729.0 // = 2^27 + 1
#define _QD_SPLIT_THRESH 6.69692879491417e+299 // = 2^996
#define m_quick_two_sum(rslt_1, a_1, b_1, err_1)\
{\
double s_1 = (a_1) + (b_1);\
(err_1) = (b_1) - (s_1 - (a_1));\
(rslt_1) = s_1; \
}
#define m_quick_two_diff(rslt_2, a_2, b_2, err_2)\
{\
double s_2 = (a_2) - (b_2);\
(err_2) = ((a_2) - s_2) - (b_2);\
(rslt_2) = s_2;\
}
#define m_two_sum(rslt_3, a_3, b_3, err_3)\
{\
double s_3 = (a_3) + (b_3);\
double v_3 = s_3 - (a_3);\
(err_3) = ((a_3) - (s_3 - v_3)) + ((b_3) - v_3);\
(rslt_3) = s_3;\
}
/* Computes fl(a-b) and err(a-b). */
#define m_two_diff(rslt_4, a_4, b_4, err_4)\
{\
double s_4 = (a_4) - (b_4);\
double bb_4 = s_4 - (a_4);\
(err_4) = ((a_4) - (s_4 - bb_4)) - ((b_4) + bb_4);\
(rslt_4) = s_4;\
}
#define m_three_sum(a_5, b_5, c_5)\
{ \
double t1_5, t2_5, t3_5; \
m_two_sum(t1_5, (a_5), (b_5), t2_5); \
m_two_sum((a_5), (c_5), t1_5, t3_5); \
m_two_sum((b_5), t2_5, t3_5, (c_5)); \
}
#define m_three_sum2(a_6, b_6, c_6)\
{\
double t1_6, t2_6, t3_6;\
m_two_sum(t1_6, (a_6), (b_6), t2_6);\
m_two_sum((a_6), (c_6), t1_6, t3_6);\
(b_6) = t2_6 + t3_6;\
}
#ifndef QD_FMS
/* Computes high word and lo word of a */
#define m_split(a_7, hi_7, lo_7)\
{\
double temp_7;\
double thi_7, tlo_7;\
if ((a_7) > _QD_SPLIT_THRESH || (a_7) < -_QD_SPLIT_THRESH) {\
(a_7) *= 3.7252902984619140625e-09;\
temp_7 = _QD_SPLITTER * (a_7);\
thi_7 = temp_7 - (temp_7 - (a_7));\
tlo_7 = (a_7) - thi_7;\
thi_7 *= 268435456.0;\
tlo_7 *= 268435456.0;\
} else {\
temp_7 = _QD_SPLITTER * (a_7);\
thi_7 = temp_7 - (temp_7 - (a_7));\
tlo_7 = (a_7) - thi_7;\
}\
(hi_7) = thi_7; \
(lo_7) = tlo_7; \
}
#endif
#ifdef QD_FMS
/* Computes fl(a*b) and err(a*b). */
#define m_two_prod(rslt_8, a_8, b_8, err_8)\
{\
double p_8 = (a_8) * (b_8);\
err_8 = QD_FMS((a_8), (b_8), p_8);\
rslt_8 = p_8; \
}
#else
#define m_two_prod(rslt_8, a_8, b_8, err_8)\
{\
double a_hi_8, a_lo_8, b_hi_8, b_lo_8;\
double p_8 = (a_8) * (b_8);\
m_split(a_8, a_hi_8, a_lo_8);\
m_split(b_8, b_hi_8, b_lo_8);\
err_8 = ((a_hi_8 * b_hi_8 - p_8) + a_hi_8 * b_lo_8 + a_lo_8 * b_hi_8) + a_lo_8 * b_lo_8;\
rslt_8 = p_8; \
}
#endif
#ifdef QD_FMS
#define m_two_sqr(rslt_9, a_9, err_9)\
{\
double p_9 = (a_9) * (a_9);\
err_9 = QD_FMS((a_9), (a_9), p_9);\
rslt_9 = p_9;\
}
#else
#define m_two_sqr(rslt_9, a_9, err_9)\
{\
double hi_9, lo_9;\
double q_9 = (a_9) * (a_9);\
m_split(a_9, hi_9, lo_9);\
err_9 = ((hi_9 * hi_9 - q_9) + 2.0 * hi_9 * lo_9) + lo_9 * lo_9;\
rslt_9 = q_9;\
}
#endif
#define m_renorm4(c0_10, c1_10, c2_10, c3_10)\
{\
double s0_10, s1_10, s2_10 = 0.0, s3_10 = 0.0;\
\
if (! isinf(c0_10)) { \
\
m_quick_two_sum(s0_10, c2_10, c3_10, c3_10);\
m_quick_two_sum(s0_10, c1_10, s0_10, c2_10);\
m_quick_two_sum(c0_10, c0_10, s0_10, c1_10);\
\
s0_10 = c0_10;\
s1_10 = c1_10;\
if (s1_10 != 0.0) {\
m_quick_two_sum(s1_10, s1_10, c2_10, s2_10);\
if (s2_10 != 0.0) {\
m_quick_two_sum(s2_10, s2_10, c3_10, s3_10);\
} else {\
m_quick_two_sum(s1_10, s1_10, c3_10, s2_10);\
}\
} else {\
m_quick_two_sum(s0_10, s0_10, c2_10, s1_10);\
if (s1_10 != 0.0) {\
m_quick_two_sum(s1_10, s1_10, c3_10, s2_10);\
} else {\
m_quick_two_sum(s0_10, s0_10, c3_10, s1_10);\
}\
}\
\
c0_10 = s0_10;\
c1_10 = s1_10;\
c2_10 = s2_10;\
c3_10 = s3_10;\
}\
}
#define m_renorm5(c0_11, c1_11, c2_11, c3_11, c4_11)\
{\
double s0_11, s1_11, s2_11 = 0.0, s3_11 = 0.0; \
\
if (! isinf(c0_11)) { \
m_quick_two_sum(s0_11, c3_11, c4_11, c4_11); \
m_quick_two_sum(s0_11, c2_11, s0_11, c3_11); \
m_quick_two_sum(s0_11, c1_11, s0_11, c2_11); \
m_quick_two_sum(c0_11, c0_11, s0_11, c1_11); \
\
s0_11 = c0_11; \
s1_11 = c1_11; \
\
m_quick_two_sum(s0_11, c0_11, c1_11, s1_11); \
if (s1_11 != 0.0) { \
m_quick_two_sum(s1_11, s1_11, c2_11, s2_11); \
if (s2_11 != 0.0) { \
m_quick_two_sum(s2_11 ,s2_11, c3_11, s3_11); \
if (s3_11 != 0.0) {\
s3_11 += c4_11; \
} else {\
s2_11 += c4_11;\
}\
} else { \
m_quick_two_sum(s1_11, s1_11, c3_11, s2_11); \
if (s2_11 != 0.0) {\
m_quick_two_sum(s2_11, s2_11, c4_11, s3_11); \
} else { \
m_quick_two_sum(s1_11, s1_11, c4_11, s2_11); \
} \
} \
} else { \
m_quick_two_sum(s0_11,s0_11, c2_11, s1_11); \
if (s1_11 != 0.0) { \
m_quick_two_sum(s1_11,s1_11, c3_11, s2_11); \
if (s2_11 != 0.0) {\
m_quick_two_sum(s2_11,s2_11, c4_11, s3_11); \
} else { \
m_quick_two_sum(s1_11 ,s1_11, c4_11, s2_11); \
} \
} else { \
m_quick_two_sum(s0_11,s0_11, c3_11, s1_11); \
if (s1_11 != 0.0) { \
m_quick_two_sum(s1_11,s1_11, c4_11, s2_11); \
} else { \
m_quick_two_sum(s0_11,s0_11, c4_11, s1_11); \
} \
} \
} \
\
c0_11 = s0_11; \
c1_11 = s1_11; \
c2_11 = s2_11; \
c3_11 = s3_11; \
} \
}
%}
! !
!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)
// renorm : 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
// qd_sqrt : square root (scalar)
// qd_pow : qd ^ n (n integer)
static INLINE void
fast_two_sum(double *s, double *e, double a, double b)
{
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 *o0, double *o1, double *o2, double *o3, double a0, double a1, double a2, double a3, double b)
{
double e,x,y,w,z;
double c0, c1, c2, c3;
c0 = 0.0; c1 = 0.0; c2 = 0.0; c3 = 0.0;
two_sum(&x, &y, a0, b); c0 = x; e = y;
two_sum(&x, &y, a1, e); c1 = x; e = y;
two_sum(&x, &y, a2, e); c2 = x; e = y;
two_sum(&x, &y, a3, e); c3 = x; e = y;
renorm(&x, &y, &w, &z, c0, c1, c2, c3, e);
o0[0] = x; o1[0] = y; o2[0] = w; o3[0] = z;
}
//--------------------------------------------------------------------------------------------
//
// addition quad-double + double-double
//
//--------------------------------------------------------------------------------------------
static INLINE void
qd_add_dd(double *o0, double *o1, double *o2, double *o3, double a0, double a1, double a2, double a3, double b0, double b1)
{
double e1,e2,x,y,w,z;
double c0, c1, c2, c3;
c0 = 0.0; c1 = 0.0; c2 = 0.0; c3 = 0.0;
two_sum(&x, &y, a0, b0); c0 = x; e1 = y;
two_sum(&x, &y, a1, b1); c1 = x; e2 = y;
two_sum(&x, &y, c1, e1); c1 = x; e1 = y;
two_sum(&x, &y, a2, e2); c2 = x; e2 = y;
two_sum(&x, &y, c2, e1); c2 = x; e1 = y;
two_sum(&x, &y, e1, e2); e1 = x; e2 = y;
two_sum(&x, &y, a3, e1); c3 = x; e1 = y;
e1 = e1 + e2;
renorm(&x, &y, &w, &z, c0, c1, c2, c3, e1);
o0[0] = x; o1[0] = y; o2[0] = w; o3[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 *o0, double *o1, double *o2, double *o3, double a, double b0, double b1, double b2, double b3)
{
double e,x,y,w,z;
double c0, c1, c2, c3;
e=0.0;
c0 = 0.0; c1 = 0.0; c2 = 0.0; c3 = 0.0;
b0=-b0; b1=-b1; b2=-b2; b3=-b3;
two_sum(&x, &y, a, b0);
c0 = x;
e = y;
two_sum(&x, &y, b1, e);
c1 = x;
e = y;
two_sum(&x, &y, b2, e);
c2 = x;
e = y;
two_sum(&x, &y, b3, e);
c3 = x;
e = y;
renorm(&x, &y, &w, &z, c0, c1, c2, c3, e);
o0[0] = x;
o1[0] = y;
o2[0] = w;
o3[0] = z;
}
//--------------------------------------------------------------------------------------------
//
// 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 *o0, double *o1, double *o2, double *o3, double b, double a0, double a1, double a2, double a3)
{
double e0,e1,e2,x,y,w,z;
double c0, c1, c2, c3;
c0 = 0.0; c1 = 0.0; c2 = 0.0; c3 = 0.0;
two_prod(&x, &y, a0, b); c0 = x; e0 = y;
two_prod(&x, &y, a1, b); c1 = x; e1 = y;
two_sum(&x, &y, c1, e0); c1 = x; e0 = y;
two_prod(&x, &y, a2, b); c2 = x; e2 = y;
three_sum(&x, &y, &z, c2, e1, e0); c2 = x; e0 = y; e1 = z;
c3 = a3*b;
three_sum2(&x, &y, c3, e2, e0); c3 = x; e0 = y;
e0 = e0 + e1;
renorm(&x, &y, &w, &z, c0, c1, c2, c3, e0);
o0[0] = x; o1[0] = y; o2[0] = w; o3[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
qd_sqrt(double *c0, double *c1, double *c2, double *c3, double a0, double a1, double a2, double a3)
{
double h0,h1,h2,h3,x0,x1,x2,x3,p,q,r,s;
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
qd_pow(double *c0, double *c1, double *c2, double *c3, double a0, double a1, double a2, double a3, int p)
{
double r0,r1,r2,r3,x,y,w,z;
int abs_p;
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 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
sin_taylor_qd(double *s0, double *s1, double *s2, double *s3, double x0, double x1, double x2, double x3)
{
double eps = 1.21543267145725e-63; // = 2^-209
double thresh = 0.5*fabs(x0)*eps;
double y0,y1,y2,y3,r0,r1,r2,r3,t0,t1,t2,t3;
int i;
if(x0==0.0) {
s0[0]=0.0; s1[0]=0.0; s2[0]=0.0; s3[0]=0.0;
return;
}
i=0;
qd_mul_qd(&y0,&y1,&y2,&y3,x0,x1,x2,x3,x0,x1,x2,x3);
y0 = -y0; y1 = -y1; y2 = -y2; y3 = -y3;
s0[0]=x0; s1[0]=x1; s2[0]=x2; s3[0]=x3;
r0=x0; r1=x1; r2=x2; r3=x3;
qd_mul_qd(&r0,&r1,&r2,&r3,r0,r1,r2,r3,y0,y1,y2,y3);
qd_mul_qd(&t0,&t1,&t2,&t3,r0,r1,r2,r3,inv_fact[i][0],inv_fact[i][1],inv_fact[i][2],inv_fact[i][3]);
qd_add_qd(&s0[0],&s1[0],&s2[0],&s3[0],s0[0],s1[0],s2[0],s3[0],t0,t1,t2,t3);
i=i+2;
while ((i<=15)||(fabs(t0)>thresh)) {
qd_mul_qd(&r0,&r1,&r2,&r3,r0,r1,r2,r3,y0,y1,y2,y3);
qd_mul_qd(&t0,&t1,&t2,&t3,r0,r1,r2,r3,inv_fact[i][0],inv_fact[i][1],inv_fact[i][2],inv_fact[i][3]);
qd_add_qd(&s0[0],&s1[0],&s2[0],&s3[0],s0[0],s1[0],s2[0],s3[0],t0,t1,t2,t3);
i=i+2;
}
}
static void
cos_taylor_qd(double *c0, double *c1, double *c2, double *c3, double x0, double x1, double x2, double x3)
{
double eps = 1.21543267145725e-63; // = 2^-209
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);
qd_sqrt(&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
qd_sin(double *s0, double *s1, double *s2, double *s3, double *a0, double *a1, double *a2, double *a3)
{
double p0,p1,p2,p3, q0,q1,q2,q3, z0,z1,z2,z3, r0,r1,r2,r3, j, t0,t1,t2,t3, k,abs_k;
double _2pi[4] = {6.283185307179586232e+00,2.449293598294706414e-16,-5.989539619436679332e-33,2.224908441726730563e-49}; // 2*pi
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
qd_cos(double *c0, double *c1, double *c2, double *c3, double *a0, double *a1, double *a2, double *a3)
{
double p0,p1,p2,p3, q0,q1,q2,q3, z0,z1,z2,z3, r0,r1,r2,r3, j, t0,t1,t2,t3, k,abs_k;
double _2pi[4] = {6.283185307179586232e+00,2.449293598294706414e-16,-5.989539619436679332e-33,2.224908441726730563e-49}; // 2*pi
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
qd_tan(double *t0, double *t1, double *t2, double *t3, double *a0, double *a1, double *a2, double *a3)
{
double sin0,sin1,sin2,sin3,cos0,cos1,cos2,cos3;
sincos_qd(&sin0,&sin1,&sin2,&sin3,&cos0,&cos1,&cos2,&cos3,a0[0],a1[0],a2[0],a3[0]);
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
qd_exp(double *e0, double *e1, double *e2, double *e3, double *x0, double *x1, double *x2, double *x3) {
double k = ldexp(1.0,16);
double inv_k = 1.0 / k;
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
quick_two_sum(double a, double b, double *errPtr) {
double s = a + b;
*errPtr = b - (s - a);
return s;
}
/* Computes fl(a-b) and err(a-b). Assumes |a| >= |b| */
INLINE double
quick_two_diff(double a, double b, double *errPtr) {
double s = a - b;
*errPtr = (a - s) - b;
return s;
}
/* Computes fl(a+b) and err(a+b). */
INLINE double
two_sum(double a, double b, double *errPtr) {
double s = a + b;
double bb = s - a;
*errPtr = (a - (s - bb)) + (b - bb);
return s;
}
/* Computes fl(a-b) and err(a-b). */
INLINE double
two_diff(double a, double b, double *errPtr) {
double s = a - b;
double bb = s - a;
*errPtr = (a - (s - bb)) - (b + bb);
return s;
}
#ifndef QD_FMS
/* Computes high word and lo word of a */
INLINE void
split(double a, double *hiPtr, double *loPtr) {
double temp;
if (a > _QD_SPLIT_THRESH || a < -_QD_SPLIT_THRESH) {
a *= 3.7252902984619140625e-09; // 2^-28
temp = _QD_SPLITTER * a;
*hiPtr = temp - (temp - a);
*loPtr = a - *hiPtr;
*hiPtr *= 268435456.0; // 2^28
*loPtr *= 268435456.0; // 2^28
} else {
temp = _QD_SPLITTER * a;
*hiPtr = temp - (temp - a);
*loPtr = a - *hiPtr;
}
}
#endif
/* Computes fl(a*b) and err(a*b). */
INLINE double
two_prod(double a, double b, double *errPtr) {
#ifdef QD_FMS
double p = a * b;
*errPtr = QD_FMS(a, b, p);
return p;
#else
double a_hi, a_lo, b_hi, b_lo;
double p = a * b;
split(a, &a_hi, &a_lo);
split(b, &b_hi, &b_lo);
*errPtr = ((a_hi * b_hi - p) + a_hi * b_lo + a_lo * b_hi) + a_lo * b_lo;
return p;
#endif
}
/* Computes fl(a*a) and err(a*a). Faster than the above method. */
INLINE double
two_sqr(double a, double *errPtr) {
#ifdef QD_FMS
double p = a * a;
*errPtr = QD_FMS(a, a, p);
return p;
#else
double hi, lo;
double q = a * a;
split(a, &hi, &lo);
*errPtr = ((hi * hi - q) + 2.0 * hi * lo) + lo * lo;
return q;
#endif
}
/* Computes the nearest integer to d. */
INLINE double
nint(double d) {
if (d == floor(d))
return d;
return floor(d + 0.5);
}
/* Computes the truncated integer. */
INLINE double
aint(double d) {
return (d >= 0.0) ? floor(d) : ceil(d);
}
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;
s0 = quick_two_sum(*c2Ptr, *c3Ptr, c3Ptr);
s0 = quick_two_sum(*c1Ptr, s0, c2Ptr);
c0 = quick_two_sum(c0, s0, c1Ptr);
s0 = c0;
s1 = *c1Ptr;
if (s1 != 0.0) {
s1 = quick_two_sum(s1, *c2Ptr, &s2);
if (s2 != 0.0)
s2 = quick_two_sum(s2, *c3Ptr, &s3);
else
s1 = quick_two_sum(s1, *c3Ptr, &s2);
} else {
s0 = quick_two_sum(s0, *c2Ptr, &s1);
if (s1 != 0.0)
s1 = quick_two_sum(s1, *c3Ptr, &s2);
else
s0 = quick_two_sum(s0, *c3Ptr, &s1);
}
*c0Ptr = s0;
*c1Ptr = s1;
*c2Ptr = s2;
*c3Ptr = s3;
}
INLINE void
renorm5(double *c0Ptr, double *c1Ptr, double *c2Ptr, double *c3Ptr, double *c4Ptr) {
double s0, s1, s2 = 0.0, s3 = 0.0;
if (isinf(*c0Ptr)) return;
s0 = quick_two_sum(*c3Ptr, *c4Ptr, c4Ptr);
s0 = quick_two_sum(*c2Ptr, s0, c3Ptr);
s0 = quick_two_sum(*c1Ptr, s0, c2Ptr);
*c0Ptr = quick_two_sum(*c0Ptr, s0, c1Ptr);
s0 = *c0Ptr;
s1 = *c1Ptr;
s0 = quick_two_sum(*c0Ptr, *c1Ptr, &s1);
if (s1 != 0.0) {
s1 = quick_two_sum(s1, *c2Ptr, &s2);
if (s2 != 0.0) {
s2 =quick_two_sum(s2, *c3Ptr, &s3);
if (s3 != 0.0)
s3 += *c4Ptr;
else
s2 += *c4Ptr;
} else {
s1 = quick_two_sum(s1, *c3Ptr, &s2);
if (s2 != 0.0)
s2 = quick_two_sum(s2, *c4Ptr, &s3);
else
s1 = quick_two_sum(s1, *c4Ptr, &s2);
}
} else {
s0 = quick_two_sum(s0, *c2Ptr, &s1);
if (s1 != 0.0) {
s1 = quick_two_sum(s1, *c3Ptr, &s2);
if (s2 != 0.0)
s2 = quick_two_sum(s2, *c4Ptr, &s3);
else
s1 = quick_two_sum(s1, *c4Ptr, &s2);
} else {
s0 = quick_two_sum(s0, *c3Ptr, &s1);
if (s1 != 0.0)
s1 = quick_two_sum(s1, *c4Ptr, &s2);
else
s0 = quick_two_sum(s0, *c4Ptr, &s1);
}
}
*c0Ptr = s0;
*c1Ptr = s1;
*c2Ptr = s2;
*c3Ptr = s3;
}
INLINE void
three_sum(double *aPtr, double *bPtr, double *cPtr) {
double t1, t2, t3;
t1 = two_sum(*aPtr, *bPtr, &t2);
*aPtr = two_sum(*cPtr, t1, &t3);
*bPtr = two_sum(t2, t3, 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);
*bPtr = t2 + t3;
}
#endif
#if 0
/* These are provided to give consistent
interface for double with double-double and quad-double. */
INLINE void
sincosh(double t, double &sinh_t, double &cosh_t) {
sinh_t = sinh(t);
cosh_t = cosh(t);
}
INLINE double
sqr(double t) {
return t * t;
}
INLINE double
to_double(double a) {
return a;
}
INLINE int
to_int(double a) {
return static_cast<int>(a);
}
#endif
%}
! !
!QDouble class methodsFor:'documentation'!
copyright
"
COPYRIGHT (c) 2017 by eXept Software AG
All Rights Reserved
This software is furnished under a license and may be used
only in accordance with the terms of that license and with the
inclusion of the above copyright notice. This software may not
be provided or otherwise made available to, or used by, any
other person. No title to or ownership of the software is
hereby transferred.
"
!
documentation
"
ATTENTION: ongoing, unfinished work.
No warranty that this works correctly...
QDoubles represent rational numbers with extended, but still limited precision.
In contrast to Floats (which use the C-compiler's native 64bit 'double' format),
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.
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
The number of bits:
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')
[author:]
Claus Gittinger
[see also:]
Number
Float ShortFloat LongFloat
Fraction FixedPoint Integer Complex
FloatArray DoubleArray
"
!
examples
"
Floats, LongFloats suffer from loosing bits:
(Float readFrom:'0.3333333333333333333333333333333333333333333333333333333333')
-(Float readFrom:'0.333333333333333333333333333333333333333333333333333333333')
-> 0.0
(Float readFrom:'0.3333333333333333333333333333333333333333333333333333333333')
= (Float readFrom:'0.333333333333333333333333333333333333333333333333333333333')
-> true
(Float readFrom:'0.33333333333333333333333333333333333333333333333333333333333333333333')
= (Float readFrom:'0.3333333333333333333333333333333333333333333333333333333333333333333')
-> true
1000 0110 1000 0101 1000 0101 1000 0101 1000 0101 1000 0101 1101 0101 0011 1111
(Float readFrom:'0.3333333333333333333333333333333333333333333333333333333333')
= (Float readFrom:'0.3333333333333333333333333333333333333333333333333333333333333333333')
(LongFloat readFrom:'0.3333333333333333333333333333333333333333333333333333333333')
-(LongFloat readFrom:'0.333333333333333333333333333333333333333333333333333333333')
-> 0.0
(LongFloat readFrom:'0.3333333333333333333333333333333333333333333333333333333333')
= (LongFloat readFrom:'0.333333333333333333333333333333333333333333333333333333333')
-> 0.0
(QDouble readFrom:'0.3333333333333333333333333333333333333333333333333333333333')
-(QDouble readFrom:'0.333333333333333333333333333333333333333333333333333333333')
(QDouble readFrom:'0.33333333333333333333333333333333333333333333333333333333333')
-(QDouble readFrom:'0.3333333333333333333333333333333333333333333333333333333333')
(QDouble readFrom:'0.33333333333333333333333333333333333333333333333333333333333333333333333333333333333333333')
-(QDouble readFrom:'0.3333333333333333333333333333333333333333333333333333333333333333333333333333333333333333')
"
! !
!QDouble class methodsFor:'instance creation'!
basicNew
"return a new quad-precision double - here we return 0.0
Notice that numbers are usually NOT created this way ...
It's implemented here to allow things like binary store & load
of floats. (but even this support will go away eventually, it's not
a good idea to store the bits of a float - the reader might have a
totally different representation - so floats should be
binary stored in a device independent format."
%{ /* NOCONTEXT */
#ifdef __SCHTEAM__
ERROR("trying to instantiate a qDouble");
#else
OBJ newQD;
__qNew_qdReal(newQD, 0.0, 0.0, 0.0, 0.0);
RETURN (newQD);
#endif /* not SCHTEAM */
%}
"
self basicNew
"
"Created: / 12-06-2017 / 16:00:38 / cg"
!
d0:d0 d1:d1 d2:d2 d3:d3
"return a new quad-precision double from individual double components"
%{ /* NOCONTEXT */
#ifdef __SCHTEAM__
ERROR("trying to instantiate a qDouble");
#else
OBJ newQD;
if (__isFloatLike(d0)
&& __isFloatLike(d1)
&& __isFloatLike(d2)
&& __isFloatLike(d3)) {
__qNew_qdReal(newQD, __floatVal(d0), __floatVal(d1),
__floatVal(d2), __floatVal(d3));
RETURN (newQD);
}
#endif
%}.
self error:'invalid argument'
"
self d0: 3.141592653589793116e+00
d1: 1.224646799147353207e-16
d2: -2.994769809718339666e-33
d3: 1.112454220863365282e-49
"
"Created: / 12-06-2017 / 20:17:14 / cg"
!
fromDoubleArray:aDoubleArray
"return a new quad-precision double from coercing a double array"
%{ /* NOCONTEXT */
#ifdef __SCHTEAM__
ERROR("trying to instantiate a qDouble");
#else
OBJ newQD;
if (__isDoubleArray(aDoubleArray)) {
double* __d__ = __DoubleArrayInstPtr(aDoubleArray)->d_element;
__qNew_qdReal(newQD, __d__[0], __d__[1], __d__[2], __d__[3]);
RETURN (newQD);
}
#endif
%}.
self error:'invalid argument'
"
self fromDoubleArray(DoubleArray
with: 3.141592653589793116e+00
with: 1.224646799147353207e-16
with: -2.994769809718339666e-33
with: 1.112454220863365282e-49)
"
"Created: / 12-06-2017 / 18:25:32 / cg"
!
fromFloat:aFloat
"return a new quad-precision double from coercing aFloat"
%{ /* NOCONTEXT */
#ifdef __SCHTEAM__
ERROR("trying to instantiate a qDouble");
#else
double dVal;
OBJ newQD;
if (__isFloatLike(aFloat)) {
dVal = __floatVal(aFloat);
} else if (__isShortFloat(aFloat)) {
dVal = __shortFloatVal(aFloat);
} else {
goto badArg;
}
__qNew_qdReal(newQD, dVal, 0.0, 0.0, 0.0);
RETURN (newQD);
badArg: ;
#endif
%}.
self argumentError:'invalid (non-float) argument'
"
self fromFloat:1.0
"
"Created: / 12-06-2017 / 16:06:54 / cg"
!
fromInteger:anInteger
"return a new quad-precision double from coercing anInteger"
%{ /* NOCONTEXT */
#ifdef __SCHTEAM__
ERROR("trying to instantiate a qDouble");
#else
OBJ newQD;
if (__isSmallInteger(anInteger)) {
INT iVal = __smallIntegerVal(anInteger);
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))) {
d[0] = (double)(iVal & ~0xFFFFFFFF);
d[1] = (double)(iVal & 0xFFFFFFFF);
renorm(&(d[0]), &(d[1]), &(d[2]), &(d[3]), d[0], d[1], 0.0, 0.0, 0.0);
// renorm4(&(a[0]), &(a[1]), &(a[2]), &(a[3]));
} else {
d[0] = (double)iVal;
}
RETURN (newQD);
}
#endif
%}.
^ super fromInteger:anInteger
"
self fromInteger:2
self fromInteger:16rFFFFFFFF -- 32bit 4294967295.0
self fromInteger:16rFFFFFFFFFFFF -- 48bit 281474976710655.0
self fromInteger:16rFFFFFFFFFFFFF -- 52bit 4503599627370495.0
self fromInteger:16rFFFFFFFFFFFFFF -- 56bit 72057594037927935.0
self fromInteger:16rFFFFFFFFFFFFFFF -- 60bit 1152921504606846975.0
self fromInteger:16r1FFFFFFFFFFFFFFF -- 61bit 2305843009213693951.0
self fromInteger:16r3FFFFFFFFFFFFFFF -- 62bit 4611686018427387903.0
self fromInteger:16r7FFFFFFFFFFFFFFF -- 63bit 9223372036854775807.0
self fromInteger:16rFFFFFFFFFFFFFFFF -- 64bit 18446744073709551615.0
"
"Created: / 12-06-2017 / 16:10:10 / cg"
"Modified: / 04-07-2017 / 12:51:52 / cg"
! !
!QDouble class methodsFor:'coercing & converting'!
coerce:aNumber
"convert the argument aNumber into an instance of the receiver's class and return it."
^ aNumber asQDouble
"Created: / 12-06-2017 / 17:13:47 / cg"
"Modified: / 12-06-2017 / 21:09:06 / cg"
! !
!QDouble class methodsFor:'constants'!
NaN
"return a QDouble which represents not-a-Number (i.e. an invalid number)"
NaN isNil ifTrue:[
NaN := self d0:(Float NaN) d1:(Float NaN) d2:(Float NaN) d3:(Float NaN)
].
^ NaN
"Created: / 21-06-2017 / 20:44:57 / cg"
!
e
"return the constant e as quad precision double.
(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 e printfPrintString:'%.61f'
-> '2.7182818284590452353602874713526624977572470936999595749669676'
Wolfram says:
2.71828182845904523536028747135266249775724709369995957496696762772407663035354759457138217852516642742746
"
"Created: / 12-06-2017 / 18:29:36 / cg"
!
fmax
"return the constant e as quad precision double.
(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
"
Float fmax
self fmax
"
"Created: / 14-06-2017 / 19:14:18 / cg"
!
fmin
"return the smallest representable instance of this class"
FMin isNil ifTrue:[
FMin := Float fmin asQDouble. "/ 1.6259745436952323e-260 asQDouble
].
^ FMin
"
QDouble fmin
Float fmin
"
"Created: / 14-06-2017 / 19:14:49 / cg"
!
infinity
^ Infinity positive
"Created: / 18-06-2017 / 23:41:07 / cg"
!
ln10
"return the constant e as quad precision double.
(returns approx. 200 bits of precision)"
Ln10 isNil ifTrue:[
Ln10 := self d0: 2.302585092994045901e+00
d1: -2.170756223382249351e-16
d2: -9.984262454465776570e-33
d3: -4.023357454450206379e-49
].
^ Ln10
"
self ln10 printfPrintString:'%.61f'
-> '2.3025850929940456840179914546843642076011014886287729760333279'
Wolfram says:
2.30258509299404568401799145468436420760110148862877297603332790096757260967735248023599720508959829834196778404228...
"
"Created: / 12-06-2017 / 18:32:29 / cg"
!
ln2
"return the constant e as quad precision double.
(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 ln2 printfPrintString:'%.61f'
-> '0.6931471805599452709398341558750792990469129794959648865081141'
Wolfram says:
0.69314718055994530941723212145817656807550013436025525412068000949339362196969471560586332699641868754200148102057...
"
"Created: / 12-06-2017 / 18:31:34 / cg"
!
negativeInfinity
^ Infinity negative
"Created: / 18-06-2017 / 23:40:47 / cg"
!
pi
"return the constant pi as quad precision double.
(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 pi printfPrintString:'%.60f'
'3.141592653589793238462643383279502884197169399375105820974945'
Wolfram says:
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068
(QDouble readFrom:'3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253')
printfPrintString:'%.60f'
"
"Created: / 12-06-2017 / 18:27:13 / cg"
!
unity
"return the neutral element for multiplication (1.0) as QDouble"
QDoubleOne isNil ifTrue:[
QDoubleOne := 1.0 asQDouble.
].
^ QDoubleOne
"
self unity
"
"Created: / 15-06-2017 / 11:45:22 / cg"
!
zero
"return the neutral element for addition (0.0) as QDouble"
QDoubleZero isNil ifTrue:[
QDoubleZero := 0.0 asQDouble
].
^ QDoubleZero
"
self zero
"
"Created: / 15-06-2017 / 11:44:13 / cg"
! !
!QDouble class methodsFor:'queries'!
defaultPrintPrecision
"return the number of decimal digits printed by default"
^ 30
"
ShortFloat defaultPrintPrecision
Float defaultPrintPrecision
LongFloat defaultPrintPrecision
QDouble defaultPrintPrecision
QuadFloat defaultPrintPrecision
OctaFloat defaultPrintPrecision
"
"Created: / 17-06-2017 / 02:58:51 / cg"
"Modified: / 21-06-2017 / 13:39:08 / cg"
!
epsilon
"return the maximum relative spacing of instances of mySelf
(i.e. the value-delta of the least significant bit)
see https://en.wikipedia.org/wiki/Machine_epsilon"
"/ ^ 1.2154326714572500565324311366323150942261000827598106963711353e-63
Epsilon isNil ifTrue:[
Epsilon := self computeEpsilon.
].
^ Epsilon
"
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)
"
"Created: / 12-06-2017 / 18:52:44 / cg"
"Modified: / 22-06-2017 / 15:34:56 / cg"
!
numBitsInExponent
"answer the number of bits in the exponent.
I use regular IEEE doubles to store the value,
thus my exponent bits are the same as double's exponent bits"
^ Float numBitsInExponent
"
1.0 asQDouble numBitsInExponent
"
"Created: / 12-06-2017 / 11:11:04 / cg"
"Modified (comment): / 28-05-2019 / 08:55:04 / Claus Gittinger"
!
numBitsInMantissa
"answer the number of bits in the mantissa.
Here, a fake number is returned"
^ (Float numBitsInMantissa - 1) * 4
"
1.0 asFloat numBitsInMantissa
1.0 asShortFloat numBitsInMantissa
1.0 asLongFloat numBitsInMantissa
1.0 asQDouble numBitsInMantissa
1.0 asQDouble class numBitsInMantissa
Float numBitsInMantissa
ShortFloat numBitsInMantissa
QDouble numBitsInMantissa
"
"Created: / 12-06-2017 / 11:13:44 / cg"
"Modified (comment): / 20-06-2017 / 11:05:26 / cg"
"Modified (comment): / 28-05-2019 / 09:07:07 / Claus Gittinger"
!
precision
"answer the number of bits in the mantissa"
"/ subtract some due to overlap in the component numbers
"/ actual precision seems to be more like:
"/ ^ (Float precision) * 4 - 3 + 1.
"/ but I am a bit conservative here:
^ (Float precision - 2) * 4
"
ShortFloat precision -> 24
Float precision -> 53
LongFloat precision -> 64
QDouble precision -> 204
QuadFloat precision -> 113
OctaFloat precision -> 237
1.0 class numBitsInMantissa
1.0 asShortFloat class numBitsInMantissa
1.0 asLongFloat class numBitsInMantissa
1.0 asQDouble class numBitsInMantissa
Float numBitsInMantissa
ShortFloat numBitsInMantissa
QDouble numBitsInMantissa
"
"Created: / 12-06-2017 / 18:49:11 / cg"
"Modified (comment): / 20-06-2017 / 12:59:00 / cg"
!
radix
"answer the radix of a QDouble's exponent
This is an IEEE float, which is represented as binary"
^ Float radix
"Created: / 12-06-2017 / 18:50:04 / cg"
"Modified (comment): / 19-07-2019 / 17:28:25 / Claus Gittinger"
! !
!QDouble methodsFor:'arithmetic'!
* 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) * 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: / 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
"
"Created: / 12-06-2017 / 16:17:46 / cg"
"Modified: / 12-06-2017 / 23:06:22 / 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_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) - 1.0
((QDouble fromFloat:1e20) - (QDouble fromFloat:1.0)) asDoubleArray
(QDouble fromFloat:1e-20) asDoubleArray
((QDouble fromFloat:1e-20) - (QDouble fromFloat:1.0)) asDoubleArray
((QDouble fromFloat:2.0) - (QDouble fromFloat:1.0)) asDoubleArray
((QDouble fromFloat:2.0) - (QDouble fromFloat:1.0) + (QDouble fromFloat:1.0)) asDoubleArray
((QDouble fromFloat:1e-20) - (QDouble fromFloat:1.0) + (QDouble fromFloat:1.0)) asDoubleArray
"
"Created: / 12-06-2017 / 23:41:39 / cg"
"Modified (comment): / 15-06-2017 / 00:34:41 / cg"
!
/ 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
"
((QDouble fromFloat:1e20) / (QDouble fromFloat:2.0)) asDoubleArray
((QDouble fromFloat:1.2345) / (QDouble fromFloat:10.0)) asDoubleArray
((QDouble fromFloat:1.2345) / 10.0) asDoubleArray
"
"Created: / 13-06-2017 / 17:59:09 / cg"
"Modified (comment): / 15-06-2017 / 00:14:26 / cg"
! !
!QDouble methodsFor:'coercing & converting'!
asDoubleArray
^ DoubleArray
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:2.0) asDoubleArray
"
"Created: / 12-06-2017 / 18:19:19 / cg"
"Modified (comment): / 13-06-2017 / 17:58:09 / cg"
!
asFloat
^ self d0 + self d1
"
(QDouble fromFloat:1.0) asFloat -> 1.0
(QDouble fromFloat:2.0) asFloat -> 2.0
(2.0 asQDouble + 1e-14) asFloat -> 2.00000000000001
(2.0 + 1e-14) - 2.0 -> 1.02140518265514E-14
(2.0 + 1e-15) - 2.0 -> 8.88178419700125E-16
(2.0 + 1e-16) - 2.0 -> 0.0
"
"Created: / 12-06-2017 / 18:15:27 / cg"
"Modified: / 13-06-2017 / 17:56:50 / cg"
!
asInteger
^ self d0 asInteger
+ self d1 asInteger
+ self d2 asInteger
+ self d3 asInteger
"Created: / 19-06-2017 / 18:07:17 / cg"
!
asLargeFloat
^ (self d0 asLargeFloat precision:self precision) + self d1 + self d2 + self d3
"
(QDouble fromFloat:1.0) asLargeFloat -> 1.000000000000000000000000000000
(QDouble fromFloat:2.0) asLargeFloat -> 2.000000000000000000000000000000
(2.0 asQDouble + 1e-14) asLargeFloat -> 2.000000000000010214051826551440
(2.0 asLargeFloat + 1e-14) - 2.0 -> 0.000000000000010214051826551440
(2.0 + 1e-14) - 2.0 -> 1.02140518265514E-14
(2.0 asLargeFloat + 1e-14) - 2.0 -> 0.000000000000010214051826551440
(2.0 asLargeFloat + 1e-15) - 2.0 -> 0.000000000000000888178419700125
(2.0 asLargeFloat + 1e-16) - 2.0 -> 0.0
(2QL + 1QL-14) - 2QL -> 0.000000000000010000000000000000
"
!
asLongFloat
^ self d0 asLongFloat + self d1
"
(QDouble fromFloat:1.0) asLongFloat -> 1.0
(QDouble fromFloat:2.0) asLongFloat -> 2.0
(2.0 asQDouble + 1e-14) asLongFloat -> 2.00000000000001
(2.0 asLongFloat + 1e-14) - 2.0 -> 1.00000303177028016E-14
(2.0 + 1e-14) - 2.0 -> 1.02140518265514E-14
(2.0 asLargeFloat + 1e-14) - 2.0 -> 0.000000000000010214051826551440
(2.0 asLargeFloat + 1e-15) - 2.0 -> 0.000000000000000888178419700125
(2.0 asLargeFloat + 1e-16) - 2.0 -> 0.0
(2QL + 1QL-14) - 2QL -> 0.000000000000010000000000000000
"
"Created: / 12-06-2017 / 18:15:27 / cg"
"Modified: / 13-06-2017 / 17:56:50 / cg"
!
asQDouble
"return a QDouble with same value as myself."
^ self
!
asTrueFraction
^ self d0 asTrueFraction
+ self d1 asTrueFraction
+ self d2 asTrueFraction
+ self d3 asTrueFraction
"
1e10 asTrueFraction -> 10000000000
1e20 asTrueFraction -> 100000000000000000000
(1e20 + 1) asTrueFraction -> 100000000000000000000 ouch!!
1e10 asQDouble asTrueFraction -> 10000000000
1e20 asQDouble asTrueFraction -> 100000000000000000000
(1e20 asQDouble + 1) asTrueFraction -> 100000000000000000001
(1e40 asQDouble + 1e20 + 1) asTrueFraction -> 10000000000000000303886028427003666890753
(1e40 asQDouble + 1e20) asTrueFraction
"
"Created: / 20-06-2017 / 11:09:03 / cg"
!
coerce:aNumber
"convert the argument aNumber into an instance of the receiver's class and return it."
^ aNumber asQDouble
"Created: / 12-06-2017 / 17:13:47 / cg"
"Modified: / 12-06-2017 / 21:09:06 / cg"
!
exponent
"extract a normalized float's (unbiased) exponent.
The returned value depends on the float-representation of
the underlying machine and is therefore highly unportable.
This is not for general use.
This assumes that the mantissa is normalized to
0.5 .. 1.0 and the float's value is: mantissa * 2^exp"
^ self d0 exponent
"Created: / 20-06-2017 / 11:06:02 / cg"
!
generality
"return the generality value - see ArithmeticValue>>retry:coercing:"
^ 95
"Created: / 12-06-2017 / 17:13:14 / cg"
!
mantissa
"extract a normalized float's mantissa.
The returned value depends on the float-representation of
the underlying machine and is therefore highly unportable.
This is not for general use.
This assumes that the mantissa is normalized to
0.5 .. 1.0 and the float's value is mantissa * 2^exp"
"/ fake it here
^ self / (2 raisedTo:self exponent)
"
1.0 exponent -> 1
1.0 mantissa -> 0.5
12345.0 exponent -> 14
12345.0 mantissa -> 0.75347900390625
-1.0 exponent -> 1
-1.0 mantissa -> -0.5
-12345.0 exponent -> 14
-12345.0 mantissa -> -0.75347900390625
(1e40 + 1e-40) exponent -> 133
(1e40 + 1e-40) mantissa -> 0.918354961579912
1.0 asQDouble exponent -> 1
1.0 asQDouble mantissa -> 0.5
12345.0 asQDouble exponent -> 14
12345.0 asQDouble mantissa -> 0.75347900390625
-1.0 asQDouble exponent -> 1
-1.0 asQDouble mantissa -> -0.5
-12345.0 asQDouble exponent -> 14
-12345.0 asQDouble mantissa -> -0.75347900390625
(1e40 + 1e-40) asQDouble exponent -> 133
(1e40 + 1e-40) asQDouble mantissa -> 0.918354961579912
"
"Created: / 20-06-2017 / 11:06:02 / cg"
! !
!QDouble methodsFor:'comparing'!
< aNumber
"return true, if the argument, aNumber is greater than the receiver"
^ aNumber lessFromQDouble:self
"Created: / 13-06-2017 / 16:58:53 / cg"
!
= aNumber
"return true, if the argument, aNumber has the same value as than the receiver"
%{
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 (__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 (__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
"
"Created: / 13-06-2017 / 17:12:09 / cg"
! !
!QDouble methodsFor:'double dispatching'!
differenceFromFloat:aFloat
%{
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)
1e20 - (QDouble fromFloat:1.0)
(1.0 - (QDouble fromFloat:1.0)) asFloat
(1e20 - (QDouble fromFloat:1.0)) asFloat
(1.0 - (QDouble fromFloat:1.0)) asDoubleArray
(1e20 - (QDouble fromFloat:1.0)) asDoubleArray
(1e20 - (QDouble fromFloat:1.0) + 1e-20) asDoubleArray
"
!
differenceFromQDouble: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_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)
((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
"
!
equalFromQDouble:aQDouble
%{
if (__Class(aQDouble) == QDouble) {
double *a = __QDoubleInstPtr(self)->d_qDoubleValue;
double *b = __QDoubleInstPtr(aQDouble)->d_qDoubleValue;
RETURN ((a[0] == b[0]
&& a[1] == b[1]
&& a[2] == b[2]
&& a[3] == b[3]) ? true : false);
}
%}.
^ (aQDouble d0 = self d0)
and:[ (aQDouble d1 = self d1)
and:[ (aQDouble d2 = self d2)
and:[ (aQDouble d3 = self d3) ]]]
"
(QDouble fromFloat:1.0) = (QDouble fromFloat:1.0)
(QDouble fromFloat:1.0) = 1.0
1.0 = (QDouble fromFloat:1.0)
"
"Created: / 13-06-2017 / 03:01:19 / cg"
"Modified: / 13-06-2017 / 18:01:52 / cg"
!
lessFromQDouble:aQDouble
"sent when aQDouble does not know how to compare to the receiver..
Return true if aQDouble < self"
%{
if (__Class(aQDouble) == QDouble) {
double *a = __QDoubleInstPtr(aQDouble)->d_qDoubleValue;
double *b = __QDoubleInstPtr(self)->d_qDoubleValue;
// now compare if a < b!
RETURN
((a[0] < b[0] ||
(a[0] == b[0] && (a[1] < b[1] ||
(a[1] == b[1] && (a[2] < b[2] ||
(a[2] == b[2] && a[3] < b[3])))))) ? true : false);
}
%}.
^ super lessFromQDouble:aQDouble
"
(1.0 + 1e-40) > 1.0
((QDouble fromFloat:1.0) + (QDouble fromFloat:1e-40)) > (QDouble fromFloat:1.0)
(QDouble fromFloat:1.0) > (QDouble fromFloat:1.0)
(QDouble fromFloat:1.1) > (QDouble fromFloat:1.0)
(QDouble fromFloat:1.0) > 1.0
(QDouble fromFloat:1.1) > 1.0
1.0 > (QDouble fromFloat:1.0)
"
"Created: / 13-06-2017 / 17:07:47 / cg"
!
productFromFloat: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_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.
"
loosing bits here:
(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:
((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)
(QDouble fromFloat:2.0) * (QDouble fromFloat:3.0)
QDouble ln2 DoubleArray(0.693147180559945 2.3190468138463E-17 5.70770843841621E-34 -3.58243221060181E-50)
2.0 * QDouble ln2 DoubleArray(1.38629436111989 4.6380936276926E-17 1.14154168768324E-33 -7.16486442120362E-50)
QDouble ln2 * 2.0
2.0 * ((QDouble fromFloat:1e20) + (QDouble fromFloat:1.0)) DoubleArray(2E+20 2.0 0.0 0.0)
((QDouble fromFloat:1e20) + (QDouble fromFloat:1.0)) * 2.0 DoubleArray(2E+20 4E+20 0.0 0.0)
((QDouble fromFloat:1e20) + (QDouble fromFloat:1.0)) * (QDouble fromFloat:2.0) DoubleArray(2E+20 4E+20 0.0 0.0)
(QDouble fromFloat:2.0) * ((QDouble fromFloat:1e20) + (QDouble fromFloat:1.0)) DoubleArray(2E+20 4E+20 0.0 0.0)
(2.0 * ((QDouble fromFloat:1e20) + (QDouble fromFloat:1.0))) - (QDouble fromFloat:1e20) - (QDouble fromFloat:1e20)
(2.0 * (QDouble fromFloat:1.0)) asFloat
(1e20 * (QDouble fromFloat:1.0)) asFloat
(1e20 * (QDouble fromFloat:1.0) * 1e-20) asDoubleArray
"
"Created: / 13-06-2017 / 00:58:56 / cg"
"Modified: / 19-06-2017 / 16:48:18 / cg"
"Modified (comment): / 19-06-2017 / 18:11:43 / cg"
!
productFromQDouble: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_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.
"
(QDouble fromFloat:1.0) * 2.0
2.0 * (QDouble fromFloat:1.0)
(QDouble fromFloat:1.0) * (QDouble fromFloat:2.0)
1e20 * (QDouble fromFloat:2.0)
2.0 * (QDouble fromFloat:1e20)
(QDouble fromFloat:1e20) * (QDouble fromFloat:1e20)
(1e20 * (QDouble fromFloat:1.0) * 1e-20) asDoubleArray
( ((QDouble fromFloat:1.0) + (QDouble fromFloat:1e20)) * (QDouble fromFloat:2.0)) asDoubleArray
"
"Created: / 13-06-2017 / 01:06:22 / cg"
"Modified: / 05-07-2017 / 11:07:16 / cg"
!
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)
1e20 / (QDouble fromFloat:2.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
"
"Created: / 13-06-2017 / 17:50:35 / cg"
"Modified (comment): / 15-06-2017 / 01:02:05 / cg"
!
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
(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
"
"Created: / 13-06-2017 / 17:50:35 / cg"
"Modified (comment): / 15-06-2017 / 01:02:05 / cg"
!
sumFromFloat:aFloat
%{
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);
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 );
}
%}.
^ super sumFromFloat:aFloat.
"
1.0 + (QDouble fromFloat:1.0)
1e20 + (QDouble fromFloat:1.0)
(1.0 + (QDouble fromFloat:1.0)) asFloat
(1e20 + (QDouble fromFloat:1.0)) asFloat
(1.0 + (QDouble fromFloat:1.0)) asDoubleArray
(1e20 + (QDouble fromFloat:1.0)) asDoubleArray
(1e20 + (QDouble fromFloat:1.0) + 1e-20) asDoubleArray
"
"Created: / 12-06-2017 / 17:16:41 / cg"
"Modified: / 14-06-2017 / 11:43:47 / cg"
!
sumFromInteger:anInteger
^ self sumFromFloat:(anInteger asFloat)
"
1 + (QDouble fromFloat:1.0)
1e20 asInteger + (QDouble fromFloat:1.0)
(1 + (QDouble fromFloat:1.0)) asFloat
(1e20 asInteger + (QDouble fromFloat:1.0)) asFloat
"
"Created: / 03-07-2017 / 10:35:46 / cg"
!
sumFromQDouble:aQDouble
%{
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
"
(QDouble fromFloat:1.0) + (QDouble fromFloat:1.0)
(QDouble fromFloat:1.0) + 1.0
1.0 + (QDouble fromFloat:1.0)
((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
"
"Created: / 12-06-2017 / 21:15:43 / cg"
"Modified: / 03-07-2017 / 23:09:11 / cg"
! !
!QDouble methodsFor:'inspecting'!
inspectorExtraAttributes
"extra (pseudo instvar) entries to be shown in an inspector."
^ super inspectorExtraAttributes
add:'-{doubles}' -> [ self asDoubleArray printString ];
yourself
"Created: / 12-06-2017 / 23:43:05 / cg"
"Modified (format): / 18-07-2017 / 19:54:48 / cg"
! !
!QDouble methodsFor:'mathematical functions'!
cos
%{
double *a = __QDoubleInstPtr(self)->d_qDoubleValue;
double q0, q1, q2, q3;
OBJ newQD;
int savedCV;
fpu_fix_start(&savedCV);
qd_cos(&q0, &q1, &q2, &q3, &a[0], &a[1], &a[2], &a[3]);
fpu_fix_end(&savedCV);
__qNew_qdReal(newQD, q0, q1, q2, q3);
RETURN( newQD );
%}.
"
1.0 cos
(QDouble fromFloat:1.0) cos
"
!
exp
%{
double *a = __QDoubleInstPtr(self)->d_qDoubleValue;
double q0, q1, q2, q3;
OBJ newQD;
int savedCV;
fpu_fix_start(&savedCV);
qd_exp(&q0, &q1, &q2, &q3, &a[0], &a[1], &a[2], &a[3]);
fpu_fix_end(&savedCV);
__qNew_qdReal(newQD, q0, q1, q2, q3);
RETURN( newQD );
%}.
"
1.0 exp
(QDouble fromFloat:1.0) exp
3.0 exp
(QDouble fromFloat:3.0) exp
"
!
ldexp:exp
"multiply the receiver by an integral power of 2.
I.e. return self * (2 ^ exp).
This is also the operation to reconstruct the original float from its
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)
"
|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
1.0 ldexp:100 -> 1.26765060022823E+30
1.0 asQDouble ldexp:100 -> 1.26765060022823E+30
"
"Created: / 19-06-2017 / 01:43:35 / cg"
!
ln
"return the natural logarithm of myself.
Raises an exception, if the receiver is less or equal to zero.
Not sure if this is really faster than using a taylor right away:
the three exp-computations at the end are done in qDouble and are tailors themself..."
|d0 x|
"/ ^ super ln.
d0 := self d0.
d0 = 1.0 ifTrue:[
"/ note: d0 checking alone is not sufficient - there could still be more in d1...
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 negated exp)) - 1.0.
x := x + (self * (x negated exp)) - 1.0.
x := x + (self * (x negated exp)) - 1.0.
^ x
].
"/ now done via trapInfinity; was:
"/ d0 = 0.0 ifTrue:[
"/ ^ Infinity negative
"/ ].
"/ 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)'
"
-1 ln
-1.0 asQDouble ln
0.0 asQDouble ln
1.0 asQDouble ln
0.5 ln
0.5 asQDouble ln
3.0 ln printfPrintString:'%60.58lf'
-> 1.0986122886681097821082175869378261268138885498046875000000'
^
3.0 asQDouble ln printfPrintString:'%60.58f'
-> 1.0986122886681096913952452369225257046474905578227494517347
3.0 asQDouble ln printfPrintString:'%70.68f'
-> 1.09861228866810969139524523692252570464749055782274945173469433364779
(3.0 asQDouble ln_withAccuracy:1e-64) printfPrintString:'%70.68f'
1.09861228866810969139524523692252570464749055782274945173469433364475
(3.0 asQDouble ln_withAccuracy:1e-100) printfPrintString:'%70.68f'
'1.098612288668109691395245236922525704647490557822749451734694333656909'
actual result:
-> 1.0986122886681096913952452369225257046474905578227494517346943336374942932186089668736157548137320887879700290659...
"
"Created: / 18-06-2017 / 23:32:54 / cg"
"Modified: / 04-07-2017 / 11:46:27 / cg"
!
negated
^ self class
d0:(self d0) negated
d1:(self d1) negated
d2:(self d2) negated
d3:(self d3) negated
"
(QDouble fromFloat:1.0) negated
((QDouble fromFloat:1e20) + (QDouble fromFloat:1.0)) negated asDoubleArray
(((QDouble fromFloat:1e20) + (QDouble fromFloat:1.0))
+ ((QDouble fromFloat:1e20) + (QDouble fromFloat:1.0))) asDoubleArray
"
"Created: / 12-06-2017 / 20:14:55 / cg"
"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);
qd_pow(&q0, &q1, &q2, &q3, a[0], a[1], a[2], a[3], __intVal(n));
fpu_fix_end(&savedCV);
__qNew_qdReal(newQD, q0, q1, q2, q3);
RETURN( newQD );
}
%}.
^ 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);
qd_sin(&q0, &q1, &q2, &q3, &a[0], &a[1], &a[2], &a[3]);
fpu_fix_end(&savedCV);
__qNew_qdReal(newQD, q0, q1, q2, q3);
RETURN( newQD );
%}.
"
1.0 sin
(QDouble fromFloat:1.0) sin
"
!
sqrt
"Return the square root of the receiver"
%{
double *a = __QDoubleInstPtr(self)->d_qDoubleValue;
double q0, q1, q2, q3;
OBJ newQD;
int savedCV;
fpu_fix_start(&savedCV);
qd_sqrt(&q0, &q1, &q2, &q3, a[0], a[1], a[2], a[3]);
fpu_fix_end(&savedCV);
__qNew_qdReal(newQD, q0, q1, q2, q3);
RETURN( newQD );
%}.
"
(QDouble fromFloat:4.0) sqrt
(QDouble fromFloat:2.0) sqrt
(QDouble fromFloat:1e20) sqrt
"
!
squared
"return receiver * receiver"
%{
double *a = __QDoubleInstPtr(self)->d_qDoubleValue;
double q0, q1, q2, q3;
OBJ newQD;
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 );
%}.
"
(QDouble fromFloat:4.0) squared
(1e20 + (QDouble fromFloat:1.0)) squared
"
"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);
qd_tan(&q0, &q1, &q2, &q3, &a[0], &a[1], &a[2], &a[3]);
fpu_fix_end(&savedCV);
__qNew_qdReal(newQD, q0, q1, q2, q3);
RETURN( newQD );
%}.
"
1.0 tan
(QDouble fromFloat:1.0) tan
"
! !
!QDouble methodsFor:'printing & storing'!
digitsWithPrecision:precision
<resource: #obsolete>
"generate digits and exponent.
if precision is >0, that many digits are generated.
If it is 0 the required number of digits is generated
(but never more than the decimalPrecision, which is 65)"
|numDigits r exp i d out str|
numDigits := precision+1. "/ number of digits
r := self abs.
self d0 = 0.0 ifTrue:[
^ { String new:(precision max:1) withAll:$0 . 0 }
].
out := WriteStream on:(String new:precision+5).
"/ determine approx. exponent
exp := self d0 abs log10 floor.
exp < -300 ifTrue:[
"/ 1e-305 asQDouble
r := r * (10.0 raisedToInteger:300).
r := r / (10.0 raisedToInteger:(exp+300)).
] ifFalse:[
exp > 300 ifTrue:[
"/ 1e305 asQDouble
"/ lexpr(x,exp) = x * 2 ^ exp
self halt.
r := r * (2 raisedTo:-53).
r := r / (10.0 asQDouble raisedTo: exp).
r := r * (2 raisedTo:53).
] ifFalse:[
r := r / (10.0 asQDouble raisedTo:exp).
]
].
"/ Fix exponent if we are off by one
(r >= 10.0) ifTrue:[
r := r / 10.0.
exp := exp + 1.
] ifFalse:[
(r < 1.0) ifTrue:[
r := r * 10.0.
exp := exp - 1.
]
].
((r >= 10.0) or:[ r < 1.0 ]) ifTrue:[
self error:'can''t compute exponent.'.
].
"/
"/ Extract the digits
"/ notice, that the d1,d2 and d3 components might
"/ be negative; therefore characters out of the 0..9 range
"/ might be produced here
"/
i := 1.
[ (precision ~~ 0 and:[ i <= numDigits ])
or:[ (precision == 0 and:[r d0 ~= 0.0]) ]] whileTrue:[
d := r d0 truncated.
r := r - d.
r := r * 10.0.
out nextPut:($0 + d).
i := i + 1.
].
numDigits := i-1.
str := out contents.
"/ Fix out-of-range digits.
numDigits to:2 by:-1 do:[:i |
(str at:i) < $0 ifTrue:[
str at:i-1 put:(str at:i-1) - 1.
str at:i put:(str at:i) + 10.
] ifFalse:[
(str at:i) > $9 ifTrue:[
str at:i-1 put:(str at:i-1) + 1.
str at:i put:(str at:i) - 10.
].
].
].
str first <= $0 ifTrue:[
self error:'non-positive leading digit'
].
"/ Round, handle carry
(str at:numDigits) >= $5 ifTrue:[
str at:numDigits-1 put:(str at:numDigits-1) + 1.
i := numDigits-1.
[i > 1 and:[(str at:i) > $9]] whileTrue:[
str at:i put:(str at:i) - 10.
i := i - 1.
str at:i put:(str at:i) + 1.
]
].
"/ If first digit is 10, shift everything.
str first > $9 ifTrue:[
exp := exp + 1.
str at:1 put:$0.
str := '1',str
].
^ { (str copyTo:numDigits-1) . exp }
"
0 asQDouble digitsWithPrecision:1 -> #('0' 0)
0 asQDouble digitsWithPrecision:0 -> #('0' 0)
1.2345 printfPrintString:'%.4f'
1.2345 asQDouble digitsWithPrecision:5 -> #('12345' 0)
--- but 1.2345 is not really what you think:
1.2345 printfPrintString:'%.20f'
1.2345 asQDouble digitsWithPrecision:20 -> #('12344999999999999307' 0)
12.345 asQDouble digitsWithPrecision:5 -> #('12345' 1)
12345 asQDouble digitsWithPrecision:5 -> #('12345' 4)
12345.1 asQDouble digitsWithPrecision:5 -> #('12345' 4)
12345.9 asQDouble digitsWithPrecision:5 -> #('12346' 4)
1.2345 asQDouble / 10.0 asQDouble
1.2345 asQDouble / 10.0
"
"Created: / 15-06-2017 / 09:10:01 / cg"
"Modified: / 16-06-2017 / 10:01:03 / cg"
!
printOn:aStream
"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"
"/ self d1 = 0.0 ifTrue:[
"/ self d0 printOn:aStream.
"/ ^ self
"/ ].
thisContext isRecursive ifTrue:[
aStream nextPutAll:'aQDouble (error while printing)'.
^ self.
].
PrintfScanf printf:'%g' on:aStream argument:self.
"/ self
"/ printOn:aStream precision:40 width:0
"/ 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 asShortFloat) 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'
"
"Created: / 15-06-2017 / 01:51:36 / cg"
"Modified (comment): / 21-06-2017 / 09:55:10 / cg"
"Modified: / 05-06-2019 / 20:38:58 / Claus Gittinger"
!
printOn:aStream precision:precisionIn width:width
fixed:fixed showPositive:showPositive uppercase:uppercase fillChar:fillChar
<resource: #obsolete>
"return a printed representation of the receiver.
This is a parametrized entry, which can be used by printf-like functions.
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
"
"
1.2345 asQDouble printString
12.345 asQDouble printString
12345 asQDouble printString
"
|sgn count delta exp precision|
"/ self d1 = 0.0 ifTrue:[
"/ self d0 printOn:aStream.
"/ ^ self.
"/ ].
count := 0.
sgn := true.
exp := 0.
precision := precisionIn.
self isInfinite ifTrue:[
self < 0 ifTrue:[
aStream nextPut:$-.
count := 1.
] ifFalse:[
showPositive ifTrue:[
aStream nextPut:$+.
count := 1.
] ifFalse:[
sgn := false.
].
].
uppercase ifTrue:[
aStream nextPutAll:'INF'
] ifFalse:[
aStream nextPutAll:'inf'
].
count := count + 3.
] ifFalse:[
self isNaN ifTrue:[
uppercase ifTrue:[
aStream nextPutAll:'NAN'
] ifFalse:[
aStream nextPutAll:'nan'
].
count := count + 3.
] ifFalse:[
self < 0 ifTrue:[
aStream nextPut:$-.
count := count + 1.
] ifFalse:[
showPositive ifTrue:[
aStream nextPut:$+.
count := count + 1.
] ifFalse:[
sgn := false.
].
].
self = 0.0 ifTrue:[
aStream nextPut:$0.
count := count + 1.
precision > 0 ifTrue:[
aStream nextPut:$..
count := count + 1.
precision timesRepeat:[ aStream nextPut:$0 ].
count := count + precision.
].
self halt.
] ifFalse:[
|off d d_width_extra|
"/ non-zero case
off := fixed
ifTrue:[ 1 + self asFloat abs log10 floor asInteger ]
ifFalse:[1].
d := precision + off.
d_width_extra := d.
fixed ifTrue:[
d_width_extra := 40 max:d.
].
"/ highly special case - fixed mode, precision is zero, abs(*this) < 1.0
"/ without this trap a number like 0.9 printed fixed with 0 precision prints as 0
"/ should be rounded to 1.
(fixed and:[ (precision == 0) and:[ (self abs < 1.0) ]]) ifTrue:[
(self abs >= 0.5) ifTrue:[
aStream nextPut:$1
] ifFalse:[
aStream nextPut:$0
].
^ self
].
"/ handle near zero to working precision (but not exactly zero)
(fixed and:[ d <= 0 ]) ifTrue:[
aStream nextPut:$0.
(precision > 0) ifTrue:[
aStream nextPut:$. .
aStream next:precision put:$0.
]
] ifFalse:[
"/ default
|t j|
t := self digitsWithPrecision:(fixed ifTrue:[d_width_extra] ifFalse:[d])+1.
exp := t second.
t := t first.
fixed ifTrue:[
"/ fix the string if it's been computed incorrectly
"/ round here in the decimal string if required
t := self round_string_qd:t at:(d + 1) offset:off.
precision := t at:3.
off := t at:2.
t := t at:1.
(off > 0) ifTrue:[
aStream next:off putAll:t startingAt:1.
(precision > 0) ifTrue:[
aStream nextPut:$. .
aStream next:precision-1 putAll:t startingAt:off+1.
]
] ifFalse:[
aStream nextPutAll:'0.'.
(off < 0) ifTrue:[
aStream next:off negated put:$0.
].
aStream next:d putAll:t startingAt:0.
]
] ifFalse:[
aStream nextPut:(t at:1).
(precision > 0) ifTrue:[
aStream nextPut:$. .
].
aStream next:precision putAll:t startingAt:2.
]
].
].
]
].
"/ trap for improper offset with large values
"/ without this trap, output of values of the for 10^j - 1 fail for j > 28
"/ and are output with the point in the wrong place, leading to a dramatically off value
"/ (fixed and:[ (precision > 0) ]) ifTrue:[
"/ "/ make sure that the value isn't dramatically larger
"/ from_string = atof(s.c_str());
"/
"/ // if this ratio is large, then we've got problems
"/ if( fabs( from_string / this->x[0] ) > 3.0 ){
"/
"/ int point_position;
"/ char temp;
"/
"/ // loop on the string, find the point, move it up one
"/ // don't act on the first character
"/ for(i=1; i < s.length(); i++){
"/ if(s[i] == '.'){
"/ s[i] = s[i-1] ;
"/ s[i-1] = '.' ;
"/ break;
"/ }
"/ }
"/
"/ from_string = atof(s.c_str());
"/ // if this ratio is large, then the string has not been fixed
"/ if( fabs( from_string / this->x[0] ) > 3.0 ){
"/ dd_real::error("Re-rounding unsuccessful in large number fixed point trap.") ;
"/ }
"/ }
"/ }
"/
fixed ifFalse:[
"/ Fill in exponent part
aStream nextPut:(uppercase ifTrue:[$E] ifFalse:[$e]).
aStream print:exp.
].
"/ fill in the blanks
(delta := width-count) > 0 ifTrue:[
self halt.
"/ if (fmt & ios_base::internal) {
"/ if (sgn)
"/ s.insert(static_cast<string::size_type>(1), delta, fill);
"/ else
"/ s.insert(static_cast<string::size_type>(0), delta, fill);
"/ } else if (fmt & ios_base::left) {
"/ s.append(delta, fill);
"/ } else {
"/ s.insert(static_cast<string::size_type>(0), delta, fill);
"/ }
"/ }
].
"Created: / 15-06-2017 / 02:37:31 / cg"
"Modified (comment): / 16-06-2017 / 14:48:30 / cg"
!
round_string_qd:str at:precisionIn offset:offsetIn
<resource: #obsolete>
"returns a triple of: { new-str . new-offset . new-precision }"
"/
"/ Input string must be all digits or errors will occur.
"/
|i numDigits offsetOut precisionOut|
numDigits := precisionIn.
offsetOut := offsetIn.
precisionOut := precisionIn.
"/ Round, handle carry
((str at:numDigits) >= $5) ifTrue:[
str at:numDigits-1 put:(str at:numDigits-1)+1.
i := numDigits-1.
[ i > 1 and:[ (str at:i) > $9] ] whileTrue:[
str at:i put:(str at:i)-10.
i := i - 1.
str at:i put:(str at:i)+1.
]
].
"/ If first digit is 10, shift everything.
(str at:1) > $9 ifTrue:[
"/ e++; // don't modify exponent here
str replaceFrom:2 with:str startingAt:1.
str at:1 put:$1.
str at:2 put:$0.
offsetOut := offsetOut + 1.
precisionOut := precisionOut + 1.
].
^ { (str copyTo:precisionOut) . offsetOut . precisionOut }
"Created: / 16-06-2017 / 10:12:39 / cg"
"Modified (comment): / 16-06-2017 / 11:22:03 / cg"
! !
!QDouble methodsFor:'private'!
nintAsFloat
"return the receiver truncated towards negative infinity"
%{
/* Computes the nearest integer to d. */
#define nint(d) (((d) == floor(d)) ? (d) : floor((d) + 0.5))
double *a = __QDoubleInstPtr(self)->d_qDoubleValue;
OBJ newQD;
double x0, x1, x2, x3;
x0 = nint(a[0]);
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;
}
}
} else {
/* First double is not an integer. */
if (abs(x0 - a[0]) == 0.5 && a[1] < 0.0) {
x0 -= 1.0;
}
}
renorm(&x0, &x1, &x2, &x3, x0, x1, x2, x3, 0.0);
// m_renorm4(x0, x1, x2, x3);
__qNew_qdReal(newQD, x0, x1, x2, x3);
RETURN( newQD );
%}.
"
(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"
"Modified (comment): / 13-06-2017 / 17:33:19 / cg"
!
renorm
"destructive renormalization"
%{
double *a = __QDoubleInstPtr(self)->d_qDoubleValue;
double c0, c1, c2, c3;
c0 = a[0];
c1 = a[1];
c2 = a[2];
c3 = a[3];
m_renorm4(c0, c1, c2, c3);
a[0] = c0;
a[1] = c1;
a[2] = c2;
a[3] = c3;
RETURN( self );
%}.
^ self error.
"
(QDouble fromFloat:1.0) renorm
"
"Created: / 13-06-2017 / 18:05:33 / cg"
"Modified: / 15-06-2017 / 00:12:59 / cg"
! !
!QDouble methodsFor:'private accessing'!
d0
"the most significant (and highest valued) 53 bits of precision"
%{
RETURN ( __MKFLOAT(__QDoubleInstPtr(self)->d_qDoubleValue[0]) );
%}
"Created: / 12-06-2017 / 20:15:12 / cg"
"Modified (comment): / 13-06-2017 / 17:59:47 / cg"
!
d1
"the next most significant (and next highest valued) 53 bits of precision"
%{
RETURN ( __MKFLOAT(__QDoubleInstPtr(self)->d_qDoubleValue[1]) );
%}
"Created: / 12-06-2017 / 20:15:12 / cg"
"Modified (comment): / 13-06-2017 / 18:00:00 / cg"
!
d2
%{
RETURN ( __MKFLOAT(__QDoubleInstPtr(self)->d_qDoubleValue[2]) );
%}
"Created: / 12-06-2017 / 20:15:29 / cg"
!
d3
"the least significant (and smallest valued) 53 bits of precision"
%{
RETURN ( __MKFLOAT(__QDoubleInstPtr(self)->d_qDoubleValue[3]) );
%}
"Created: / 12-06-2017 / 20:15:32 / cg"
"Modified (comment): / 13-06-2017 / 18:00:18 / cg"
! !
!QDouble methodsFor:'testing'!
isFinite
"return true, if the receiver is a finite float (not NaN and not +/-INF)"
^ self d0 isFinite
"Created: / 17-06-2017 / 03:40:30 / cg"
!
isInfinite
"return true, if the receiver is an infinite float (+Inf or -Inf)."
^ self d0 isInfinite
"Created: / 15-06-2017 / 01:57:57 / cg"
!
isNaN
"return true, if the receiver is an invalid float (NaN - not a number)"
^ self d0 isNaN
"Created: / 15-06-2017 / 01:57:35 / cg"
!
isOne
^ self d0 = 1.0
and:[ self d1 = 0.0
and:[ self d2 = 0.0
and:[ self d3 = 0.0 ]]]
"Created: / 18-06-2017 / 23:29:07 / cg"
!
isZero
^ self d0 = 0.0
"Created: / 18-06-2017 / 23:29:25 / cg"
!
negative
^ self d0 negative
"
(QDouble fromFloat:0.0) negative
(QDouble fromFloat:1.0) negative
(QDouble fromFloat:-1.0) negative
"
"Created: / 13-06-2017 / 01:57:39 / cg"
"Modified: / 13-06-2017 / 17:58:26 / cg"
!
positive
"return true, if the receiver is greater or equal to zero (not negative)"
^ self d0 positive
"
(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
"
"Created: / 13-06-2017 / 01:56:53 / cg"
"Modified: / 13-06-2017 / 17:58:41 / cg"
"Modified (comment): / 28-05-2019 / 05:55:55 / Claus Gittinger"
!
sign
"return the sign of the receiver"
^ 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
"
! !
!QDouble methodsFor:'truncation & rounding'!
ceiling
"return the smallest integer which is greater or equal to the receiver."
|f|
f := self ceilingAsFloat.
^ 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
(1e20 + (QDouble fromFloat:1.0)) ceiling
(1e20 + (QDouble fromFloat:1.1)) ceiling
(QDouble fromFloat:1.5) ceiling
(QDouble fromFloat:0.5) ceiling
(QDouble fromFloat:-0.5) ceiling
(QDouble fromFloat:-1.5) ceiling
"
!
ceilingAsFloat
"return the smallest integer-valued float greater or equal to the receiver.
This is much like #ceiling, but avoids a (possibly expensive) conversion
of the result to an integer.
It may be useful, if the result is to be further used in another float-operation."
%{
double *a = __QDoubleInstPtr(self)->d_qDoubleValue;
OBJ newQD;
int savedCV;
double x0, x1, x2, x3;
x1 = x2 = x3 = 0.0;
x0 = ceil(a[0]);
if (x0 == a[0]) {
x1 = ceil(a[1]);
if (x1 == a[1]) {
x2 = ceil(a[2]);
if (x2 == a[2]) {
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 );
%}.
"
(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
(QDouble fromFloat:1.5) ceiling
(QDouble fromFloat:0.5) ceiling
(QDouble fromFloat:-0.5) ceiling
(QDouble fromFloat:-1.5) ceiling
"
!
floor
"return the receiver truncated towards negative infinity"
|f|
f := self floorAsFloat.
^ f d0 asInteger + f d1 asInteger + f d2 asInteger + f d3 asInteger
"
(QDouble fromFloat:4.0) floor
(QDouble fromFloat:4.1) floor
(QDouble fromFloat:0.1) floor
(0.1 + (QDouble fromFloat:1.0)) floor
(1e20 + (QDouble fromFloat:1.0)) floor
(QDouble fromFloat:1.5) floor
(QDouble fromFloat:0.5) floor
(QDouble fromFloat:-0.5) floor
(QDouble fromFloat:-1.5) floor
"
"Created: / 13-06-2017 / 01:52:44 / cg"
"Modified (comment): / 13-06-2017 / 17:33:19 / cg"
!
floorAsFloat
"return the receiver truncated towards negative infinity"
%{
double *a = __QDoubleInstPtr(self)->d_qDoubleValue;
OBJ newQD;
int savedCV;
double x0, x1, x2, x3;
x1 = x2 = x3 = 0.0;
x0 =floor(a[0]);
if (x0 == a[0]) {
x1 = floor(a[1]);
if (x1 == a[1]) {
x2 = floor(a[2]);
if (x2 == a[2]) {
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 );
%}.
"
(QDouble fromFloat:4.0) floor
(QDouble fromFloat:4.1) floor
(QDouble fromFloat:0.1) floor
(0.1 + (QDouble fromFloat:1.0)) floor
(1e20 + (QDouble fromFloat:1.0)) floor
(QDouble fromFloat:1.5) floor
(QDouble fromFloat:0.5) floor
(QDouble fromFloat:-0.5) floor
(QDouble fromFloat:-1.5) floor
"
"Created: / 13-06-2017 / 01:52:44 / cg"
"Modified (comment): / 13-06-2017 / 17:33:19 / cg"
!
rounded
"return the smallest integer which is greater or equal to the receiver."
|f|
f := self roundedAsFloat.
"/ ^ (f d0 + f d1 + f d2 + f d3) asInteger
^ 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
"
!
roundedAsFloat
"return the receiver truncated towards negative infinity"
self positive ifTrue:[
^ 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
"
"Created: / 13-06-2017 / 01:52:44 / cg"
"Modified (comment): / 13-06-2017 / 17:33:19 / cg"
! !
!QDouble class methodsFor:'documentation'!
version
^ '$Header$'
!
version_CVS
^ '$Header$'
! !