--- a/LargeFloat.st Tue May 28 17:44:04 2019 +0200
+++ b/LargeFloat.st Tue May 28 17:46:28 2019 +0200
@@ -548,13 +548,14 @@
!
naiveRaisedToInteger: n
- "Very naive algorithm: use full precision.
- Use only for small n"
- | m e |
- m _ mantissa raisedToInteger: n.
- e _ biasedExponent * n.
- ^(m asLargeFloatPrecision: precision) timesTwoPower: e
-
+ "Very naive algorithm: use full precision.
+ Use only for small n"
+ | m e |
+ m := mantissa raisedToInteger: n.
+ e := biasedExponent * n.
+ ^(m asLargeFloatPrecision: precision) timesTwoPower: e
+
+ "Modified (comment): / 28-05-2019 / 17:43:40 / Claus Gittinger"
!
negated
@@ -583,32 +584,34 @@
!
pi
- "answer the value of pi rounded to precision.
- Note: use the Brent-Salamin Arithmetic Geometric Mean algorithm"
-
- | a b c k pi oldpi oldExpo expo |
- a _ self one asLargeFloatPrecision: precision + 16.
- b _ (a timesTwoPower: 1) sqrt reciprocal.
- c _ a timesTwoPower: -1.
- k _ 1.
- oldpi _ Float pi.
- oldExpo _ 2.
-
- [| am gm a2 |
- am _ a + b timesTwoPower: -1.
- gm _ (a * b) sqrt.
- a _ am.
- b _ gm.
- a2 _ a squared.
- c inPlaceSubtract: (a2 - b squared timesTwoPower: k).
- pi _ (a2 timesTwoPower: 1) / c.
- expo _ (oldpi - pi) exponent.
- expo isZero or: [expo > oldExpo or: [expo < (-1 - precision)]]]
- whileFalse:
- [oldpi _ pi.
- oldExpo _ expo.
- k _ k + 1].
- ^pi asLargeFloatPrecision: precision
+ "answer the value of pi rounded to precision.
+ Note: use the Brent-Salamin Arithmetic Geometric Mean algorithm"
+
+ | a b c k pi oldpi oldExpo expo |
+ a := self one asLargeFloatPrecision: precision + 16.
+ b := (a timesTwoPower: 1) sqrt reciprocal.
+ c := a timesTwoPower: -1.
+ k := 1.
+ oldpi := Float pi.
+ oldExpo := 2.
+
+ [| am gm a2 |
+ am := a + b timesTwoPower: -1.
+ gm := (a * b) sqrt.
+ a := am.
+ b := gm.
+ a2 := a squared.
+ c inPlaceSubtract: (a2 - b squared timesTwoPower: k).
+ pi := (a2 timesTwoPower: 1) / c.
+ expo := (oldpi - pi) exponent.
+ expo isZero or: [expo > oldExpo or: [expo < (-1 - precision)]]]
+ whileFalse:
+ [oldpi := pi.
+ oldExpo := expo.
+ k := k + 1].
+ ^pi asLargeFloatPrecision: precision
+
+ "Modified (comment): / 28-05-2019 / 17:43:46 / Claus Gittinger"
!
piDoublePrecision
@@ -1098,83 +1101,92 @@
!LargeFloat methodsFor:'initialization'!
setPrecisionTo: n
- precision _ n.
+ precision := n.
self roundToPrecision
"Modified: / 28-05-2019 / 11:22:21 / Claus Gittinger"
+ "Modified (format): / 28-05-2019 / 17:45:16 / Claus Gittinger"
! !
!LargeFloat methodsFor:'mathematical functions'!
agm: aNumber
- "Answer the arithmetic geometric mean of self and aNumber"
-
- | a b am gm |
- a _ self.
- b _ aNumber.
-
- [am _ a + b timesTwoPower: -1. "am is arithmetic mean"
- gm _ (a * b) sqrt. "gm is geometric mean"
- a = am or: [b = gm]]
- whileFalse:
- [a _ am.
- b _ gm].
- ^am
+ "Answer the arithmetic geometric mean of self and aNumber"
+
+ | a b am gm |
+ a := self.
+ b := aNumber.
+
+ [am := a + b timesTwoPower: -1. "am is arithmetic mean"
+ gm := (a * b) sqrt. "gm is geometric mean"
+ a = am or: [b = gm]]
+ whileFalse:
+ [a := am.
+ b := gm].
+ ^am
+
+ "Modified (format): / 28-05-2019 / 17:42:34 / Claus Gittinger"
!
arCosh
- "Evaluate the area hyperbolic cosine of the receiver."
-
- | arCosh x one y two |
- x _ self asLargeFloatPrecision: 16 + precision.
- one _ x one.
- x < one ifTrue: [DomainError signal: 'cannot compute arCosh of a number less than 1'].
- x = one ifTrue: [^self zero].
- y _ x - one.
- y < one
- ifTrue:
- [y exponent * -4 >= precision
- ifTrue: [arCosh _ (y powerExpansionArCoshp1Precision: y precision) * (y timesTwoPower: 1) sqrt]
- ifFalse:
- [two _ one timesTwoPower: 1.
- arCosh _ ((y * (y + two)) sqrt + y + one) ln]]
- ifFalse: [arCosh _ ((x squared - one) sqrt + x) ln].
- ^arCosh asLargeFloatPrecision: precision
+ "Evaluate the area hyperbolic cosine of the receiver."
+
+ | arCosh x one y two |
+ x := self asLargeFloatPrecision: 16 + precision.
+ one := x one.
+ x < one ifTrue: [DomainError signal: 'cannot compute arCosh of a number less than 1'].
+ x = one ifTrue: [^self zero].
+ y := x - one.
+ y < one
+ ifTrue:
+ [y exponent * -4 >= precision
+ ifTrue: [arCosh := (y powerExpansionArCoshp1Precision: y precision) * (y timesTwoPower: 1) sqrt]
+ ifFalse:
+ [two := one timesTwoPower: 1.
+ arCosh := ((y * (y + two)) sqrt + y + one) ln]]
+ ifFalse: [arCosh := ((x squared - one) sqrt + x) ln].
+ ^arCosh asLargeFloatPrecision: precision
+
+ "Modified (format): / 28-05-2019 / 17:42:40 / Claus Gittinger"
!
arSinh
- "Evaluate the area hyperbolic sine of the receiver."
-
- | arSinh x one |
- self isZero ifTrue: [^self].
- self exponent negated > precision ifTrue: [^self].
- x _ self asLargeFloatPrecision: 16 + precision.
- x inPlaceAbs.
- self exponent * -4 >= precision
- ifTrue: [arSinh _ x powerExpansionArSinhPrecision: x precision]
- ifFalse:
- [one _ x one.
- arSinh _ ((x squared + one) sqrt + x) ln].
- self negative ifTrue: [arSinh inPlaceNegated].
- ^arSinh asLargeFloatPrecision: precision
+ "Evaluate the area hyperbolic sine of the receiver."
+
+ | arSinh x one |
+ self isZero ifTrue: [^self].
+ self exponent negated > precision ifTrue: [^self].
+ x := self asLargeFloatPrecision: 16 + precision.
+ x inPlaceAbs.
+ self exponent * -4 >= precision
+ ifTrue: [arSinh := x powerExpansionArSinhPrecision: x precision]
+ ifFalse:
+ [one := x one.
+ arSinh := ((x squared + one) sqrt + x) ln].
+ self negative ifTrue: [arSinh inPlaceNegated].
+ ^arSinh asLargeFloatPrecision: precision
+
+ "Modified (format): / 28-05-2019 / 17:42:45 / Claus Gittinger"
!
arTanh
- "Evaluate the area hyperbolic tangent of the receiver."
-
- | arTanh x one |
- self isZero ifTrue: [^self].
- x _ self asLargeFloatPrecision: 16 + precision.
- x inPlaceAbs.
- one _ x one.
- x >= one ifTrue: [DomainError signal: 'cannot evaluate arTanh of number of magnitude >= 1'].
- self exponent * -4 >= precision
- ifTrue: [arTanh _ x powerExpansionArTanhPrecision: x precision]
- ifFalse:
- [arTanh _ ((one + x) / (one - x)) ln.
- arTanh inPlaceTimesTwoPower: -1].
- self negative ifTrue: [arTanh inPlaceNegated].
- ^arTanh asLargeFloatPrecision: precision
+ "Evaluate the area hyperbolic tangent of the receiver."
+
+ | arTanh x one |
+ self isZero ifTrue: [^self].
+ x := self asLargeFloatPrecision: 16 + precision.
+ x inPlaceAbs.
+ one := x one.
+ x >= one ifTrue: [DomainError signal: 'cannot evaluate arTanh of number of magnitude >= 1'].
+ self exponent * -4 >= precision
+ ifTrue: [arTanh := x powerExpansionArTanhPrecision: x precision]
+ ifFalse:
+ [arTanh := ((one + x) / (one - x)) ln.
+ arTanh inPlaceTimesTwoPower: -1].
+ self negative ifTrue: [arTanh inPlaceNegated].
+ ^arTanh asLargeFloatPrecision: precision
+
+ "Modified (format): / 28-05-2019 / 17:42:49 / Claus Gittinger"
!
arcCos
@@ -1303,81 +1315,85 @@
!
exp
- "Answer the exponential of the receiver."
-
- | ln2 x q r ri res n maxIter p one two |
- one _ self one.
- two _ one timesTwoPower: 1.
- "Use following decomposition:
- x exp = (2 ln * q + r) exp.
- x exp = (2**q * r exp)"
- ln2 _ two ln.
- x _ self / ln2.
- q _ x truncated.
- r _ (x - q) * ln2.
-
- "now compute r exp by power series expansion
- we compute (r/(2**p)) exp ** (2**p) in order to have faster convergence"
- p _ 10 min: precision // 2.
- r _ r timesTwoPower: p negated.
- ri _ one asLargeFloatPrecision: precision + 16.
- res _ ri copy.
- n _ 0.
- maxIter _ 1 + ((precision + 16) / p) ceiling.
- [n <= maxIter] whileTrue:
- [n _ n + 1.
- ri inPlaceMultiplyBy: r / n.
- res inPlaceAdd: ri].
- p timesRepeat: [res inPlaceMultiplyBy: res].
- res inPlaceTimesTwoPower: q.
-
- "now use a Newton iteration to refine the result
- res = res * (self - res ln + 1)"
- [| oldres delta |
- oldres _ res.
- res _ res asLargeFloatPrecision: res precision + 32.
- res inPlaceMultiplyBy: self - res ln + 1.
- delta _ (res - oldres) exponent.
- delta = 0 or: [delta <= (res exponent - precision - 8)]]
- whileFalse.
-
- ^res asLargeFloatPrecision: precision
+ "Answer the exponential of the receiver."
+
+ | ln2 x q r ri res n maxIter p one two |
+ one := self one.
+ two := one timesTwoPower: 1.
+ "Use following decomposition:
+ x exp = (2 ln * q + r) exp.
+ x exp = (2**q * r exp)"
+ ln2 := two ln.
+ x := self / ln2.
+ q := x truncated.
+ r := (x - q) * ln2.
+
+ "now compute r exp by power series expansion
+ we compute (r/(2**p)) exp ** (2**p) in order to have faster convergence"
+ p := 10 min: precision // 2.
+ r := r timesTwoPower: p negated.
+ ri := one asLargeFloatPrecision: precision + 16.
+ res := ri copy.
+ n := 0.
+ maxIter := 1 + ((precision + 16) / p) ceiling.
+ [n <= maxIter] whileTrue:
+ [n := n + 1.
+ ri inPlaceMultiplyBy: r / n.
+ res inPlaceAdd: ri].
+ p timesRepeat: [res inPlaceMultiplyBy: res].
+ res inPlaceTimesTwoPower: q.
+
+ "now use a Newton iteration to refine the result
+ res = res * (self - res ln + 1)"
+ [| oldres delta |
+ oldres := res.
+ res := res asLargeFloatPrecision: res precision + 32.
+ res inPlaceMultiplyBy: self - res ln + 1.
+ delta := (res - oldres) exponent.
+ delta = 0 or: [delta <= (res exponent - precision - 8)]]
+ whileFalse.
+
+ ^res asLargeFloatPrecision: precision
+
+ "Modified (comment): / 28-05-2019 / 17:43:06 / Claus Gittinger"
!
ln
- "Answer the neperian logarithm of the receiver."
-
- | x4 one two p res selfHighRes prec e |
- self <= self zero ifTrue: [DomainError signal: 'ln is only defined for x > 0.0'].
-
- one _ self one.
- self = one ifTrue: [^self zero].
- two _ one timesTwoPower: 1.
-
- "Use Salamin algorithm (approximation is good if x is big enough)
- x ln = Pi / (2 * (1 agm: 4/x) ).
- If x not big enough, compute (x timesTwoPower: p) ln - (2 ln * p)
- if x is close to 1, better use a power expansion"
- prec _ precision + 16.
- e _ self exponent.
- e < 0 ifTrue: [e _ -1 - e].
- e > prec
- ifTrue: [p _ 0]
- ifFalse:
- [p _ prec - e.
- prec _ prec + p highBit].
- selfHighRes _ self asLargeFloatPrecision: prec.
- (selfHighRes - one) exponent * -4 >= precision ifTrue: [^(selfHighRes powerExpansionLnPrecision: prec) asLargeFloatPrecision: precision].
- self < one ifTrue: [selfHighRes inPlaceReciprocal]. "Use ln(1/x) => - ln(x)"
- x4 _ (4 asLargeFloatPrecision: prec)
- inPlaceDivideBy: selfHighRes;
- inPlaceTimesTwoPower: p negated.
- res _ selfHighRes pi / (one agm: x4) timesTwoPower: -1.
- res _ selfHighRes = two
- ifTrue: [res / (p + 1)]
- ifFalse: [p = 0 ifTrue: [res] ifFalse: [res - ((two asLargeFloatPrecision: prec) ln * p)]].
- self < one ifTrue: [res inPlaceNegated].
- ^res asLargeFloatPrecision: precision
+ "Answer the neperian logarithm of the receiver."
+
+ | x4 one two p res selfHighRes prec e |
+ self <= self zero ifTrue: [DomainError signal: 'ln is only defined for x > 0.0'].
+
+ one := self one.
+ self = one ifTrue: [^self zero].
+ two := one timesTwoPower: 1.
+
+ "Use Salamin algorithm (approximation is good if x is big enough)
+ x ln = Pi / (2 * (1 agm: 4/x) ).
+ If x not big enough, compute (x timesTwoPower: p) ln - (2 ln * p)
+ if x is close to 1, better use a power expansion"
+ prec := precision + 16.
+ e := self exponent.
+ e < 0 ifTrue: [e := -1 - e].
+ e > prec
+ ifTrue: [p := 0]
+ ifFalse:
+ [p := prec - e.
+ prec := prec + p highBit].
+ selfHighRes := self asLargeFloatPrecision: prec.
+ (selfHighRes - one) exponent * -4 >= precision ifTrue: [^(selfHighRes powerExpansionLnPrecision: prec) asLargeFloatPrecision: precision].
+ self < one ifTrue: [selfHighRes inPlaceReciprocal]. "Use ln(1/x) => - ln(x)"
+ x4 := (4 asLargeFloatPrecision: prec)
+ inPlaceDivideBy: selfHighRes;
+ inPlaceTimesTwoPower: p negated.
+ res := selfHighRes pi / (one agm: x4) timesTwoPower: -1.
+ res := selfHighRes = two
+ ifTrue: [res / (p + 1)]
+ ifFalse: [p = 0 ifTrue: [res] ifFalse: [res - ((two asLargeFloatPrecision: prec) ln * p)]].
+ self < one ifTrue: [res inPlaceNegated].
+ ^res asLargeFloatPrecision: precision
+
+ "Modified (comment): / 28-05-2019 / 17:43:30 / Claus Gittinger"
!
sin
@@ -1403,25 +1419,27 @@
!
sincos
- "Evaluate the sine and cosine of the receiver"
-
- | pi halfPi quarterPi x sincos sinneg cosneg |
- x _ self moduloNegPiToPi.
- sinneg _ x negative.
- x inPlaceAbs.
- pi _ self piDoublePrecision.
- halfPi _ pi timesTwoPower: -1.
- (cosneg _ x > halfPi) ifTrue: [x inPlaceSubtract: pi; inPlaceNegated].
- quarterPi _ halfPi timesTwoPower: -1.
- x > quarterPi
- ifTrue:
- [x inPlaceSubtract: halfPi; inPlaceNegated.
- sincos _ (self sincos: x) reversed]
- ifFalse:
- [sincos _ self sincos: x].
- sinneg ifTrue: [sincos first inPlaceNegated].
- cosneg ifTrue: [sincos last inPlaceNegated].
- ^sincos collect: [:e| e asLargeFloatPrecision: precision]
+ "Evaluate the sine and cosine of the receiver"
+
+ | pi halfPi quarterPi x sincos sinneg cosneg |
+ x := self moduloNegPiToPi.
+ sinneg := x negative.
+ x inPlaceAbs.
+ pi := self piDoublePrecision.
+ halfPi := pi timesTwoPower: -1.
+ (cosneg := x > halfPi) ifTrue: [x inPlaceSubtract: pi; inPlaceNegated].
+ quarterPi := halfPi timesTwoPower: -1.
+ x > quarterPi
+ ifTrue:
+ [x inPlaceSubtract: halfPi; inPlaceNegated.
+ sincos := (self sincos: x) reversed]
+ ifFalse:
+ [sincos := self sincos: x].
+ sinneg ifTrue: [sincos first inPlaceNegated].
+ cosneg ifTrue: [sincos last inPlaceNegated].
+ ^sincos collect: [:e| e asLargeFloatPrecision: precision]
+
+ "Modified (format): / 28-05-2019 / 17:45:28 / Claus Gittinger"
!
sinh
@@ -1450,24 +1468,24 @@
self isZero ifTrue: [^self].
"use additional bits"
- decimalPlaces _ precision + 16.
- n _ self asLargeFloatPrecision: decimalPlaces.
+ decimalPlaces := precision + 16.
+ n := self asLargeFloatPrecision: decimalPlaces.
"constants"
- one _ n one.
+ one := n one.
"normalize n"
- norm _ n exponent quo: 2.
- n _ n timesTwoPower: norm * -2.
+ norm := n exponent quo: 2.
+ n := n timesTwoPower: norm * -2.
"Initial guess for sqrt(1/n)"
- previousGuess _ self class
+ previousGuess := self class
mantissa: 3
exponent: -2 - (n exponent quo: 2)
precision: decimalPlaces.
- guess _ previousGuess copy.
-
- "use iterations x(k+1) _ x*( 1 + (1-x*x*n)/2) to guess sqrt(1/n)"
+ guess := previousGuess copy.
+
+ "use iterations x(k+1) := x*( 1 + (1-x*x*n)/2) to guess sqrt(1/n)"
[guess inPlaceMultiplyNoRoundBy: guess.
guess inPlaceMultiplyBy: n.
@@ -1475,7 +1493,7 @@
guess inPlaceAddNoRound: one.
"stop when no evolution of precision + 12 first bits"
- stopIteration _ guess isZero or: [guess exponent < (decimalPlaces - 4) negated].
+ stopIteration := guess isZero or: [guess exponent < (decimalPlaces - 4) negated].
guess inPlaceTimesTwoPower: -1.
guess inPlaceAddNoRound: one.
guess inPlaceMultiplyNoRoundBy: previousGuess.
@@ -1492,6 +1510,7 @@
^guess asLargeFloatPrecision: precision
"Modified: / 28-05-2019 / 11:22:33 / Claus Gittinger"
+ "Modified (comment): / 28-05-2019 / 17:45:39 / Claus Gittinger"
!
tan
@@ -1767,40 +1786,46 @@
!LargeFloat methodsFor:'private'!
cos: x
- "Evaluate the cosine of x by recursive cos(2x) formula and power series expansion.
- Note that it is better to use this method with x <= pi/4."
-
- | one cos fraction power |
- x isZero ifTrue: [^x one].
- power _ ((precision bitShift: -1) + x exponent max: 0) highBit.
- fraction _ x timesTwoPower: power negated.
- cos _ fraction powerExpansionCosPrecision: precision + (1 bitShift: 1 + power).
- one _ x one.
- power timesRepeat:
- ["Evaluate cos(2x)=2 cos(x)^2-1"
- cos inPlaceMultiplyBy: cos; inPlaceTimesTwoPower: 1; inPlaceSubtract: one].
- ^cos
+ "Evaluate the cosine of x by recursive cos(2x) formula and power series expansion.
+ Note that it is better to use this method with x <= pi/4."
+
+ | one cos fraction power |
+ x isZero ifTrue: [^x one].
+ power := ((precision bitShift: -1) + x exponent max: 0) highBit.
+ fraction := x timesTwoPower: power negated.
+ cos := fraction powerExpansionCosPrecision: precision + (1 bitShift: 1 + power).
+ one := x one.
+ power timesRepeat:
+ ["Evaluate cos(2x)=2 cos(x)^2-1"
+ cos inPlaceMultiplyBy: cos; inPlaceTimesTwoPower: 1; inPlaceSubtract: one].
+ ^cos
+
+ "Modified (comment): / 28-05-2019 / 17:42:54 / Claus Gittinger"
!
digitCompare: b
- "both are positive or negative.
- answer +1 if i am of greater magnitude, -1 if i am of smaller magnitude, 0 if equal magnitude"
-
- | compare |
- self isZero
- ifTrue: [b isZero
- ifTrue: [^ 0]
- ifFalse: [^ -1]].
- b isZero
- ifTrue: [^ 1].
- compare _ (self exponent - b exponent) sign.
- ^ compare = 0
- ifTrue: [(self abs - b abs) sign]
- ifFalse: [compare]
+ "both are positive or negative.
+ answer +1 if i am of greater magnitude, -1 if i am of smaller magnitude, 0 if equal magnitude"
+
+ | compare |
+ self isZero
+ ifTrue: [b isZero
+ ifTrue: [^ 0]
+ ifFalse: [^ -1]].
+ b isZero
+ ifTrue: [^ 1].
+ compare := (self exponent - b exponent) sign.
+ ^ compare = 0
+ ifTrue: [(self abs - b abs) sign]
+ ifFalse: [compare]
+
+ "Modified (comment): / 28-05-2019 / 17:42:59 / Claus Gittinger"
!
inPlaceAbs
- mantissa _ mantissa abs
+ mantissa := mantissa abs
+
+ "Modified (format): / 28-05-2019 / 17:43:09 / Claus Gittinger"
!
inPlaceAdd: b
@@ -1808,26 +1833,26 @@
b isZero ifTrue: [^self roundToPrecision].
self isZero ifTrue:[
- mantissa _ b mantissa.
- biasedExponent _ b biasedExponent
+ mantissa := b mantissa.
+ biasedExponent := b biasedExponent
] ifFalse:[
biasedExponent = b biasedExponent ifTrue: [
- mantissa _ mantissa + b mantissa
+ mantissa := mantissa + b mantissa
] ifFalse:[
"check for early truncation. beware, keep 2 bits for rounding"
- delta _ biasedExponent - b biasedExponent.
+ delta := biasedExponent - b biasedExponent.
delta - 2 > (precision max: self precisionInMantissa) ifFalse:[
delta negated - 2 > (precision max: b precisionInMantissa) ifTrue:[
- mantissa _ b mantissa.
- biasedExponent _ b exponent
+ mantissa := b mantissa.
+ biasedExponent := b exponent
] ifFalse:[
- delta _ biasedExponent - b exponent.
+ delta := biasedExponent - b exponent.
delta > 0 ifTrue:[
- mantissa _ (self shift: mantissa by: delta) + b mantissa.
- biasedExponent _ biasedExponent - delta
+ mantissa := (self shift: mantissa by: delta) + b mantissa.
+ biasedExponent := biasedExponent - delta
] ifFalse: [
- mantissa _ mantissa + (self shift: b mantissa by: delta negated)
+ mantissa := mantissa + (self shift: b mantissa by: delta negated)
]
]
]
@@ -1837,6 +1862,7 @@
"Created: / 26-05-2019 / 03:44:50 / Claus Gittinger"
"Modified: / 28-05-2019 / 08:49:15 / Claus Gittinger"
+ "Modified (format): / 28-05-2019 / 17:43:15 / Claus Gittinger"
!
inPlaceAddNoRound: b
@@ -2009,31 +2035,32 @@
b isZero ifTrue: [^self roundToPrecision].
self isZero
ifTrue:
- [mantissa _ b mantissa negated.
- biasedExponent _ b biasedExponent]
+ [mantissa := b mantissa negated.
+ biasedExponent := b biasedExponent]
ifFalse:
[biasedExponent = b biasedExponent
- ifTrue: [mantissa _ mantissa - b mantissa]
+ ifTrue: [mantissa := mantissa - b mantissa]
ifFalse:
["check for early truncation. beware, keep 2 bits for rounding"
- delta _ biasedExponent - b biasedExponent.
+ delta := biasedExponent - b biasedExponent.
delta - 2 > (precision max: self precisionInMantissa)
ifFalse:
[delta negated - 2 > (precision max: b precisionInMantissa)
ifTrue:
- [mantissa _ b mantissa negated.
- biasedExponent _ b biasedExponent]
+ [mantissa := b mantissa negated.
+ biasedExponent := b biasedExponent]
ifFalse:
- [delta _ biasedExponent - b biasedExponent.
+ [delta := biasedExponent - b biasedExponent.
delta >= 0
ifTrue:
- [mantissa _ (self shift: mantissa by: delta) - b mantissa.
- biasedExponent _ biasedExponent - delta]
- ifFalse: [mantissa _ mantissa - (self shift: b mantissa by: delta negated)]]]]].
+ [mantissa := (self shift: mantissa by: delta) - b mantissa.
+ biasedExponent := biasedExponent - delta]
+ ifFalse: [mantissa := mantissa - (self shift: b mantissa by: delta negated)]]]]].
self roundToPrecision
"Modified: / 28-05-2019 / 08:48:46 / Claus Gittinger"
+ "Modified (format): / 28-05-2019 / 17:43:21 / Claus Gittinger"
!
inPlaceSubtractNoRound: b
@@ -2058,8 +2085,10 @@
!
inPlaceTimesTwoPower: n
- self isZero
- ifFalse: [biasedExponent _ biasedExponent + n]
+ self isZero
+ ifFalse: [biasedExponent := biasedExponent + n]
+
+ "Modified (format): / 28-05-2019 / 17:43:24 / Claus Gittinger"
!
mantissa:mantissaArg exponent:exponentArg
@@ -2089,23 +2118,25 @@
!
moduloNegPiToPi
- "answer a copy of the receiver modulo 2*pi, with doubled precision"
-
- | x pi twoPi quo |
- x _ (self asLargeFloatPrecision: precision * 2).
- self negative ifTrue: [x inPlaceNegated].
- pi _ x pi.
- twoPi _ pi timesTwoPower: 1.
- x > pi ifTrue:
- [quo _ x + pi quo: twoPi.
- quo highBitOfMagnitude > precision ifTrue:
- [x _ (self abs asLargeFloatPrecision: precision + quo highBitOfMagnitude).
- pi _ x pi.
- twoPi _ pi timesTwoPower: 1.
- quo _ x + pi quo: twoPi].
- x inPlaceSubtract: twoPi * quo.
- self negative ifTrue: [x inPlaceNegated]].
- ^x asLargeFloatPrecision: precision * 2
+ "answer a copy of the receiver modulo 2*pi, with doubled precision"
+
+ | x pi twoPi quo |
+ x := (self asLargeFloatPrecision: precision * 2).
+ self negative ifTrue: [x inPlaceNegated].
+ pi := x pi.
+ twoPi := pi timesTwoPower: 1.
+ x > pi ifTrue:
+ [quo := x + pi quo: twoPi.
+ quo highBitOfMagnitude > precision ifTrue:
+ [x := (self abs asLargeFloatPrecision: precision + quo highBitOfMagnitude).
+ pi := x pi.
+ twoPi := pi timesTwoPower: 1.
+ quo := x + pi quo: twoPi].
+ x inPlaceSubtract: twoPi * quo.
+ self negative ifTrue: [x inPlaceNegated]].
+ ^x asLargeFloatPrecision: precision * 2
+
+ "Modified (format): / 28-05-2019 / 17:43:36 / Claus Gittinger"
!
normalize
@@ -2136,298 +2167,324 @@
!
powerExpansionArCoshp1Precision: precBits
- "Evaluate arcosh(x+1)/sqrt(2*x) for the receiver x by power series expansion.
- The algorithm is interesting when the receiver is close to zero"
-
- | one two count count2 sum term term1 term2 |
- one _ self one.
- two _ one timesTwoPower: 1.
- count _ one copy.
- count2 _ one copy.
- sum _ one copy.
- term1 _ one copy.
- term2 _ one copy.
-
- [term1 inPlaceMultiplyBy: self.
- term1 inPlaceNegated.
- term2 inPlaceMultiplyBy: count2.
- term2 inPlaceMultiplyBy: count2.
- term2 inPlaceDivideBy: count.
- count inPlaceAdd: one.
- count2 inPlaceAdd: two.
- term2 inPlaceDivideBy: count2.
- term2 inPlaceTimesTwoPower: -2.
- term _ term1 * term2.
- sum inPlaceAdd: term.
- term exponent + precBits < sum exponent] whileFalse.
- ^sum
+ "Evaluate arcosh(x+1)/sqrt(2*x) for the receiver x by power series expansion.
+ The algorithm is interesting when the receiver is close to zero"
+
+ | one two count count2 sum term term1 term2 |
+ one := self one.
+ two := one timesTwoPower: 1.
+ count := one copy.
+ count2 := one copy.
+ sum := one copy.
+ term1 := one copy.
+ term2 := one copy.
+
+ [term1 inPlaceMultiplyBy: self.
+ term1 inPlaceNegated.
+ term2 inPlaceMultiplyBy: count2.
+ term2 inPlaceMultiplyBy: count2.
+ term2 inPlaceDivideBy: count.
+ count inPlaceAdd: one.
+ count2 inPlaceAdd: two.
+ term2 inPlaceDivideBy: count2.
+ term2 inPlaceTimesTwoPower: -2.
+ term := term1 * term2.
+ sum inPlaceAdd: term.
+ term exponent + precBits < sum exponent] whileFalse.
+ ^sum
+
+ "Modified (comment): / 28-05-2019 / 17:43:51 / Claus Gittinger"
!
powerExpansionArSinhPrecision: precBits
- "Evaluate the area hypebolic sine of the receiver by power series expansion.
- The algorithm is interesting when the receiver is close to zero"
-
- | one x2 two count sum term |
- one _ self one.
- two _ one timesTwoPower: 1.
- count _ one copy.
- sum _ one copy.
- term _ one copy.
- x2 _ self squared.
-
- [term inPlaceMultiplyBy: x2.
- term inPlaceMultiplyBy: count.
- term inPlaceDivideBy: count + one.
- term inPlaceNegated.
- count inPlaceAdd: two.
- sum inPlaceAdd: term / count.
- term exponent + precBits < sum exponent] whileFalse.
- sum inPlaceMultiplyBy: self.
- ^sum
+ "Evaluate the area hypebolic sine of the receiver by power series expansion.
+ The algorithm is interesting when the receiver is close to zero"
+
+ | one x2 two count sum term |
+ one := self one.
+ two := one timesTwoPower: 1.
+ count := one copy.
+ sum := one copy.
+ term := one copy.
+ x2 := self squared.
+
+ [term inPlaceMultiplyBy: x2.
+ term inPlaceMultiplyBy: count.
+ term inPlaceDivideBy: count + one.
+ term inPlaceNegated.
+ count inPlaceAdd: two.
+ sum inPlaceAdd: term / count.
+ term exponent + precBits < sum exponent] whileFalse.
+ sum inPlaceMultiplyBy: self.
+ ^sum
+
+ "Modified (comment): / 28-05-2019 / 17:43:57 / Claus Gittinger"
!
powerExpansionArTanhPrecision: precBits
- "Evaluate the area hyperbolic tangent of the receiver by power series expansion.
- arTanh (x) = x (1 + x^2/3 + x^4/5 + ... ) for -1 < x < 1
- The algorithm is interesting when the receiver is close to zero"
-
- | one x2 two count sum term |
- one _ self one.
- two _ one timesTwoPower: 1.
- count _ one copy.
- sum _ one copy.
- term _ one copy.
- x2 _ self squared.
-
- [term inPlaceMultiplyBy: x2.
- count inPlaceAdd: two.
- sum inPlaceAdd: term / count.
- term exponent + precBits < sum exponent] whileFalse.
- sum inPlaceMultiplyBy: self.
- ^sum
+ "Evaluate the area hyperbolic tangent of the receiver by power series expansion.
+ arTanh (x) = x (1 + x^2/3 + x^4/5 + ... ) for -1 < x < 1
+ The algorithm is interesting when the receiver is close to zero"
+
+ | one x2 two count sum term |
+ one := self one.
+ two := one timesTwoPower: 1.
+ count := one copy.
+ sum := one copy.
+ term := one copy.
+ x2 := self squared.
+
+ [term inPlaceMultiplyBy: x2.
+ count inPlaceAdd: two.
+ sum inPlaceAdd: term / count.
+ term exponent + precBits < sum exponent] whileFalse.
+ sum inPlaceMultiplyBy: self.
+ ^sum
+
+ "Modified (comment): / 28-05-2019 / 17:44:04 / Claus Gittinger"
!
powerExpansionArcSinPrecision: precBits
- "Evaluate the arc sine of the receiver by power series expansion.
- The algorithm is interesting when the receiver is close to zero"
-
- | one x2 two count sum term |
- one _ self one.
- two _ one timesTwoPower: 1.
- count _ one copy.
- sum _ one copy.
- term _ one copy.
- x2 _ self squared.
-
- [term inPlaceMultiplyBy: x2.
- term inPlaceMultiplyBy: count.
- term inPlaceDivideBy: count + one.
- count inPlaceAdd: two.
- sum inPlaceAdd: term / count.
- term exponent + precBits < sum exponent] whileFalse.
- sum inPlaceMultiplyBy: self.
- ^sum
+ "Evaluate the arc sine of the receiver by power series expansion.
+ The algorithm is interesting when the receiver is close to zero"
+
+ | one x2 two count sum term |
+ one := self one.
+ two := one timesTwoPower: 1.
+ count := one copy.
+ sum := one copy.
+ term := one copy.
+ x2 := self squared.
+
+ [term inPlaceMultiplyBy: x2.
+ term inPlaceMultiplyBy: count.
+ term inPlaceDivideBy: count + one.
+ count inPlaceAdd: two.
+ sum inPlaceAdd: term / count.
+ term exponent + precBits < sum exponent] whileFalse.
+ sum inPlaceMultiplyBy: self.
+ ^sum
+
+ "Modified (comment): / 28-05-2019 / 17:44:09 / Claus Gittinger"
!
powerExpansionArcTanPrecision: precBits
- "Evaluate the arc tangent of the receiver by power series expansion.
- arcTan (x) = x (1 - x^2/3 + x^4/5 - ... ) for -1 < x < 1
- The algorithm is interesting when the receiver is close to zero"
-
- | count one sum term two x2 |
- one _ self one.
- two _ one timesTwoPower: 1.
- count _ one copy.
- sum _ one copy.
- term _ one copy.
- x2 _ self squared.
-
- [term inPlaceMultiplyBy: x2.
- term inPlaceNegated.
- count inPlaceAdd: two.
- sum inPlaceAdd: term / count.
- term exponent + precBits < sum exponent] whileFalse.
- sum inPlaceMultiplyBy: self.
- ^sum
+ "Evaluate the arc tangent of the receiver by power series expansion.
+ arcTan (x) = x (1 - x^2/3 + x^4/5 - ... ) for -1 < x < 1
+ The algorithm is interesting when the receiver is close to zero"
+
+ | count one sum term two x2 |
+ one := self one.
+ two := one timesTwoPower: 1.
+ count := one copy.
+ sum := one copy.
+ term := one copy.
+ x2 := self squared.
+
+ [term inPlaceMultiplyBy: x2.
+ term inPlaceNegated.
+ count inPlaceAdd: two.
+ sum inPlaceAdd: term / count.
+ term exponent + precBits < sum exponent] whileFalse.
+ sum inPlaceMultiplyBy: self.
+ ^sum
+
+ "Modified (comment): / 28-05-2019 / 17:44:15 / Claus Gittinger"
!
powerExpansionCosPrecision: precBits
- "Evaluate the cosine of the receiver by power series expansion.
- The algorithm is interesting when the receiver is close to zero"
-
- | count one sum term two x2 |
- one _ self one.
- two _ one timesTwoPower: 1.
- count _ one copy.
- sum _ one copy.
- term _ one copy.
- x2 _ self squared.
-
- [term inPlaceMultiplyBy: x2.
- term inPlaceDivideBy: count * (count + one).
- term inPlaceNegated.
- count inPlaceAdd: two.
- sum inPlaceAdd: term.
- term exponent + precBits < sum exponent] whileFalse.
- ^sum
+ "Evaluate the cosine of the receiver by power series expansion.
+ The algorithm is interesting when the receiver is close to zero"
+
+ | count one sum term two x2 |
+ one := self one.
+ two := one timesTwoPower: 1.
+ count := one copy.
+ sum := one copy.
+ term := one copy.
+ x2 := self squared.
+
+ [term inPlaceMultiplyBy: x2.
+ term inPlaceDivideBy: count * (count + one).
+ term inPlaceNegated.
+ count inPlaceAdd: two.
+ sum inPlaceAdd: term.
+ term exponent + precBits < sum exponent] whileFalse.
+ ^sum
+
+ "Modified (comment): / 28-05-2019 / 17:44:25 / Claus Gittinger"
!
powerExpansionCoshPrecision: precBits
- "Evaluate the hyperbolic cosine of the receiver by power series expansion.
- The algorithm is interesting when the receiver is close to zero"
-
- | count one sum term two x2 |
- one _ self one.
- two _ one timesTwoPower: 1.
- count _ one copy.
- sum _ one copy.
- term _ one copy.
- x2 _ self squared.
-
- [term inPlaceMultiplyBy: x2.
- term inPlaceDivideBy: count * (count + one).
- count inPlaceAdd: two.
- sum inPlaceAdd: term.
- term exponent + precBits < sum exponent] whileFalse.
- ^sum
+ "Evaluate the hyperbolic cosine of the receiver by power series expansion.
+ The algorithm is interesting when the receiver is close to zero"
+
+ | count one sum term two x2 |
+ one := self one.
+ two := one timesTwoPower: 1.
+ count := one copy.
+ sum := one copy.
+ term := one copy.
+ x2 := self squared.
+
+ [term inPlaceMultiplyBy: x2.
+ term inPlaceDivideBy: count * (count + one).
+ count inPlaceAdd: two.
+ sum inPlaceAdd: term.
+ term exponent + precBits < sum exponent] whileFalse.
+ ^sum
+
+ "Modified (comment): / 28-05-2019 / 17:44:30 / Claus Gittinger"
!
powerExpansionLnPrecision: precBits
- "Evaluate the neperian logarithm of the receiver by power series expansion.
- For quadratic convergence, use:
- ln ((1+y)/(1-y)) = 2 y (1 + y^2/3 + y^4/5 + ... ) = 2 ar tanh( y )
- (1+y)/(1-y) = self => y = (self-1)/(self+1)
- This algorithm is interesting when the receiver is close to 1"
-
- | one |
- one _ self one.
- ^((self - one)/(self + one) powerExpansionArTanhPrecision: precBits) timesTwoPower: 1
+ "Evaluate the neperian logarithm of the receiver by power series expansion.
+ For quadratic convergence, use:
+ ln ((1+y)/(1-y)) = 2 y (1 + y^2/3 + y^4/5 + ... ) = 2 ar tanh( y )
+ (1+y)/(1-y) = self => y = (self-1)/(self+1)
+ This algorithm is interesting when the receiver is close to 1"
+
+ | one |
+ one := self one.
+ ^((self - one)/(self + one) powerExpansionArTanhPrecision: precBits) timesTwoPower: 1
+
+ "Modified (comment): / 28-05-2019 / 17:44:38 / Claus Gittinger"
!
powerExpansionSinCosPrecision: precBits
- "Evaluate the sine and cosine of the receiver by power series expansion.
- The algorithm is interesting when the receiver is close to zero"
-
- | count one sin cos term |
- one _ self one.
- count _ one copy.
- cos _ one copy.
- sin _ self copy.
- term _ self copy.
-
- [count inPlaceAdd: one.
- term
- inPlaceMultiplyBy: self;
- inPlaceDivideBy: count;
- inPlaceNegated.
- cos inPlaceAdd: term.
-
- count inPlaceAdd: one.
- term
- inPlaceMultiplyBy: self;
- inPlaceDivideBy: count.
- sin inPlaceAdd: term.
-
- term exponent + precBits < sin exponent] whileFalse.
- ^Array with: sin with: cos
+ "Evaluate the sine and cosine of the receiver by power series expansion.
+ The algorithm is interesting when the receiver is close to zero"
+
+ | count one sin cos term |
+ one := self one.
+ count := one copy.
+ cos := one copy.
+ sin := self copy.
+ term := self copy.
+
+ [count inPlaceAdd: one.
+ term
+ inPlaceMultiplyBy: self;
+ inPlaceDivideBy: count;
+ inPlaceNegated.
+ cos inPlaceAdd: term.
+
+ count inPlaceAdd: one.
+ term
+ inPlaceMultiplyBy: self;
+ inPlaceDivideBy: count.
+ sin inPlaceAdd: term.
+
+ term exponent + precBits < sin exponent] whileFalse.
+ ^Array with: sin with: cos
+
+ "Modified (comment): / 28-05-2019 / 17:44:43 / Claus Gittinger"
!
powerExpansionSinPrecision: precBits
- "Evaluate the sine of the receiver by power series expansion.
- The algorithm is interesting when the receiver is close to zero"
-
- | count one sum term two x2 |
- one _ self one.
- two _ one timesTwoPower: 1.
- count _ two copy.
- sum _ self copy.
- term _ self copy.
- x2 _ self squared.
-
- [term inPlaceMultiplyBy: x2.
- term inPlaceDivideBy: count * (count + one).
- term inPlaceNegated.
- count inPlaceAdd: two.
- sum inPlaceAdd: term.
- term exponent + precBits < sum exponent] whileFalse.
- ^sum
+ "Evaluate the sine of the receiver by power series expansion.
+ The algorithm is interesting when the receiver is close to zero"
+
+ | count one sum term two x2 |
+ one := self one.
+ two := one timesTwoPower: 1.
+ count := two copy.
+ sum := self copy.
+ term := self copy.
+ x2 := self squared.
+
+ [term inPlaceMultiplyBy: x2.
+ term inPlaceDivideBy: count * (count + one).
+ term inPlaceNegated.
+ count inPlaceAdd: two.
+ sum inPlaceAdd: term.
+ term exponent + precBits < sum exponent] whileFalse.
+ ^sum
+
+ "Modified (comment): / 28-05-2019 / 17:44:50 / Claus Gittinger"
!
powerExpansionSinhPrecision: precBits
- "Evaluate the hyperbolic sine of the receiver by power series expansion.
- The algorithm is interesting when the receiver is close to zero"
-
- | count one sum term two x2 |
- one _ self one.
- two _ one timesTwoPower: 1.
- count _ two copy.
- sum _ self copy.
- term _ self copy.
- x2 _ self squared.
-
- [term inPlaceMultiplyBy: x2.
- term inPlaceDivideBy: count * (count + one).
- count inPlaceAdd: two.
- sum inPlaceAdd: term.
- term exponent + precBits < sum exponent] whileFalse.
- ^sum
+ "Evaluate the hyperbolic sine of the receiver by power series expansion.
+ The algorithm is interesting when the receiver is close to zero"
+
+ | count one sum term two x2 |
+ one := self one.
+ two := one timesTwoPower: 1.
+ count := two copy.
+ sum := self copy.
+ term := self copy.
+ x2 := self squared.
+
+ [term inPlaceMultiplyBy: x2.
+ term inPlaceDivideBy: count * (count + one).
+ count inPlaceAdd: two.
+ sum inPlaceAdd: term.
+ term exponent + precBits < sum exponent] whileFalse.
+ ^sum
+
+ "Modified (comment): / 28-05-2019 / 17:44:56 / Claus Gittinger"
!
powerExpansionTanPrecision: precBits
- "Evaluate the tangent of the receiver by power series expansion.
- The algorithm is interesting when the receiver is close to zero"
-
- | count one sum term pow two x2 seidel |
- one _ self one.
- two _ one timesTwoPower: 1.
- count _ two copy.
- sum _ one copy.
- pow _ one copy.
- x2 _ self squared.
- seidel _ OrderedCollection new: 256.
- seidel add: 1.
-
- [pow inPlaceMultiplyBy: x2.
- pow inPlaceDivideBy: count * (count + one).
- count inPlaceAdd: two.
- 2 to: seidel size do: [:i | seidel at: i put: (seidel at: i-1) + (seidel at: i)].
- seidel addLast: seidel last.
- seidel size to: 2 by: -1 do: [:i | seidel at: i - 1 put: (seidel at: i-1) + (seidel at: i)].
- seidel addFirst: seidel first.
- term _ pow * seidel first.
- sum inPlaceAdd: term.
- term exponent + precBits < sum exponent] whileFalse.
- sum inPlaceMultiplyBy: self.
- ^sum
+ "Evaluate the tangent of the receiver by power series expansion.
+ The algorithm is interesting when the receiver is close to zero"
+
+ | count one sum term pow two x2 seidel |
+ one := self one.
+ two := one timesTwoPower: 1.
+ count := two copy.
+ sum := one copy.
+ pow := one copy.
+ x2 := self squared.
+ seidel := OrderedCollection new: 256.
+ seidel add: 1.
+
+ [pow inPlaceMultiplyBy: x2.
+ pow inPlaceDivideBy: count * (count + one).
+ count inPlaceAdd: two.
+ 2 to: seidel size do: [:i | seidel at: i put: (seidel at: i-1) + (seidel at: i)].
+ seidel addLast: seidel last.
+ seidel size to: 2 by: -1 do: [:i | seidel at: i - 1 put: (seidel at: i-1) + (seidel at: i)].
+ seidel addFirst: seidel first.
+ term := pow * seidel first.
+ sum inPlaceAdd: term.
+ term exponent + precBits < sum exponent] whileFalse.
+ sum inPlaceMultiplyBy: self.
+ ^sum
+
+ "Modified (comment): / 28-05-2019 / 17:45:02 / Claus Gittinger"
!
powerExpansionTanhPrecision: precBits
- "Evaluate the hyperbolic tangent of the receiver by power series expansion.
- The algorithm is interesting when the receiver is close to zero"
-
- | count one sum term pow two x2 seidel |
- one _ self one.
- two _ one timesTwoPower: 1.
- count _ two copy.
- sum _ one copy.
- pow _ one copy.
- x2 _ self squared.
- seidel _ OrderedCollection new: 256.
- seidel add: 1.
-
- [pow inPlaceMultiplyBy: x2.
- pow inPlaceDivideBy: count * (count + one).
- pow inPlaceNegated.
- count inPlaceAdd: two.
- 2 to: seidel size do: [:i | seidel at: i put: (seidel at: i-1) + (seidel at: i)].
- seidel addLast: seidel last.
- seidel size to: 2 by: -1 do: [:i | seidel at: i - 1 put: (seidel at: i-1) + (seidel at: i)].
- seidel addFirst: seidel first.
- term _ pow * seidel first.
- sum inPlaceAdd: term.
- term exponent + precBits < sum exponent] whileFalse.
- sum inPlaceMultiplyBy: self.
- ^sum
+ "Evaluate the hyperbolic tangent of the receiver by power series expansion.
+ The algorithm is interesting when the receiver is close to zero"
+
+ | count one sum term pow two x2 seidel |
+ one := self one.
+ two := one timesTwoPower: 1.
+ count := two copy.
+ sum := one copy.
+ pow := one copy.
+ x2 := self squared.
+ seidel := OrderedCollection new: 256.
+ seidel add: 1.
+
+ [pow inPlaceMultiplyBy: x2.
+ pow inPlaceDivideBy: count * (count + one).
+ pow inPlaceNegated.
+ count inPlaceAdd: two.
+ 2 to: seidel size do: [:i | seidel at: i put: (seidel at: i-1) + (seidel at: i)].
+ seidel addLast: seidel last.
+ seidel size to: 2 by: -1 do: [:i | seidel at: i - 1 put: (seidel at: i-1) + (seidel at: i)].
+ seidel addFirst: seidel first.
+ term := pow * seidel first.
+ sum inPlaceAdd: term.
+ term exponent + precBits < sum exponent] whileFalse.
+ sum inPlaceMultiplyBy: self.
+ ^sum
+
+ "Modified (comment): / 28-05-2019 / 17:45:07 / Claus Gittinger"
!
precision:precisionArg
@@ -2451,15 +2508,17 @@
!
reduce
- "remove trailing zero bits from mantissa so that we can do arithmetic on smaller integer
- (that will un-normalize self)"
-
- | trailing |
- trailing _ mantissa abs lowBit - 1.
- trailing > 0
- ifFalse: [ ^ self ].
- mantissa _ self shift: mantissa by: trailing negated.
- biasedExponent _ biasedExponent + trailing
+ "remove trailing zero bits from mantissa so that we can do arithmetic on smaller integer
+ (that will un-normalize self)"
+
+ | trailing |
+ trailing := mantissa abs lowBit - 1.
+ trailing > 0
+ ifFalse: [ ^ self ].
+ mantissa := self shift: mantissa by: trailing negated.
+ biasedExponent := biasedExponent + trailing
+
+ "Modified (comment): / 28-05-2019 / 17:45:13 / Claus Gittinger"
!
roundToPrecision
@@ -2499,42 +2558,46 @@
!
sin: x
- "Evaluate the sine of x by sin(5x) formula and power series expansion."
-
- | sin sin2 sin4 fifth five |
- x isZero ifTrue: [^x zero].
- five _ 5 asLargeFloatPrecision: x precision.
- fifth _ x / five.
- sin _ fifth powerExpansionSinPrecision: precision + 8.
- sin2 _ sin squared.
- sin2 inPlaceTimesTwoPower: 2.
- sin4 _ sin2 squared.
- sin2 inPlaceMultiplyBy: five.
- ^sin4
- inPlaceSubtract: sin2;
- inPlaceAdd: five;
- inPlaceMultiplyBy: sin;
- yourself
+ "Evaluate the sine of x by sin(5x) formula and power series expansion."
+
+ | sin sin2 sin4 fifth five |
+ x isZero ifTrue: [^x zero].
+ five := 5 asLargeFloatPrecision: x precision.
+ fifth := x / five.
+ sin := fifth powerExpansionSinPrecision: precision + 8.
+ sin2 := sin squared.
+ sin2 inPlaceTimesTwoPower: 2.
+ sin4 := sin2 squared.
+ sin2 inPlaceMultiplyBy: five.
+ ^sin4
+ inPlaceSubtract: sin2;
+ inPlaceAdd: five;
+ inPlaceMultiplyBy: sin;
+ yourself
+
+ "Modified (format): / 28-05-2019 / 17:45:21 / Claus Gittinger"
!
sincos: x
- "Evaluate the sine and cosine of x by recursive sin(2x) and cos(2x) formula and power series expansion.
- Note that it is better to use this method with x <= pi/4."
-
- | one sin cos sincos fraction power |
- x isZero ifTrue: [^Array with: x zero with: x one].
- power _ ((precision bitShift: -1) + x exponent max: 0) highBit.
- fraction _ x timesTwoPower: power negated.
- sincos _ fraction powerExpansionSinCosPrecision: precision + (1 bitShift: 1 + power).
- sin _ sincos first.
- cos _ sincos last.
- one _ x one.
- power timesRepeat:
- ["Evaluate sin(2x)=2 sin(x) cos(x)"
- sin inPlaceMultiplyBy: cos; inPlaceTimesTwoPower: 1.
- "Evaluate cos(2x)=2 cos(x)^2-1"
- cos inPlaceMultiplyBy: cos; inPlaceTimesTwoPower: 1; inPlaceSubtract: one].
- ^sincos
+ "Evaluate the sine and cosine of x by recursive sin(2x) and cos(2x) formula and power series expansion.
+ Note that it is better to use this method with x <= pi/4."
+
+ | one sin cos sincos fraction power |
+ x isZero ifTrue: [^Array with: x zero with: x one].
+ power := ((precision bitShift: -1) + x exponent max: 0) highBit.
+ fraction := x timesTwoPower: power negated.
+ sincos := fraction powerExpansionSinCosPrecision: precision + (1 bitShift: 1 + power).
+ sin := sincos first.
+ cos := sincos last.
+ one := x one.
+ power timesRepeat:
+ ["Evaluate sin(2x)=2 sin(x) cos(x)"
+ sin inPlaceMultiplyBy: cos; inPlaceTimesTwoPower: 1.
+ "Evaluate cos(2x)=2 cos(x)^2-1"
+ cos inPlaceMultiplyBy: cos; inPlaceTimesTwoPower: 1; inPlaceSubtract: one].
+ ^sincos
+
+ "Modified (comment): / 28-05-2019 / 17:45:33 / Claus Gittinger"
! !
!LargeFloat methodsFor:'queries'!