--- a/QDouble.st Thu Jun 06 11:29:19 2019 +0200
+++ b/QDouble.st Thu Jun 06 11:29:30 2019 +0200
@@ -176,35 +176,35 @@
#define m_two_sum(rslt_3, a_3, b_3, err_3)\
{\
- double s_3 = a_3 + b_3;\
- double bb_3 = s_3 - a_3;\
- err_3 = (a_3 - (s_3 - bb_3)) + (b_3 - bb_3);\
- rslt_3 = s_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;\
+ 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); \
+ 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;\
+ 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
@@ -213,18 +213,21 @@
#define m_split(a_7, hi_7, lo_7)\
{\
double temp_7;\
- if (a_7 > _QD_SPLIT_THRESH || a_7 < -_QD_SPLIT_THRESH) {\
- a_7 *= 3.7252902984619140625e-09;\
- temp_7 = _QD_SPLITTER * a_7;\
- hi_7 = temp_7 - (temp_7 - a_7);\
- lo_7 = a_7 - hi_7;\
- hi_7 *= 268435456.0;\
- lo_7 *= 268435456.0;\
+ 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;\
- hi_7 = temp_7 - (temp_7 - a_7);\
- lo_7 = a_7 - hi_7;\
+ 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
@@ -236,8 +239,8 @@
#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);\
+ double p_8 = (a_8) * (b_8);\
+ err_8 = QD_FMS((a_8), (b_8), p_8);\
rslt_8 = p_8; \
}
@@ -246,7 +249,7 @@
#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;\
+ 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;\
@@ -259,8 +262,8 @@
#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);\
+ double p_9 = (a_9) * (a_9);\
+ err_9 = QD_FMS((a_9), (a_9), p_9);\
rslt_9 = p_9;\
}
@@ -269,7 +272,7 @@
#define m_two_sqr(rslt_9, a_9, err_9)\
{\
double hi_9, lo_9;\
- double q_9 = a_9 * a_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;\
@@ -636,74 +639,74 @@
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, the final precision is less than 53*4)
+ QDoubles use 4 IEEE doubles, each keeping 53 bits of precision.
+ A qDouble's value is the sum of those 4 doubles,
+ and a qDouble keeps this unevaluated sum as its state.
+ (due to overlap, the final precision is less than 53*4)
Range and Precision of Storage Formats: see LimitedPrecisionReal >> documentation
the number of decmal digits:
- QDouble decimalPrecision
- LongFloat decimalPrecision
- Float decimalPrecision
- ShortFloat decimalPrecision
-
+ QDouble decimalPrecision
+ LongFloat decimalPrecision
+ Float decimalPrecision
+ ShortFloat decimalPrecision
+
the number of bits:
- QDouble precision
- LongFloat precision
- Float precision
- ShortFloat precision
-
- 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 then, yet):
- qd := QDouble readFrom:'0.1'.
-
- To see the error of the double precision version, compute:
- (0.1 asQDouble) - (QDouble readFrom:'0.1')
-
+ QDouble precision
+ LongFloat precision
+ Float precision
+ ShortFloat precision
+
+ 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 then, yet):
+ qd := QDouble readFrom:'0.1'.
+
+ To see the error of the double precision version, compute:
+ (0.1 asQDouble) - (QDouble readFrom:'0.1')
+
[author:]
- Claus Gittinger
+ Claus Gittinger
[see also:]
- Number
- Float ShortFloat LongFloat
- Fraction FixedPoint Integer Complex
- FloatArray DoubleArray
+ Number
+ Float ShortFloat LongFloat
+ Fraction FixedPoint Integer Complex
+ FloatArray DoubleArray
"
!
examples
"
Floats, LongFloats suffer from loosing bits:
-
+
(Float readFrom:'0.3333333333333333333333333333333333333333333333333333333333')
-(Float readFrom:'0.333333333333333333333333333333333333333333333333333333333')
- -> 0.0
-
+ -> 0.0
+
(Float readFrom:'0.3333333333333333333333333333333333333333333333333333333333')
= (Float readFrom:'0.333333333333333333333333333333333333333333333333333333333')
- -> true
+ -> true
(Float readFrom:'0.33333333333333333333333333333333333333333333333333333333333333333333')
= (Float readFrom:'0.3333333333333333333333333333333333333333333333333333333333333333333')
- -> true
+ -> 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
+ -> 0.0
(LongFloat readFrom:'0.3333333333333333333333333333333333333333333333333333333333')
= (LongFloat readFrom:'0.333333333333333333333333333333333333333333333333333333333')
- -> 0.0
+ -> 0.0
(QDouble readFrom:'0.3333333333333333333333333333333333333333333333333333333333')
-(QDouble readFrom:'0.333333333333333333333333333333333333333333333333333333333')
@@ -861,28 +864,28 @@
OBJ newFloat;
if (__isSmallInteger(anInteger)) {
- INT iVal = __smallIntegerVal(anInteger);
- double *a;
-
- __qNew(newFloat, sizeof(struct __quadDoubleStruct));
- __stx_setClass(newFloat, QDouble);
-
- a = __QuadDoubleInstPtr(newFloat)->d_quadDoubleValue;
- a[1] = 0.0;
- a[2] = 0.0;
- a[3] = 0.0;
-
- // need more than 52bits?
- if ((sizeof(INT) > 52)
- && ((iVal > 0xFFFFFFFFFFFFF)
- || (iVal < -0xFFFFFFFFFFFFF))) {
- a[0] = (double)(iVal & ~0xFFFFFFFF);
- a[1] = (double)(iVal & 0xFFFFFFFF);
- // m_renorm4(a[0], a[1], a[2], a[3]);
- } else {
- a[0] = (double)iVal;
- }
- RETURN (newFloat);
+ INT iVal = __smallIntegerVal(anInteger);
+ double *a;
+
+ __qNew(newFloat, sizeof(struct __quadDoubleStruct));
+ __stx_setClass(newFloat, QDouble);
+
+ a = __QuadDoubleInstPtr(newFloat)->d_quadDoubleValue;
+ a[1] = 0.0;
+ a[2] = 0.0;
+ a[3] = 0.0;
+
+ // need more than 52bits?
+ if ((sizeof(INT) > 52)
+ && ((iVal > 0xFFFFFFFFFFFFF)
+ || (iVal < -0xFFFFFFFFFFFFF))) {
+ a[0] = (double)(iVal & ~0xFFFFFFFF);
+ a[1] = (double)(iVal & 0xFFFFFFFF);
+ // m_renorm4(a[0], a[1], a[2], a[3]);
+ } else {
+ a[0] = (double)iVal;
+ }
+ RETURN (newFloat);
}
#endif
%}.
@@ -922,7 +925,7 @@
"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 := self d0:(Float NaN) d1:(Float NaN) d2:(Float NaN) d3:(Float NaN)
].
^ NaN
@@ -995,45 +998,45 @@
invFact
"table returning 1/n!!
(for taylor series)"
-
+
InvFact isNil ifTrue:[
- InvFact := Array new:15.
- InvFact at:1 put:(self d0:1.66666666666666657e-01 d1:9.25185853854297066e-18 d2:5.13581318503262866e-34 d3:2.85094902409834186e-50).
- InvFact at:2 put:(self d0:4.16666666666666644e-02 d1:2.31296463463574266e-18 d2:1.28395329625815716e-34 d3:7.12737256024585466e-51).
- InvFact at:3 put:(self d0:8.33333333333333322e-03 d1:1.15648231731787138e-19 d2:1.60494162032269652e-36 d3:2.22730392507682967e-53).
- InvFact at:4 put:(self d0:1.38888888888888894e-03 d1:-5.30054395437357706e-20 d2:-1.73868675534958776e-36 d3:-1.63335621172300840e-52).
- InvFact at:5 put:(self d0:1.98412698412698413e-04 d1:1.72095582934207053e-22 d2:1.49269123913941271e-40 d3:1.29470326746002471e-58).
- InvFact at:6 put:(self d0:2.48015873015873016e-05 d1:2.15119478667758816e-23 d2:1.86586404892426588e-41 d3:1.61837908432503088e-59).
- InvFact at:7 put:(self d0:2.75573192239858925e-06 d1:-1.85839327404647208e-22 d2:8.49175460488199287e-39 d3:-5.72661640789429621e-55).
- InvFact at:8 put:(self d0:2.75573192239858883e-07 d1:2.37677146222502973e-23 d2:-3.26318890334088294e-40 d3:1.61435111860404415e-56).
- InvFact at:9 put:(self d0:2.50521083854417202e-08 d1:-1.44881407093591197e-24 d2:2.04267351467144546e-41 d3:-8.49632672007163175e-58).
- InvFact at:10 put:(self d0:2.08767569878681002e-09 d1:-1.20734505911325997e-25 d2:1.70222792889287100e-42 d3:1.41609532150396700e-58).
- InvFact at:11 put:(self d0:1.60590438368216133e-10 d1:1.25852945887520981e-26 d2:-5.31334602762985031e-43 d3:3.54021472597605528e-59).
- InvFact at:12 put:(self d0:1.14707455977297245e-11 d1:2.06555127528307454e-28 d2:6.88907923246664603e-45 d3:5.72920002655109095e-61).
- InvFact at:13 put:(self d0:7.64716373181981641e-13 d1:7.03872877733453001e-30 d2:-7.82753927716258345e-48 d3:1.92138649443790242e-64).
- InvFact at:14 put:(self d0:4.77947733238738525e-14 d1:4.39920548583408126e-31 d2:-4.89221204822661465e-49 d3:1.20086655902368901e-65).
- InvFact at:15 put:(self d0:2.81145725434552060e-15 d1:1.65088427308614326e-31 d2:-2.87777179307447918e-50 d3:4.27110689256293549e-67).
+ InvFact := Array new:15.
+ InvFact at:1 put:(self d0:1.66666666666666657e-01 d1:9.25185853854297066e-18 d2:5.13581318503262866e-34 d3:2.85094902409834186e-50).
+ InvFact at:2 put:(self d0:4.16666666666666644e-02 d1:2.31296463463574266e-18 d2:1.28395329625815716e-34 d3:7.12737256024585466e-51).
+ InvFact at:3 put:(self d0:8.33333333333333322e-03 d1:1.15648231731787138e-19 d2:1.60494162032269652e-36 d3:2.22730392507682967e-53).
+ InvFact at:4 put:(self d0:1.38888888888888894e-03 d1:-5.30054395437357706e-20 d2:-1.73868675534958776e-36 d3:-1.63335621172300840e-52).
+ InvFact at:5 put:(self d0:1.98412698412698413e-04 d1:1.72095582934207053e-22 d2:1.49269123913941271e-40 d3:1.29470326746002471e-58).
+ InvFact at:6 put:(self d0:2.48015873015873016e-05 d1:2.15119478667758816e-23 d2:1.86586404892426588e-41 d3:1.61837908432503088e-59).
+ InvFact at:7 put:(self d0:2.75573192239858925e-06 d1:-1.85839327404647208e-22 d2:8.49175460488199287e-39 d3:-5.72661640789429621e-55).
+ InvFact at:8 put:(self d0:2.75573192239858883e-07 d1:2.37677146222502973e-23 d2:-3.26318890334088294e-40 d3:1.61435111860404415e-56).
+ InvFact at:9 put:(self d0:2.50521083854417202e-08 d1:-1.44881407093591197e-24 d2:2.04267351467144546e-41 d3:-8.49632672007163175e-58).
+ InvFact at:10 put:(self d0:2.08767569878681002e-09 d1:-1.20734505911325997e-25 d2:1.70222792889287100e-42 d3:1.41609532150396700e-58).
+ InvFact at:11 put:(self d0:1.60590438368216133e-10 d1:1.25852945887520981e-26 d2:-5.31334602762985031e-43 d3:3.54021472597605528e-59).
+ InvFact at:12 put:(self d0:1.14707455977297245e-11 d1:2.06555127528307454e-28 d2:6.88907923246664603e-45 d3:5.72920002655109095e-61).
+ InvFact at:13 put:(self d0:7.64716373181981641e-13 d1:7.03872877733453001e-30 d2:-7.82753927716258345e-48 d3:1.92138649443790242e-64).
+ InvFact at:14 put:(self d0:4.77947733238738525e-14 d1:4.39920548583408126e-31 d2:-4.89221204822661465e-49 d3:1.20086655902368901e-65).
+ InvFact at:15 put:(self d0:2.81145725434552060e-15 d1:1.65088427308614326e-31 d2:-2.87777179307447918e-50 d3:4.27110689256293549e-67).
].
^ InvFact
"
- 1.0 / (3 factorial) 0.166666666666667
-
- 1 asQDouble / (3 factorial) - (self invFact at:1)
- 1 asQDouble / (4 factorial) - (self invFact at:2)
- 1 asQDouble / (5 factorial) - (self invFact at:3)
- 1 asQDouble / (6 factorial) - (self invFact at:4)
- 1 asQDouble / (7 factorial) - (self invFact at:5)
- 1 asQDouble / (8 factorial) - (self invFact at:6)
- 1 asQDouble / (9 factorial) - (self invFact at:7)
- 1 asQDouble / (10 factorial) - (self invFact at:8)
- 1 asQDouble / (11 factorial) - (self invFact at:9)
- 1 asQDouble / (12 factorial) - (self invFact at:10)
- 1 asQDouble / (13 factorial) - (self invFact at:11)
- 1 asQDouble / (14 factorial) - (self invFact at:12)
- 1 asQDouble / (15 factorial) - (self invFact at:13)
- 1 asQDouble / (16 factorial) - (self invFact at:14)
- 1 asQDouble / (17 factorial) - (self invFact at:15)
+ 1.0 / (3 factorial) 0.166666666666667
+
+ 1 asQDouble / (3 factorial) - (self invFact at:1)
+ 1 asQDouble / (4 factorial) - (self invFact at:2)
+ 1 asQDouble / (5 factorial) - (self invFact at:3)
+ 1 asQDouble / (6 factorial) - (self invFact at:4)
+ 1 asQDouble / (7 factorial) - (self invFact at:5)
+ 1 asQDouble / (8 factorial) - (self invFact at:6)
+ 1 asQDouble / (9 factorial) - (self invFact at:7)
+ 1 asQDouble / (10 factorial) - (self invFact at:8)
+ 1 asQDouble / (11 factorial) - (self invFact at:9)
+ 1 asQDouble / (12 factorial) - (self invFact at:10)
+ 1 asQDouble / (13 factorial) - (self invFact at:11)
+ 1 asQDouble / (14 factorial) - (self invFact at:12)
+ 1 asQDouble / (15 factorial) - (self invFact at:13)
+ 1 asQDouble / (16 factorial) - (self invFact at:14)
+ 1 asQDouble / (17 factorial) - (self invFact at:15)
"
"Created: / 19-06-2017 / 02:22:23 / cg"
@@ -1335,8 +1338,8 @@
asInteger
^ self d0 asInteger
- + self d1 asInteger
- + self d2 asInteger
+ + self d1 asInteger
+ + self d2 asInteger
+ self d3 asInteger
"Created: / 19-06-2017 / 18:07:17 / cg"
@@ -1366,8 +1369,8 @@
1e20 asQDouble asTrueFraction -> 100000000000000000000
(1e20 asQDouble + 1) asTrueFraction -> 100000000000000000001
- (1e40 asQDouble + 1e20 + 1) asTrueFraction -> 10000000000000000303886028427003666890753
- (1e40 asQDouble + 1e20) asTrueFraction
+ (1e40 asQDouble + 1e20 + 1) asTrueFraction -> 10000000000000000303886028427003666890753
+ (1e40 asQDouble + 1e20) asTrueFraction
"
"Created: / 20-06-2017 / 11:09:03 / cg"
@@ -1545,44 +1548,44 @@
productFromFloat:aFloat
%{
if (__isFloatLike(aFloat)) {
- double *a = __QuadDoubleInstPtr(self)->d_quadDoubleValue;
- double b = __floatVal(aFloat);
- double p0, p1, p2, p3;
- double q0, q1, q2;
- double s0, s1, s2, s3, s4;
- OBJ newQD;
- int savedCV;
-
- fpu_fix_start(&savedCV);
-
- m_two_prod(p0, a[0], b, q0);
- m_two_prod(p1, a[1], b, q1);
- m_two_prod(p2, a[2], b, q2);
- p3 = a[3] * b;
-
- s0 = p0;
-
- m_two_sum(s1, q0, p1, s2);
-
- m_three_sum(s2, q1, p2);
-
- m_three_sum2(q1, q2, p3);
- s3 = q1;
-
- s4 = q2 + p2;
-
- m_renorm5(s0, s1, s2, s3, s4);
- fpu_fix_end(&savedCV);
-
- __qNew_qdReal(newQD, s0, s1, s2, s3);
- RETURN( newQD );
+ double *a = __QuadDoubleInstPtr(self)->d_quadDoubleValue;
+ double b = __floatVal(aFloat);
+ double p0, p1, p2, p3;
+ double q0, q1, q2;
+ double s0, s1, s2, s3, s4;
+ OBJ newQD;
+ int savedCV;
+
+ fpu_fix_start(&savedCV);
+
+ m_two_prod(p0, a[0], b, q0);
+ m_two_prod(p1, a[1], b, q1);
+ m_two_prod(p2, a[2], b, q2);
+ p3 = a[3] * b;
+
+ s0 = p0;
+
+ m_two_sum(s1, q0, p1, s2);
+
+ m_three_sum(s2, q1, p2);
+
+ m_three_sum2(q1, q2, p3);
+ s3 = q1;
+
+ s4 = q2 + p2;
+
+ m_renorm5(s0, s1, s2, s3, s4);
+ fpu_fix_end(&savedCV);
+
+ __qNew_qdReal(newQD, s0, s1, s2, s3);
+ RETURN( newQD );
}
%}.
^ super productFromFloat:aFloat.
"
loosing bits here:
-
+
(1e20+1.0)*2.0
(1e20+1.0)*1e20
(1e40+1.0)*2.0
@@ -1590,24 +1593,24 @@
but not here:
(1.0 asQDouble) * 2.0
- ((1e20 asQDouble) + (1.0)) * 2.0
+ ((1e20 asQDouble) + (1.0)) * 2.0
((1e20 asQDouble) + (1.0)) * 100.0 10000000000000000000100.0
- ((1e20 asQDouble) + (1.0)) * 1000.0
- ((1e40 asQDouble) + (1.0)) * 2.0
+ ((1e20 asQDouble) + (1.0)) * 1000.0
+ ((1e40 asQDouble) + (1.0)) * 2.0
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
-
+ 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
@@ -1624,136 +1627,136 @@
productFromQDouble:aQDouble
%{
if (__Class(aQDouble) == QDouble) {
- double *a = __QuadDoubleInstPtr(self)->d_quadDoubleValue;
- double *b = __QuadDoubleInstPtr(aQDouble)->d_quadDoubleValue;
- OBJ newQD;
- int savedCV;
+ double *a = __QuadDoubleInstPtr(self)->d_quadDoubleValue;
+ double *b = __QuadDoubleInstPtr(aQDouble)->d_quadDoubleValue;
+ OBJ newQD;
+ int savedCV;
#define QD_IEEE_MUL
#ifndef QD_IEEE_MUL
- // sloppy
- double p0, p1, p2, p3, p4, p5;
- double q0, q1, q2, q3, q4, q5;
- double t0, t1;
- double s0, s1, s2;
-
- fpu_fix_start(&savedCV);
-
- m_two_prod(p0, a[0], b[0], q0);
- // fprintf(stderr, "%f * %f -> %f, %f\n", a[0], b[0], p0, q0);
-
- m_two_prod(p1, a[0], b[1], q1);
- m_two_prod(p2, a[1], b[0], q2);
- // fprintf(stderr, "%f * %f -> %f, %f\n", a[0], b[1], p1, q1);
- // fprintf(stderr, "%f * %f -> %f, %f\n", a[1], b[0], p2, q2);
-
- m_two_prod(p3, a[0], b[2], q3);
- m_two_prod(p4, a[1], b[1], q4);
- m_two_prod(p5, a[2], b[0], q5);
-
- // fprintf(stderr, "%f * %f -> %f, %f\n", a[0], b[2], p3, q3);
- // fprintf(stderr, "%f * %f -> %f, %f\n", a[1], b[1], p4, q4);
- // fprintf(stderr, "%f * %f -> %f, %f\n", a[2], b[0], p5, q5);
-
- /* Start Accumulation */
- m_three_sum(p1, p2, q0);
- // fprintf(stderr, "mul7: after three_sum: p1:%e p2:%e q0:%e\n", p1, p2, q0);
-
- /* Six-Three Sum of p2, q1, q2, p3, p4, p5. */
- m_three_sum(p2, q1, q2);
- // fprintf(stderr, "mul8: after three_sum: p2:%e q1:%e q2:%e\n", p2, q1, q2);
- m_three_sum(p3, p4, p5);
- // fprintf(stderr, "mul9: after three_sum: p3:%e p4:%e p5:%e\n", p3, p4, p5);
-
- /* compute (s0, s1, s2) = (p2, q1, q2) + (p3, p4, p5). */
- m_two_sum(s0, p2, p3, t0);
- // fprintf(stderr, "mul10: after two_sum: s0:%e p2:%e p3:%e t0:%e\n", s0, p2, p3, t0);
- m_two_sum(s1, q1, p4, t1);
- // fprintf(stderr, "mul11: after two_sum: s1:%e q1:%e p4:%e t1:%e\n", s1, q1, p4, t1);
- s2 = q2 + p5;
- m_two_sum(s1, s1, t0, t0);
- // fprintf(stderr, "mul12: after two_sum: s1:%e s1:%e t0:%e t0:%e\n", s1, s1, t0, t0);
- s2 += (t0 + t1);
-
- /* O(eps^3) order terms */
- s1 += a[0]*b[3] + a[1]*b[2] + a[2]*b[1] + a[3]*b[0] + q0 + q3 + q4 + q5;
-
- // fprintf(stderr, "before renorm5: p0:%e p2:%e s0:%e s1:%e s2:%e\n", p0, p1, s0, s1, s2);
- m_renorm5(p0, p1, s0, s1, s2);
- // fprintf(stderr, "after renorm5: p0:%e p2:%e s0:%e s1:%e s2:%e\n", p0, p1, s0, s1, s2);
- fpu_fix_end(&savedCV);
-
- __qNew_qdReal(newQD, p0, p1, s0, s1);
+ // sloppy
+ double p0, p1, p2, p3, p4, p5;
+ double q0, q1, q2, q3, q4, q5;
+ double t0, t1;
+ double s0, s1, s2;
+
+ fpu_fix_start(&savedCV);
+
+ m_two_prod(p0, a[0], b[0], q0);
+ // fprintf(stderr, "%f * %f -> %f, %f\n", a[0], b[0], p0, q0);
+
+ m_two_prod(p1, a[0], b[1], q1);
+ m_two_prod(p2, a[1], b[0], q2);
+ // fprintf(stderr, "%f * %f -> %f, %f\n", a[0], b[1], p1, q1);
+ // fprintf(stderr, "%f * %f -> %f, %f\n", a[1], b[0], p2, q2);
+
+ m_two_prod(p3, a[0], b[2], q3);
+ m_two_prod(p4, a[1], b[1], q4);
+ m_two_prod(p5, a[2], b[0], q5);
+
+ // fprintf(stderr, "%f * %f -> %f, %f\n", a[0], b[2], p3, q3);
+ // fprintf(stderr, "%f * %f -> %f, %f\n", a[1], b[1], p4, q4);
+ // fprintf(stderr, "%f * %f -> %f, %f\n", a[2], b[0], p5, q5);
+
+ /* Start Accumulation */
+ m_three_sum(p1, p2, q0);
+ // fprintf(stderr, "mul7: after three_sum: p1:%e p2:%e q0:%e\n", p1, p2, q0);
+
+ /* Six-Three Sum of p2, q1, q2, p3, p4, p5. */
+ m_three_sum(p2, q1, q2);
+ // fprintf(stderr, "mul8: after three_sum: p2:%e q1:%e q2:%e\n", p2, q1, q2);
+ m_three_sum(p3, p4, p5);
+ // fprintf(stderr, "mul9: after three_sum: p3:%e p4:%e p5:%e\n", p3, p4, p5);
+
+ /* compute (s0, s1, s2) = (p2, q1, q2) + (p3, p4, p5). */
+ m_two_sum(s0, p2, p3, t0);
+ // fprintf(stderr, "mul10: after two_sum: s0:%e p2:%e p3:%e t0:%e\n", s0, p2, p3, t0);
+ m_two_sum(s1, q1, p4, t1);
+ // fprintf(stderr, "mul11: after two_sum: s1:%e q1:%e p4:%e t1:%e\n", s1, q1, p4, t1);
+ s2 = q2 + p5;
+ m_two_sum(s1, s1, t0, t0);
+ // fprintf(stderr, "mul12: after two_sum: s1:%e s1:%e t0:%e t0:%e\n", s1, s1, t0, t0);
+ s2 += (t0 + t1);
+
+ /* O(eps^3) order terms */
+ s1 += a[0]*b[3] + a[1]*b[2] + a[2]*b[1] + a[3]*b[0] + q0 + q3 + q4 + q5;
+
+ // fprintf(stderr, "before renorm5: p0:%e p2:%e s0:%e s1:%e s2:%e\n", p0, p1, s0, s1, s2);
+ m_renorm5(p0, p1, s0, s1, s2);
+ // fprintf(stderr, "after renorm5: p0:%e p2:%e s0:%e s1:%e s2:%e\n", p0, p1, s0, s1, s2);
+ fpu_fix_end(&savedCV);
+
+ __qNew_qdReal(newQD, p0, p1, s0, s1);
#else
- // accurate
- double p0, p1, p2, p3, p4, p5;
- double q0, q1, q2, q3, q4, q5;
- double p6, p7, p8, p9;
- double q6, q7, q8, q9;
- double r0, r1;
- double t0, t1;
- double s0, s1, s2;
-
- fpu_fix_start(&savedCV);
-
- m_two_prod(p0, a[0], b[0], q0);
-
- m_two_prod(p1, a[0], b[1], q1);
- m_two_prod(p2, a[1], b[0], q2);
-
- m_two_prod(p3, a[0], b[2], q3);
- m_two_prod(p4, a[1], b[1], q4);
- m_two_prod(p5, a[2], b[0], q5);
-
- /* Start Accumulation */
- m_three_sum(p1, p2, q0);
-
- /* Six-Three Sum of p2, q1, q2, p3, p4, p5. */
- m_three_sum(p2, q1, q2);
- m_three_sum(p3, p4, p5);
- /* compute (s0, s1, s2) = (p2, q1, q2) + (p3, p4, p5). */
- m_two_sum(s0, p2, p3, t0);
- m_two_sum(s1, q1, p4, t1);
- s2 = q2 + p5;
- m_two_sum(s1, s1, t0, t0);
- s2 += (t0 + t1);
-
- /* O(eps^3) order terms */
- m_two_prod(p6, a[0], b[3], q6);
- m_two_prod(p7, a[1], b[2], q7);
- m_two_prod(p8, a[2], b[1], q8);
- m_two_prod(p9, a[3], b[0], q9);
-
- /* Nine-Two-Sum of q0, s1, q3, q4, q5, p6, p7, p8, p9. */
- m_two_sum(q0, q0, q3, q3);
- m_two_sum(q4, q4, q5, q5);
- m_two_sum(p6, p6, p7, p7);
- m_two_sum(p8, p8, p9, p9);
- /* Compute (t0, t1) = (q0, q3) + (q4, q5). */
- m_two_sum(t0, q0, q4, t1);
- t1 += (q3 + q5);
- /* Compute (r0, r1) = (p6, p7) + (p8, p9). */
- m_two_sum(r0, p6, p8, r1);
- r1 += (p7 + p9);
- /* Compute (q3, q4) = (t0, t1) + (r0, r1). */
- m_two_sum(q3, t0, r0, q4);
- q4 += (t1 + r1);
- /* Compute (t0, t1) = (q3, q4) + s1. */
- m_two_sum(t0, q3, s1, t1);
- t1 += q4;
-
- /* O(eps^4) terms -- Nine-One-Sum */
- t1 += a[1] * b[3] + a[2] * b[2] + a[3] * b[1] + q6 + q7 + q8 + q9 + s2;
-
- // fprintf(stderr, "before accur-renorm5: p0:%e p1:%e s0:%e t0:%e t1:%e\n", p0, p1, s0, t0, t1);
- m_renorm5(p0, p1, s0, t0, t1);
- // fprintf(stderr, "after accur-renorm5: p0:%e p1:%e s0:%e t0:%e t1:%e\n", p0, p1, s0, t0, t1);
- fpu_fix_end(&savedCV);
-
- __qNew_qdReal(newQD, p0, p1, s0, t0);
-#endif
- RETURN( newQD );
+ // accurate
+ double p0, p1, p2, p3, p4, p5;
+ double q0, q1, q2, q3, q4, q5;
+ double p6, p7, p8, p9;
+ double q6, q7, q8, q9;
+ double r0, r1;
+ double t0, t1;
+ double s0, s1, s2;
+
+ fpu_fix_start(&savedCV);
+
+ m_two_prod(p0, a[0], b[0], q0);
+
+ m_two_prod(p1, a[0], b[1], q1);
+ m_two_prod(p2, a[1], b[0], q2);
+
+ m_two_prod(p3, a[0], b[2], q3);
+ m_two_prod(p4, a[1], b[1], q4);
+ m_two_prod(p5, a[2], b[0], q5);
+
+ /* Start Accumulation */
+ m_three_sum(p1, p2, q0);
+
+ /* Six-Three Sum of p2, q1, q2, p3, p4, p5. */
+ m_three_sum(p2, q1, q2);
+ m_three_sum(p3, p4, p5);
+ /* compute (s0, s1, s2) = (p2, q1, q2) + (p3, p4, p5). */
+ m_two_sum(s0, p2, p3, t0);
+ m_two_sum(s1, q1, p4, t1);
+ s2 = q2 + p5;
+ m_two_sum(s1, s1, t0, t0);
+ s2 += (t0 + t1);
+
+ /* O(eps^3) order terms */
+ m_two_prod(p6, a[0], b[3], q6);
+ m_two_prod(p7, a[1], b[2], q7);
+ m_two_prod(p8, a[2], b[1], q8);
+ m_two_prod(p9, a[3], b[0], q9);
+
+ /* Nine-Two-Sum of q0, s1, q3, q4, q5, p6, p7, p8, p9. */
+ m_two_sum(q0, q0, q3, q3);
+ m_two_sum(q4, q4, q5, q5);
+ m_two_sum(p6, p6, p7, p7);
+ m_two_sum(p8, p8, p9, p9);
+ /* Compute (t0, t1) = (q0, q3) + (q4, q5). */
+ m_two_sum(t0, q0, q4, t1);
+ t1 += (q3 + q5);
+ /* Compute (r0, r1) = (p6, p7) + (p8, p9). */
+ m_two_sum(r0, p6, p8, r1);
+ r1 += (p7 + p9);
+ /* Compute (q3, q4) = (t0, t1) + (r0, r1). */
+ m_two_sum(q3, t0, r0, q4);
+ q4 += (t1 + r1);
+ /* Compute (t0, t1) = (q3, q4) + s1. */
+ m_two_sum(t0, q3, s1, t1);
+ t1 += q4;
+
+ /* O(eps^4) terms -- Nine-One-Sum */
+ t1 += a[1] * b[3] + a[2] * b[2] + a[3] * b[1] + q6 + q7 + q8 + q9 + s2;
+
+ // fprintf(stderr, "before accur-renorm5: p0:%e p1:%e s0:%e t0:%e t1:%e\n", p0, p1, s0, t0, t1);
+ m_renorm5(p0, p1, s0, t0, t1);
+ // fprintf(stderr, "after accur-renorm5: p0:%e p1:%e s0:%e t0:%e t1:%e\n", p0, p1, s0, t0, t1);
+ fpu_fix_end(&savedCV);
+
+ __qNew_qdReal(newQD, p0, p1, s0, t0);
+#endif
+ RETURN( newQD );
}
%}.
^ super productFromQDouble:aQDouble.
@@ -1922,153 +1925,153 @@
q3_outRef = s;\
} else {\
if (!zb) {\
- q3_b = q3_a;\
- q3_a = s;\
+ q3_b = q3_a;\
+ q3_a = s;\
} else {\
- q3_a = s;\
+ q3_a = s;\
}\
q3_outRef = 0.0;\
}\
}
if (__Class(aQDouble) == QDouble) {
- double *a = __QuadDoubleInstPtr(self)->d_quadDoubleValue;
- double *b = __QuadDoubleInstPtr(aQDouble)->d_quadDoubleValue;
- OBJ newQD;
- int savedCV;
+ double *a = __QuadDoubleInstPtr(self)->d_quadDoubleValue;
+ double *b = __QuadDoubleInstPtr(aQDouble)->d_quadDoubleValue;
+ OBJ newQD;
+ int savedCV;
#define QD_IEEE_ADD
#ifndef xQD_IEEE_ADD
- // sloppy_add...
+ // sloppy_add...
# if 0
- double s0, s1, s2, s3;
- double t0, t1, t2, t3;
-
- fpu_fix_start(&savedCV);
- m_two_sum(s0, a[0], b[0], t0);
- m_two_sum(s1, a[1], b[1], t1);
- m_two_sum(s2, a[2], b[2], t2);
- m_two_sum(s3, a[3], b[3], t3);
-
- m_two_sum(s1, s1, t0, t0);
- m_three_sum(s2, t0, t1);
- m_three_sum2(s3, t0, t2);
- t0 = t0 + t1 + t3;
-
- m_renorm5(s0, s1, s2, s3, t0);
- fpu_fix_end(&savedCV);
-
- __qNew_qdReal(newQD, s0, s1, s2, s3);
+ double s0, s1, s2, s3;
+ double t0, t1, t2, t3;
+
+ fpu_fix_start(&savedCV);
+ m_two_sum(s0, a[0], b[0], t0);
+ m_two_sum(s1, a[1], b[1], t1);
+ m_two_sum(s2, a[2], b[2], t2);
+ m_two_sum(s3, a[3], b[3], t3);
+
+ m_two_sum(s1, s1, t0, t0);
+ m_three_sum(s2, t0, t1);
+ m_three_sum2(s3, t0, t2);
+ t0 = t0 + t1 + t3;
+
+ m_renorm5(s0, s1, s2, s3, t0);
+ fpu_fix_end(&savedCV);
+
+ __qNew_qdReal(newQD, s0, s1, s2, s3);
# else
- /* Same as above, but addition re-organized to minimize
- data dependency ... unfortunately some compilers are
- not very smart to do this automatically */
- double s0, s1, s2, s3;
- double t0, t1, t2, t3;
- double v0, v1, v2, v3;
- double u0, u1, u2, u3;
- double w0, w1, w2, w3;
-
- fpu_fix_start(&savedCV);
- s0 = a[0] + b[0];
- s1 = a[1] + b[1];
- s2 = a[2] + b[2];
- s3 = a[3] + b[3];
-
- v0 = s0 - a[0];
- v1 = s1 - a[1];
- v2 = s2 - a[2];
- v3 = s3 - a[3];
-
- u0 = s0 - v0;
- u1 = s1 - v1;
- u2 = s2 - v2;
- u3 = s3 - v3;
-
- w0 = a[0] - u0;
- w1 = a[1] - u1;
- w2 = a[2] - u2;
- w3 = a[3] - u3;
-
- u0 = b[0] - v0;
- u1 = b[1] - v1;
- u2 = b[2] - v2;
- u3 = b[3] - v3;
-
- t0 = w0 + u0;
- t1 = w1 + u1;
- t2 = w2 + u2;
- t3 = w3 + u3;
-
- m_two_sum(s1, s1, t0, t0);
- m_three_sum(s2, t0, t1);
- m_three_sum2(s3, t0, t2);
- t0 = t0 + t1 + t3;
-
- /* renormalize */
- m_renorm5(s0, s1, s2, s3, t0);
- fpu_fix_end(&savedCV);
- __qNew_qdReal(newQD, s0, s1, s2, s3);
+ /* Same as above, but addition re-organized to minimize
+ data dependency ... unfortunately some compilers are
+ not very smart to do this automatically */
+ double s0, s1, s2, s3;
+ double t0, t1, t2, t3;
+ double v0, v1, v2, v3;
+ double u0, u1, u2, u3;
+ double w0, w1, w2, w3;
+
+ fpu_fix_start(&savedCV);
+ s0 = a[0] + b[0];
+ s1 = a[1] + b[1];
+ s2 = a[2] + b[2];
+ s3 = a[3] + b[3];
+
+ v0 = s0 - a[0];
+ v1 = s1 - a[1];
+ v2 = s2 - a[2];
+ v3 = s3 - a[3];
+
+ u0 = s0 - v0;
+ u1 = s1 - v1;
+ u2 = s2 - v2;
+ u3 = s3 - v3;
+
+ w0 = a[0] - u0;
+ w1 = a[1] - u1;
+ w2 = a[2] - u2;
+ w3 = a[3] - u3;
+
+ u0 = b[0] - v0;
+ u1 = b[1] - v1;
+ u2 = b[2] - v2;
+ u3 = b[3] - v3;
+
+ t0 = w0 + u0;
+ t1 = w1 + u1;
+ t2 = w2 + u2;
+ t3 = w3 + u3;
+
+ m_two_sum(s1, s1, t0, t0);
+ m_three_sum(s2, t0, t1);
+ m_three_sum2(s3, t0, t2);
+ t0 = t0 + t1 + t3;
+
+ /* renormalize */
+ m_renorm5(s0, s1, s2, s3, t0);
+ fpu_fix_end(&savedCV);
+ __qNew_qdReal(newQD, s0, s1, s2, s3);
# endif
#else
- // ieee_add...
- int i, j, k;
- double s, t;
- double u, v; /* double-length accumulator */
- double x[4] = {0.0, 0.0, 0.0, 0.0};
-
- fpu_fix_start(&savedCV);
-
- i = j = k = 0;
- if (abs(a[i]) > abs(b[j]))
- u = a[i++];
- else
- u = b[j++];
- if (abs(a[i]) > abs(b[j]))
- v = a[i++];
- else
- v = b[j++];
-
- m_quick_two_sum(u, u, v, v);
-
- while (k < 4) {
- if (i >= 4 && j >= 4) {
- x[k] = u;
- if (k < 3)
- x[++k] = v;
- break;
- }
-
- if (i >= 4)
- t = b[j++];
- else if (j >= 4)
- t = a[i++];
- else if (abs(a[i]) > abs(b[j])) {
- t = a[i++];
- } else
- t = b[j++];
-
- m_quick_three_accum(s, u, v, t);
-
- if (s != 0.0) {
- x[k++] = s;
- }
- }
-
- /* add the rest. */
- for (k = i; k < 4; k++)
- x[3] += a[k];
- for (k = j; k < 4; k++)
- x[3] += b[k];
-
- m_renorm4(x[0], x[1], x[2], x[3]);
- fpu_fix_end(&savedCV);
- __qNew_qdReal(newQD, x[0], x[1], x[2], x[3]);
+ // ieee_add...
+ int i, j, k;
+ double s, t;
+ double u, v; /* double-length accumulator */
+ double x[4] = {0.0, 0.0, 0.0, 0.0};
+
+ fpu_fix_start(&savedCV);
+
+ i = j = k = 0;
+ if (abs(a[i]) > abs(b[j]))
+ u = a[i++];
+ else
+ u = b[j++];
+ if (abs(a[i]) > abs(b[j]))
+ v = a[i++];
+ else
+ v = b[j++];
+
+ m_quick_two_sum(u, u, v, v);
+
+ while (k < 4) {
+ if (i >= 4 && j >= 4) {
+ x[k] = u;
+ if (k < 3)
+ x[++k] = v;
+ break;
+ }
+
+ if (i >= 4)
+ t = b[j++];
+ else if (j >= 4)
+ t = a[i++];
+ else if (abs(a[i]) > abs(b[j])) {
+ t = a[i++];
+ } else
+ t = b[j++];
+
+ m_quick_three_accum(s, u, v, t);
+
+ if (s != 0.0) {
+ x[k++] = s;
+ }
+ }
+
+ /* add the rest. */
+ for (k = i; k < 4; k++)
+ x[3] += a[k];
+ for (k = j; k < 4; k++)
+ x[3] += b[k];
+
+ m_renorm4(x[0], x[1], x[2], x[3]);
+ fpu_fix_end(&savedCV);
+ __qNew_qdReal(newQD, x[0], x[1], x[2], x[3]);
#endif
- RETURN(newQD);
+ RETURN(newQD);
}
%}.
^ super sumFromQDouble:aQDouble
@@ -2095,8 +2098,8 @@
"extra (pseudo instvar) entries to be shown in an inspector."
^ super inspectorExtraAttributes
- add:'-{doubles}' -> [ self asDoubleArray printString ];
- yourself
+ add:'-{doubles}' -> [ self asDoubleArray printString ];
+ yourself
"Created: / 12-06-2017 / 23:43:05 / cg"
"Modified (format): / 18-07-2017 / 19:54:48 / cg"
@@ -2112,15 +2115,15 @@
^ super exp
"
- 2.0 exp -> 7.38905609893065
+ 2.0 exp -> 7.38905609893065
2.0 asQDouble exp -> 7.38905609893065022723
- 2 asQDouble exp printfPrintString:'%70.68f'
- -> 7.389056098930650227230427460575007813180315570551847324087127822432669
+ 2 asQDouble exp printfPrintString:'%70.68f'
+ -> 7.389056098930650227230427460575007813180315570551847324087127822432669
actual result:
- -> 7.3890560989306502272304274605750078131803155705518473240871278225225737960790577633843124850791217947737531612654...
-
+ -> 7.3890560989306502272304274605750078131803155705518473240871278225225737960790577633843124850791217947737531612654...
+
"
"Created: / 19-06-2017 / 01:49:32 / cg"
@@ -2132,7 +2135,7 @@
"/ however, it is inexact in the 37th digit
"/
"/ the inherited exp code is much slower, but more precise.
-
+
"/ Strategy: We first reduce the size of x by noting that
"/
"/ exp(kr + m * log(2)) = 2^m * exp(r)^k
@@ -2149,19 +2152,19 @@
d0 := self d0.
(d0 <= -709.0) ifTrue:[
- ^ 0.0 asQDouble.
+ ^ 0.0 asQDouble.
].
(d0 >= 709.0) ifTrue:[
- ^ Infinity positive
+ ^ Infinity positive
].
(d0 = 0.0) ifTrue:[
- ^ 1.0 asQDouble.
+ ^ 1.0 asQDouble.
].
(d0 = 1.0) ifTrue:[
- self isOne ifTrue:[
- ^ self class e
- ].
+ self isOne ifTrue:[
+ ^ self class e
+ ].
].
k := 1.0 ldexp:16.
@@ -2179,11 +2182,11 @@
s := r + (mul_pwr2 value:p value:0.5).
i := 1.
[
- p := p * r.
- t := p * (self class invFact at:i).
- i := i+1.
- s := s + t.
- ] doWhile:[ (t asFloat abs > thresh) and:[i < 10] ].
+ p := p * r.
+ t := p * (self class invFact at:i).
+ i := i+1.
+ s := s + t.
+ ] doWhile:[ (t asFloat abs > thresh) and:[i < 10] ].
s := (mul2 value:s) + s squared.
s := (mul2 value:s) + s squared.
@@ -2233,36 +2236,36 @@
approx := self + 1 + (num / den).
[
- facN := facN + 1.
- den := den * facN.
- num := num * self.
-
- delta := num / den.
- "/ Transcript showCR:delta.
- delta isNaN ifTrue:[self halt:'nan when dividing for delta'. num / den].
- delta = 0 ifTrue:[self halt:'zero delta'].
- approx := approx + delta.
+ facN := facN + 1.
+ den := den * facN.
+ num := num * self.
+
+ delta := num / den.
+ "/ Transcript showCR:delta.
+ delta isNaN ifTrue:[self halt:'nan when dividing for delta'. num / den].
+ delta = 0 ifTrue:[self halt:'zero delta'].
+ approx := approx + delta.
] doUntil:[delta abs <= epsilon].
^ approx
"
wolfram:
- 7.389056098930650227230427460575007813180315570551847324087
-
- 2 asQDouble exp_withAccuracy:QDouble epsilon 7.38905609893065022723
- 2 asQDouble exp_withAccuracy:LongFloat epsilon 7.38905609893065022723
+ 7.389056098930650227230427460575007813180315570551847324087
+
+ 2 asQDouble exp_withAccuracy:QDouble epsilon 7.38905609893065022723
+ 2 asQDouble exp_withAccuracy:LongFloat epsilon 7.38905609893065022723
2 asQDouble exp_withAccuracy:Float epsilon 7.38905609893065022489
- 2.718281828459045235360287471352662497757247093699959574966...
-
- 1 asQDouble exp_withAccuracy:QDouble epsilon 2.71828182845904523536
- 1 asQDouble exp_withAccuracy:LongFloat epsilon 2.71828182845904523536
+ 2.718281828459045235360287471352662497757247093699959574966...
+
+ 1 asQDouble exp_withAccuracy:QDouble epsilon 2.71828182845904523536
+ 1 asQDouble exp_withAccuracy:LongFloat epsilon 2.71828182845904523536
1 asQDouble exp_withAccuracy:Float epsilon 2.71828182845904522671
- 1.7392749415205010473946813036112352261479840577250084... × 10^18
-
- 42 asQDouble exp_withAccuracy:QDouble epsilon NAN
+ 1.7392749415205010473946813036112352261479840577250084... × 10^18
+
+ 42 asQDouble exp_withAccuracy:QDouble epsilon NAN
42 asQDouble exp_withAccuracy:LongFloat epsilon 1.73927494152050104739e18
42 asQDouble exp_withAccuracy:Float epsilon 1.73927494152050104739e18
@@ -2372,22 +2375,22 @@
|d0 x|
^ super ln.
-
+
d0 := self d0.
d0 = 1.0 ifTrue:[
- self isOne ifTrue:[ ^ self class zero ].
+ self isOne ifTrue:[ ^ self class zero ].
].
d0 > 0.0 ifTrue:[
- "/ initial approx.
- x := d0 ln asQDouble.
- "/ three more iterations of newton...
- x := x + (self / (x exp)) - 1.0.
- x := x + (self / (x exp)) - 1.0.
- x := x + (self / (x exp)) - 1.0.
-
- ^ x
+ "/ initial approx.
+ x := d0 ln asQDouble.
+ "/ three more iterations of newton...
+ x := x + (self / (x exp)) - 1.0.
+ x := x + (self / (x exp)) - 1.0.
+ x := x + (self / (x exp)) - 1.0.
+
+ ^ x
].
-
+
"/ now done via trapInfinity; was:
"/ d0 = 0.0 ifTrue:[
"/ ^ Infinity negative
@@ -2395,11 +2398,11 @@
"/ if you need -INF for a zero receiver, try Number trapInfinity:[...]
^ self class
- raise:(self = 0 ifTrue:[#infiniteResultSignal] ifFalse:[#domainErrorSignal])
- receiver:self
- selector:#ln
- arguments:#()
- errorString:'bad receiver in ln (not strictly positive)'
+ raise:(self = 0 ifTrue:[#infiniteResultSignal] ifFalse:[#domainErrorSignal])
+ receiver:self
+ selector:#ln
+ arguments:#()
+ errorString:'bad receiver in ln (not strictly positive)'
"
-1 ln
@@ -2408,23 +2411,23 @@
0.0 asQDouble ln
1.0 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 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
+ 1.09861228866810969139524523692252570464749055782274945173469433364475
(3.0 asQDouble ln_withAccuracy:1e-100) printfPrintString:'%70.68f'
- '1.098612288668109691395245236922525704647490557822749451734694333656909'
+ '1.098612288668109691395245236922525704647490557822749451734694333656909'
actual result:
- -> 1.0986122886681096913952452369225257046474905578227494517346943336374942932186089668736157548137320887879700290659...
+ -> 1.0986122886681096913952452369225257046474905578227494517346943336374942932186089668736157548137320887879700290659...
"
"Created: / 18-06-2017 / 23:32:54 / cg"
@@ -2457,28 +2460,28 @@
using sqrt from the double as an initial guess"
|guess|
-
+
guess := self asFloat sqrt asQDouble.
^ self sqrt_withAccuracy:(self epsilon) fromInitialGuess:guess
- "FloatD is only correct in roughly 17 digits
-
- 2 sqrt printfPrintString:'%50.48f'
- -> 1.414213562373095145474621858738828450441360473633
- ^
+ "FloatD is only correct in roughly 17 digits
+
+ 2 sqrt printfPrintString:'%50.48f'
+ -> 1.414213562373095145474621858738828450441360473633
+ ^
QDouble gives you much more:
- 2 asQDouble sqrt printfPrintString:'%50.48f'
- -> 1.414213562373095048801688724209698078569671875377
-
- 2 asQDouble sqrt printfPrintString:'%60.58f'
- -> 1.4142135623730950488016887242096980785696718753769480731767'
-
- 2 asQDouble sqrt printfPrintString:'%70.68f'
- -> 1.41421356237309504880168872420969807856967187537694807317667973799602
+ 2 asQDouble sqrt printfPrintString:'%50.48f'
+ -> 1.414213562373095048801688724209698078569671875377
+
+ 2 asQDouble sqrt printfPrintString:'%60.58f'
+ -> 1.4142135623730950488016887242096980785696718753769480731767'
+
+ 2 asQDouble sqrt printfPrintString:'%70.68f'
+ -> 1.41421356237309504880168872420969807856967187537694807317667973799602
actual digits:
- -> 1.4142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276415727350138462309...
+ -> 1.4142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276415727350138462309...
"
"Created: / 22-06-2017 / 13:53:47 / cg"
@@ -2499,7 +2502,7 @@
OBJ newQD;
// code makes use of some symetries to avoid a few multiplications...
-
+
m_two_sqr(p0, a[0], q0);
t = 2.0 * a[0];
@@ -2687,16 +2690,16 @@
"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"
-
- thisContext isRecursive ifTrue:[
- aStream nextPutAll:'aQDouble (error while printing)'.
- ^ self.
+ 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"
+
+ thisContext isRecursive ifTrue:[
+ aStream nextPutAll:'aQDouble (error while printing)'.
+ ^ self.
].
-
+
PrintfScanf printf:'%g' on:aStream argument:self.
"/ self
@@ -3081,4 +3084,3 @@
version_CVS
^ '$Header$'
! !
-