--- a/QDouble.st Mon Nov 25 16:01:46 2019 +0100
+++ b/QDouble.st Tue Nov 26 00:55:46 2019 +0100
@@ -939,17 +939,16 @@
(returns approx. 200 bits of precision)"
E isNil ifTrue:[
- E := self fromDoubleArray:(DoubleArray
- with: 2.718281828459045091e+00
- with: 1.445646891729250158e-16
- with: -2.127717108038176765e-33
- with: 1.515630159841218954e-49)
+ E := self d0: 2.718281828459045091e+00
+ d1: 1.445646891729250158e-16
+ d2: -2.127717108038176765e-33
+ d3: 1.515630159841218954e-49
].
^ E
"
- self e printfPrintString:'%.60f'
- -> '2.718281828459045235360287471352662497757247093699959574966968'
+ self e printfPrintString:'%.61f'
+ -> '2.7182818284590452353602874713526624977572470936999595749669676'
Wolfram says:
2.71828182845904523536028747135266249775724709369995957496696762772407663035354759457138217852516642742746
"
@@ -962,11 +961,10 @@
(returns approx. 200 bits of precision)"
FMax isNil ifTrue:[
- FMax := self fromDoubleArray:(DoubleArray
- with: 1.797693134862314E+308
- with: 9.97920154767359795037e+291
- with: 5.53956966280111259858e+275
- with: 3.07507889307840487279e+259)
+ FMax := self d0: 1.797693134862314E+308
+ d1: 9.97920154767359795037e+291
+ d2: 5.53956966280111259858e+275
+ d3: 3.07507889307840487279e+259
].
^ FMax
@@ -1053,16 +1051,18 @@
(returns approx. 200 bits of precision)"
Ln10 isNil ifTrue:[
- Ln10 := self fromDoubleArray:(DoubleArray
- with: 2.302585092994045901e+00
- with: -2.170756223382249351e-16
- with: -9.984262454465776570e-33
- with: -4.023357454450206379e-49)
+ Ln10 := self d0: 2.302585092994045901e+00
+ d1: -2.170756223382249351e-16
+ d2: -9.984262454465776570e-33
+ d3: -4.023357454450206379e-49
].
^ Ln10
"
- self ln10
+ self ln10 printfPrintString:'%.61f'
+ -> '2.3025850929940456840179914546843642076011014886287729760333279'
+ Wolfram says:
+ 2.30258509299404568401799145468436420760110148862877297603332790096757260967735248023599720508959829834196778404228...
"
"Created: / 12-06-2017 / 18:32:29 / cg"
@@ -1073,16 +1073,18 @@
(returns approx. 200 bits of precision)"
Ln2 isNil ifTrue:[
- Ln2 := self fromDoubleArray:(DoubleArray
- with: 6.931471805599452862e-01
- with: 2.319046813846299558e-17
- with: 5.707708438416212066e-34
- with: -3.582432210601811423e-50)
+ Ln2 := self d0: 6.931471805599452862e-01
+ d1: 2.319046813846299558e-17
+ d2: 5.707708438416212066e-34
+ d3: -3.582432210601811423e-50
].
^ Ln2
"
- self ln2
+ self ln2 printfPrintString:'%.61f'
+ -> '0.6931471805599452709398341558750792990469129794959648865081141'
+ Wolfram says:
+ 0.69314718055994530941723212145817656807550013436025525412068000949339362196969471560586332699641868754200148102057...
"
"Created: / 12-06-2017 / 18:31:34 / cg"
@@ -1099,17 +1101,16 @@
(returns approx. 200 bits of precision)"
Pi isNil ifTrue:[
- Pi := self fromDoubleArray:(DoubleArray
- with: 3.141592653589793116e+00
- with: 1.224646799147353207e-16
- with: -2.994769809718339666e-33
- with: 1.112454220863364332e-49 "1.112454220863365282e-49")
+ Pi := self d0: 3.141592653589793116e+00
+ d1: 1.224646799147353207e-16
+ d2: -2.994769809718339666e-33
+ d3: 1.112454220863365282e-49
].
^ Pi
"
- self pi printfPrintString:'%.60f'
- -> '3.141592653589793238462643383279502884197169399375447934857836'
+ self pi printfPrintString:'%.60f'
+ '3.141592653589793238462643383279502884197169399375105820974945'
Wolfram says:
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068
@@ -1935,6 +1936,195 @@
"Modified (comment): / 15-06-2017 / 01:02:05 / cg"
!
+quotientFromQDouble_accurate:aQDouble
+ |q0 q1 q2 q3 q4 r|
+
+%{
+/* Computes fl(a+b) and err(a+b). Assumes |a| >= |b|. */
+#define quick_two_sum(s, a, b, err) \
+ { \
+ s = a + b; \
+ err = b - (s - a); \
+ }
+
+#define renorm5(c0, c1, c2, c3, c4) \
+ { \
+ double s0, s1, s2 = 0.0, s3 = 0.0; \
+ \
+ quick_two_sum(s0, c3, c4, c4); \
+ quick_two_sum(s0, c2, s0, c3); \
+ quick_two_sum(s0, c1, s0, c2); \
+ quick_two_sum(c0, c0, s0, c1); \
+ \
+ s0 = c0; \
+ s1 = c1; \
+ \
+ quick_two_sum(s0, c0, c1, s1); \
+ if (s1 != 0.0) { \
+ quick_two_sum(s1, s1, c2, s2); \
+ if (s2 != 0.0) { \
+ quick_two_sum(s2, s2, c3, s3); \
+ if (s3 != 0.0) { \
+ s3 += c4; \
+ } else { \
+ s2 += c4; \
+ } \
+ } else { \
+ quick_two_sum(s1, s1, c3, s2); \
+ if (s2 != 0.0) { \
+ quick_two_sum(s2, s2, c4, s3); \
+ } else { \
+ quick_two_sum(s1, s1, c4, s2); \
+ } \
+ } \
+ } else { \
+ quick_two_sum(s0, s0, c2, s1); \
+ if (s1 != 0.0) { \
+ quick_two_sum(s1, s1, c3, s2); \
+ if (s2 != 0.0) { \
+ quick_two_sum(s2, s2, c4, s3); \
+ } else { \
+ quick_two_sum(s1, s1, c4, s2); \
+ } \
+ } else { \
+ quick_two_sum(s0, s0, c3, s1); \
+ if (s1 != 0.0) { \
+ quick_two_sum(s1, s1, c4, s2); \
+ } else { \
+ quick_two_sum(s0, s0, c4, s1); \
+ } \
+ } \
+ } \
+ \
+ c0 = s0; \
+ c1 = s1; \
+ c2 = s2; \
+ c3 = s3; \
+ }
+%}.
+ q0 := aQDouble d0 / self d0.
+ r := aQDouble - (self * q0).
+
+ q1 := r d0 / self d0.
+ r := r - (self * q1).
+
+ q2 := r d0 / self d0.
+ r := r - (self * q2).
+
+ q3 := r d0 / self d0.
+ r := r - (self * q3).
+
+ q4 := r d0 / self d0.
+ q0 isFinite ifTrue:[
+%{
+ double cq0, cq1, cq2, cq3, cq4;
+ cq0 = __floatVal(q0);
+ cq1 = __floatVal(q1);
+ cq2 = __floatVal(q2);
+ cq3 = __floatVal(q3);
+ cq4 = __floatVal(q4);
+ renorm5(cq0, cq1, cq2, cq3, cq4)
+ {
+ OBJ newQD;
+ __qNew_qdReal(newQD, cq0, cq1, cq2, cq3);
+ RETURN (newQD);
+ }
+#undef quick_two_sum
+#undef renorm5
+%}.
+ ].
+ ^ self class d0:q0 d1:0.0 d2:0.0 d3:0.0
+
+ "
+ 2.0 / (QDouble fromFloat:2.0)
+ 2.0 / (QDouble fromFloat:1.0)
+ 1e20 / (QDouble fromFloat:1.0)
+ (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
+ "
+!
+
+quotientFromQDouble_sloppy:aQDouble
+ "sloppy"
+
+ |q0 q1 q2 q3 r|
+
+ q0 := aQDouble d0 / self d0.
+ "/ Stdout showCR:('q0: %1 (a[0]=%2; b[0]=%3)\n' bindWith:q0 with:self d0 with:aQDouble d0).
+ r := aQDouble - (self * q0).
+ "/ Stdout showCR:('r: %1\n' bindWith:r asDoubleArray).
+
+ q1 := r d0 / self d0.
+ "/ Stdout showCR:('q1: %1 (r[0]=%2; b[0]=%3)\n' bindWith:q1 with:r d0 with:aQDouble d0).
+ r := r - (self * q1).
+ "/ Stdout showCR:('r: %1\n' bindWith:r asDoubleArray).
+
+ q2 := r d0 / self d0.
+ "/ Stdout showCR:('q2: %1 (r[0]=%2; b[0]=%3)\n' bindWith:q2 with:r d0 with:aQDouble d0).
+ r := r - (self * q2).
+ "/ Stdout showCR:('r: %1\n' bindWith:r asDoubleArray).
+
+ q3 := r d0 / self d0.
+ "/ Stdout showCR:('q3: %1 (r[0]=%2; b[0]=%3)\n' bindWith:q3 with:r d0 with:aQDouble d0).
+
+ r := self class d0:q0 d1:q1 d2:q2 d3:q3.
+ "/ Stdout showCR:('before renorm: %1\n' bindWith:r asDoubleArray).
+ r renorm.
+ "/ Stdout showCR:('after renorm: %1\n' bindWith:r asDoubleArray).
+ ^ r
+
+ "
+ 2.0 / (QDouble fromFloat:2.0)
+ 2.0 / (QDouble fromFloat:1.0)
+ 1e20 / (QDouble fromFloat:1.0)
+ (2.0 / (QDouble fromFloat:1.0)) asFloat
+ (1e20 / (QDouble fromFloat:1.0)) asFloat
+
+ (QDouble fromFloat:2.0) / 2.0
+ (QDouble fromFloat:1e20) / 2.0
+ ((QDouble fromFloat:1.0) / 2.0) asFloat
+ ((QDouble fromFloat:1e20 / 2.0)) asFloat
+
+ ((1e20 + (QDouble fromFloat:1.0) + 1e-20) / 2.0) asDoubleArray
+
+ ((QDouble fromFloat:10.0) quotientFromQDouble: (QDouble fromFloat:1.234)) asDoubleArray
+ ((QDouble fromFloat:1.234) / (QDouble fromFloat:10.0)) asDoubleArray
+
+q0: 1.234000e-01 (a[0]=1.234000e+00; b[0]=1.000000e+01)
+a: 1.234000e+00/0.000000e+00/0.000000e+00/0.000000e+00
+b: 1.000000e+01/0.000000e+00/0.000000e+00/0.000000e+00
+(b * q0): 1.234000e+00/-2.775558e-17/0.000000e+00/0.000000e+00
+r: 2.775558e-17/0.000000e+00/0.000000e+00/0.000000e+00
+
+q1: 2.775558e-18 (r[0]=2.775558e-17; b[0]=1.000000e+01)
+r: -1.540744e-33/0.000000e+00/0.000000e+00/0.000000e+00
+
+q2: -1.540744e-34 (r[0]=-1.540744e-33; b[0]=1.000000e+01)
+r: 8.552847e-50/0.000000e+00/0.000000e+00/0.000000e+00
+
+q3: 8.552847e-51 (r[0]=8.552847e-50; b[0]=1.000000e+01)
+
+before renorm: 1.234000e-01/2.775558e-18/-1.540744e-34/8.552847e-51
+after renorm: 1.234000e-01/2.775558e-18/-1.540744e-34/8.552847e-51
+1.234/10.0 is: 0.123400 / 0.000000 / -0.000000 / 0.000000
+
+
+ "
+
+ "Created: / 13-06-2017 / 17:50:35 / cg"
+ "Modified (comment): / 15-06-2017 / 01:02:05 / cg"
+!
+
sumFromFloat:aFloat
%{
if (__isFloatLike(aFloat)) {
@@ -2214,22 +2404,22 @@
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.
+ k := 65536.0. "/ 1.0 ldexp:16.
inv_k := 1.0 / k.
m := (d0 / Float ln2 + 0.5) floor.
@@ -2244,10 +2434,10 @@
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.
+ 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.
@@ -2272,7 +2462,16 @@
"
1.0 exp -> 2.71828182845905
- 1.0 asQDouble exp -> 2.71828182845905
+ 1.0 asQDouble exp -> 2.718281828459045235360287471353
+
+ 10.0 exp -> 22026.4657948067
+ 10.0 asQDouble exp -> 22026.46579480671651695790064528
+
+ Wolfram:
+ 22026.46579480671651695790064528424436635351261855678107423...
+
+ 1000.0 exp -> INF
+ 1000.0 asQDouble exp -> 22026.46579480671651695790064528
"
"Created: / 21-06-2017 / 13:48:27 / cg"
@@ -2476,7 +2675,7 @@
2 asQDouble sqrt printfPrintString:'%70.68f'
-> 1.41421356237309504880168872420969807856967187537694807317667973799602
- actual digits:
+ actual digits (wolfram):
-> 1.4142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276415727350138462309...
"