diff -r 663ccfb1701b -r f61106094adc QDouble.st --- 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... "