--- a/Complex.st Sat May 05 22:56:23 2018 +0200
+++ b/Complex.st Sat May 05 23:44:08 2018 +0200
@@ -1,3 +1,5 @@
+"{ Encoding: utf8 }"
+
"
This is a Manchester Goodie. It is distributed freely on condition
that you observe these conditions in respect of the whole Goodie, and on
@@ -65,11 +67,59 @@
[Instance variables:]
real <Number> the part of the number which can be expressed as a Real number
imaginary <Number> the part of the number which, in terms of how the number behaves,
- has been multiplied by 'i' (-1 sqrt)
+ has been multiplied by 'i' (-1 sqrt)
+
+ [Constructors:]
+ 5 i
+ 6 + 7 i.
+ 5.6 - 8 i.
+ Complex real: 10 imaginary: 5.
+ Complex abs: 5 arg: (Float pi / 4)
+
+ NOTE:
+ Although Complex something similiar to the Smalltalk''s Number class,
+ it would not be a good idea to make a Complex to be a subclass of a Number because:
+ - Number is subclass of Magnitude and Complex is certainly not a magnitude.
+ Complex does not behave very well as a Magnitude. Operations such as
+ <
+ >
+ <=
+ >=
+ do not make sense in case of complex numbers.
+ - Methods in the following Number methods'' categories do not make sense for a Complex numbers
+ truncation and round off
+ testing
+ intervals
+ comparing
+
+ However the following Number methods'' categories do have sense for a Complex number
+ arithmetic (with the exception of operation
+ //
+ \\
+ quo:
+ rem:
+ mathematical functions
+
+ Thus Complex is somewhat similar to a Number but it is not a subclass of it.
+ Some operations we would like to inherit (e.g. #abs, #negated, #reciprocal)
+ but some of the Number operations do not have sens to inherit or to overload.
+ Classes are not always neat mechanism.
+
+ !!!!!! We had to COPY the implementation of some methods
+ abs
+ negated
+ reciprocal
+ log:
+ isZero
+ reciprocal
+ ...
+ methods from the Number class to the Complex class.
+ Awful solution. Now I begin to appreciate Self.
[Author:]
- Kurt Hebel (hebel@uinova.cerl.uiuc.edu)
- minor changes and double dispatching code by cg.
+ Kurt Hebel (hebel@uinova.cerl.uiuc.edu)
+ minor changes and double dispatching code by cg.
+ additions (trigonometric) and fixes by cg.
"
!
@@ -131,17 +181,20 @@
!Complex class methodsFor:'instance creation'!
abs:aNumber1 arg:aNumber2
- |real imaginary|
+ ^ self
+ real:(aNumber1 * aNumber2 cos)
+ imaginary:(aNumber1 * aNumber2 sin).
- real := aNumber1 * aNumber2 cos.
- imaginary := aNumber1 * aNumber2 sin.
- ^ real + imaginary i
+ "
+ self abs:10 arg:2 -> (-4.16146836547142387+9.092974268256816954i)
+ "
!
fromReal: aNumber
"Create a new complex number from the given real number."
- ^ self basicNew setReal: aNumber setImaginary: 0
+ ^ self basicNew
+ setReal:aNumber setImaginary:0
"
Complex fromReal:1.0
@@ -150,12 +203,12 @@
"Modified (comment): / 12-06-2017 / 20:42:56 / cg"
!
-imaginary: v
+imaginary:v
"Create a new complex number with 0 as real and given imaginary parts.
If the imaginary part is zero, return the real part of the number."
v = 0 ifTrue: [^ 0].
- ^ self basicNew setReal: 0 setImaginary: v
+ ^ self basicNew setReal:0 setImaginary:v
"
Complex imaginary:1.0
@@ -165,10 +218,10 @@
"Modified (comment): / 12-06-2017 / 20:44:51 / cg"
!
-real: aNumber
+real:aNumber
"Create a new complex number from the given real number."
- ^ self basicNew setReal: aNumber setImaginary: 0
+ ^ self basicNew setReal:aNumber setImaginary:0
"
Complex real:1.0
@@ -177,12 +230,12 @@
"Modified (comment): / 12-06-2017 / 20:42:14 / cg"
!
-real: u imaginary: v
- "Create a new complex number with the given real and imaginary parts. If the
- imaginary part is zero, return the real part of the number."
+real:u imaginary:v
+ "Create a new complex number with the given real and imaginary parts.
+ If the imaginary part is zero, return the real part of the number."
v = 0 ifTrue: [^ u].
- ^ self basicNew setReal: u setImaginary: v
+ ^ self basicNew setReal:u setImaginary:v
"
Complex real:1.0 imaginary:2.0
@@ -204,26 +257,27 @@
unity
"Answer the value which allows, for any given arithmetic value, the following to be true:
- aNumber * aNumber class unity = aNumber
+ aNumber * aNumber class unity = aNumber
This must be true regardless of how a given subclass chooses to define #*"
ComplexOne isNil ifTrue:[
- ComplexOne := self fromReal: 1
+ ComplexOne := self basicNew setReal:1 setImaginary:0
].
^ ComplexOne
!
zero
"Answer the value which allows, for any given arithmetic value, the following to be true:
- aNumber + aNumber class zero = aNumber
+ aNumber + aNumber class zero = aNumber
This must be true regardless of how a given subclass chooses to define #+"
ComplexZero isNil ifTrue:[
- ComplexZero := self fromReal: 0
+ ComplexZero := self basicNew setReal:0 setImaginary:0
].
^ ComplexZero
! !
+
!Complex methodsFor:'accessing'!
imaginary
@@ -257,6 +311,10 @@
* aNumber
"Return the product of the receiver and the argument."
+"/ calling productFromComplex:self
+"/ is just as fast than an isComplex call
+"/ (unless the VM cheats)
+
"/ | u v r i |
"/
"/ aNumber isComplex ifTrue:[
@@ -275,6 +333,10 @@
+ aNumber
"Return the sum of the receiver and the argument."
+"/ calling sumFromComplex:self
+"/ is just as fast than an isComplex call
+"/ (unless the VM cheats)
+
"/ | r i |
"/
"/ aNumber isComplex ifTrue: [
@@ -291,6 +353,10 @@
- aNumber
"Return the difference of the receiver and the argument."
+"/ calling differenceFromComplex:self
+"/ is just as fast than an isComplex call
+"/ (unless the VM cheats)
+
"/ | r i |
"/
"/ aNumber isComplex ifTrue: [
@@ -307,6 +373,10 @@
/ aNumber
"Return the quotient of the receiver and the argument."
+"/ calling quotientFromComplex:self
+"/ is just as fast than an isComplex call
+"/ (unless the VM cheats)
+
"/ | denom u v r i |
"/
"/ aNumber isComplex ifTrue:[
@@ -325,12 +395,31 @@
abs
"Return the magnitude (or absolute value) of the complex number
- (that's the distance from the origin in the complex plane)."
+ that's the distance from zero (the origin in the complex plane)."
^ (real * real + (imaginary * imaginary)) sqrt
"
- (1 % 1) abs
+ (1 + 1i) abs
+ "
+!
+
+absSecure
+ "Answer the distance of the receiver from zero (0 + 0 i).
+ Try avoiding overflow and/or underflow"
+
+ | scale |
+
+ scale := real abs max: imaginary abs.
+ scale isZero ifTrue: [^ scale].
+ ^ (self class
+ real: real / scale
+ imaginary: imaginary / scale)
+ squaredNorm sqrt * scale
+
+ "
+ (10 + 4i) abs -> 10.770329614269
+ (10 + 4i) absSecure -> 10.770329614269
"
!
@@ -343,6 +432,73 @@
imaginary: imaginary negated
!
+divideFastAndSecureBy: anObject
+ "Answer the result of dividing receiver by aNumber"
+ " Both operands are scaled to avoid arithmetic overflow.
+ This algorithm works for a wide range of values, and it needs only three divisions.
+ Note: #reciprocal uses #/ for devision "
+
+ | r d newReal newImaginary |
+
+ anObject isComplex ifTrue:[
+ (anObject real abs) > (anObject imaginary abs) ifTrue:[
+ r := anObject imaginary / anObject real.
+ d := r*anObject imaginary + anObject real.
+ newReal := r*imaginary + real/d.
+ newImaginary := r negated * real + imaginary/d.
+ ] ifFalse:[
+ r := anObject real / anObject imaginary.
+ d := r*anObject real + anObject imaginary.
+ newReal := r*real + imaginary/d.
+ newImaginary := r*imaginary - real/d.
+ ].
+
+ ^ Complex real: newReal imaginary: newImaginary
+ ].
+ ^ anObject adaptToComplex: self andSend: #/.
+!
+
+divideSecureBy:anObject
+ "Answer the result of dividing receiver by aNumber"
+
+ "Both operands are scaled to avoid arithmetic overflow.
+ This algorithm works for a wide range of values,
+ but it requires six divisions.
+ #divideFastAndSecureBy: is also quite good,
+ but it uses only 3 divisions.
+ Note: #reciprocal uses #/ for devision"
+
+ | s ars ais brs bis newReal newImaginary |
+
+ anObject isComplex ifTrue:[
+ s := anObject real abs + anObject imaginary abs.
+ ars := self real / s.
+ ais := self imaginary / s.
+ brs := anObject real / s.
+ bis := anObject imaginary / s.
+ s := brs squared + bis squared.
+
+ newReal := ars*brs + (ais*bis) /s.
+ newImaginary := ais*brs - (ars*bis)/s.
+ ^ Complex real:newReal imaginary:newImaginary
+ ].
+ ^ anObject adaptToComplex:self andSend:#/.
+!
+
+i
+ "Answer the result of multiplying the receiver with pure imaginary.
+ ^ self * 1 i
+ This is an obvious extension of method i implemented in Number."
+
+ ^ self class
+ real: imaginary negated
+ imaginary: real
+
+ "
+ (10+4i) i -> (-4+10i)
+ "
+!
+
modulus
| absReal absImag multiplicand quotient |
@@ -367,18 +523,56 @@
imaginary: imaginary negated
! !
-!Complex methodsFor:'coercing'!
+!Complex methodsFor:'coercing & converting'!
+
+asComplex
+ "I am a complex - so return the receiver"
+
+ ^ self
+!
-generality
- ^ 150
-! !
+asFloat
+ imaginary = 0 ifTrue: [^ real asFloat].
+ ^ Number
+ raise: #coercionErrorSignal
+ receiver: self
+ selector: #asFloat
+ errorString: 'Can''t coerce an instance of Complex to a Float'
+!
-!Complex methodsFor:'coercing & converting'!
+asInteger
+ imaginary = 0 ifTrue: [^real asInteger].
+ ^ Number
+ raise: #coercionErrorSignal
+ receiver: self
+ selector: #asInteger
+ errorString: 'Can''t coerce an instance of Complex to an Integer'
+!
+
+asPoint
+ "Return the complex number as a point."
+
+ ^ real @ imaginary
+!
coerce: aNumber
"convert the argument aNumber into an instance of the receiver's class and return it."
^ aNumber asComplex
+!
+
+generality
+ ^ 150
+!
+
+reduceGeneralityIfPossible
+ "Answer the receiver transformed to a lower generality, if such a
+ transformation is possible without losing information.
+ If not, answer the receiver"
+
+ imaginary isZero
+ ifTrue: [^ real]
+ ifFalse: [^ self]
! !
!Complex methodsFor:'comparing'!
@@ -414,7 +608,7 @@
!
hash
- "Hash is implemented because equals is implemented."
+ "Hash is implemented because = is implemented."
^ (real hash) bitXor:(imaginary hash bitShift:16)
@@ -424,48 +618,6 @@
"
! !
-!Complex methodsFor:'converting'!
-
-asComplex
- "I am a complex - so return the receiver"
-
- ^ self
-!
-
-asFloat
- imaginary = 0 ifTrue: [^ real asFloat].
- ^ Number
- raise: #coercionErrorSignal
- receiver: self
- selector: #asFloat
- errorString: 'Can''t coerce an instance of Complex to a Float'
-!
-
-asInteger
- imaginary = 0 ifTrue: [^real asInteger].
- ^ Number
- raise: #coercionErrorSignal
- receiver: self
- selector: #asInteger
- errorString: 'Can''t coerce an instance of Complex to an Integer'
-!
-
-asPoint
- "Return the complex number as a point."
-
- ^ real @ imaginary
-!
-
-reduceGeneralityIfPossible
- "Answer the receiver transformed to a lower generality, if such a
- transformation is possible without losing information.
- If not, answer the receiver"
-
- imaginary isZero
- ifTrue: [^ real]
- ifFalse: [^ self]
-! !
-
!Complex methodsFor:'double dispatching'!
differenceFromComplex:aComplex
@@ -731,10 +883,35 @@
"
!
+arg
+ "Answer the argument of the receiver."
+
+ self isZero ifTrue: [self error: 'zero has no argument.'].
+ ^ imaginary arcTan:real
+!
+
exp
"Return the complex exponential of the receiver."
- ^ (imaginary cos % imaginary sin) * real exp
+ ^ (Complex
+ real:imaginary cos
+ imaginary:imaginary sin) * real exp
+
+ "
+ (10+4i) exp -> (-14397.45885694595768-16669.6842763244522i)
+ "
+!
+
+ln
+ "Answer the natural log of the receiver."
+
+ ^ self abs ln + self arg i
+!
+
+log:base
+ "Answer the log base aNumber of the receiver."
+
+ ^ self ln / base ln
!
sqrt
@@ -782,6 +959,12 @@
"
-4 asComplex sqrt squared
"
+!
+
+squaredNorm
+ "Answer the square of receiver's norm."
+
+ ^real * real + (imaginary * imaginary)
! !
!Complex methodsFor:'printing'!
@@ -854,10 +1037,301 @@
!
sign
- "return a new complex, consisting of the signs of the real and imaginary parts.
- Q: is this a good thing to do ?"
+ imaginary isZero ifTrue:[^ real sign].
+ ^ self / self abs
+
+ "/ squeak does:
+ "/ ^ real sign
+ "/
+ "/ the old code was:
+ "/ return a new complex, consisting of the signs of the real and imaginary parts.
+ "/ Q: is this a good thing to do ?
+ "/ cg: no
+
+ "/ ^ self class real:(real sign) imaginary:(imaginary sign)
+
+ " this is consistent with Wolfram's output:
+
+ (Complex new setReal:1 setImaginary:0) sign -> 1
+ (Complex new setReal:-1 setImaginary:0) sign -> -1
+ (Complex new setReal:0 setImaginary:0) sign -> 0
+
+ (Complex new setReal:1 setImaginary:1) sign -> (0.707106781186547+0.707106781186547i)
+ (Complex new setReal:0 setImaginary:1) sign -> (0.0+1.0i)
+ (Complex new setReal:-1 setImaginary:1) sign -> (-0.707106781186547+0.707106781186547i)
+
+ (Complex new setReal:1 setImaginary:-1) sign -> (0.707106781186547-0.707106781186547i)
+ (Complex new setReal:0 setImaginary:-1) sign -> (0.0-1.0i)
+ (Complex new setReal:-1 setImaginary:-1) sign -> (-0.707106781186547-0.707106781186547i)
+ "
+! !
+
+!Complex methodsFor:'trigonometric functions'!
+
+arcCos
+ "Answer the arc cosine of the receiver.
+ This is the inverse function of cos."
+
+ | x y tmp sh2y shy delta ch2y chy |
+
+ imaginary = 0 ifTrue: [
+ real abs > 1 ifTrue:[
+ x := real < 0 ifTrue: [Float pi]
+ ifFalse: [0].
+ y := real copySignTo: real abs arcCosh.
+ ^ self class real: x imaginary: y
+ ] ifFalse: [
+ ^self class real: real arcCos imaginary: 0
+ ]
+ ].
+ tmp := self squaredNorm - 1 / 2.
+ delta := tmp squared + imaginary squared.
+ sh2y := tmp + delta sqrt.
+ shy := sh2y sqrt.
+ ch2y := 1 + sh2y.
+ chy := ch2y sqrt.
+ y := imaginary copySignTo: shy arSinh.
+ x := (real / chy) arcCos.
+ ^ self class real: x imaginary: y negated
+
+ "
+ (10+4i) cos (-22.91356068209214107-14.84629106966035727i)
+ (10+4i) cos arcCos
+ "
+!
+
+arcCosh
+ "Answer the receiver's area hyperbolic cosine.
+ That is the inverse function of cosh.
+ Some possible implementations:
+ ^imaginary > 0
+ ifTrue: [(self + (self * self - 1) sqrt) ln]
+ ifFalse: [(self + (self * self - 1) sqrt) ln negated]
+ ^self arcCos i
+ This implementation provides an answer with a positive real part.
+ It also avoids creating intermediate Complex."
+
+ | x y tmp sh2x shx delta ch2x chx |
+
+ imaginary = 0 ifTrue: [
+ real abs > 1 ifTrue:[
+ y := real < 0 ifTrue: [Float pi]
+ ifFalse: [0].
+ x := real abs arcCosh.
+ ^ self class real: x imaginary: y
+ ] ifFalse: [
+ ^ self class real: 0 imaginary: real arcCos
+ ]
+ ].
+ tmp := self squaredNorm - 1 / 2.
+ delta := tmp squared + imaginary squared.
+ sh2x := tmp + delta sqrt.
+ shx := sh2x sqrt.
+ ch2x := 1 + sh2x.
+ chx := ch2x sqrt.
+ x := shx arSinh.
+ y := imaginary copySignTo: (real / chx) arcCos.
+ ^ self class real: x imaginary: y
+
+ "
+ (10 + 4i) arcCosh
+ "
+!
+
+arcSin
+ "Answer the arc sine of the receiver.
+ This is the inverse function of sin."
+
+ | x y tmp delta sh2y shy ch2y chy |
+
+ imaginary = 0 ifTrue:[
+ real abs > 1 ifTrue:[
+ x := Float pi / 2 * real sign.
+ y := (real copySignTo: real abs arcCosh) negated.
+ ^ self class real: x imaginary: y
+ ] ifFalse: [
+ ^ self class real: real arcSin imaginary: 0
+ ]
+ ].
+ tmp := (self squaredNorm - 1) / 2.
+ delta := tmp squared + imaginary squared.
+ sh2y := tmp + delta sqrt.
+ shy := sh2y sqrt.
+ ch2y := 1 + sh2y.
+ chy := ch2y sqrt.
+ y := imaginary copySignTo: shy arcSinh.
+ x := (real / chy) arcSin.
+ ^self class real: x imaginary: y
+
+ "
+ (10 + 4i) sin -> (-14.85625516387525498-22.89819255096375875i)
+ "
+!
+
+arcSinh
+ "Answer receiver's area hyperbolic sine.
+ That is the inverse function of sinh."
+
+ "Some possible implementation:
+
+ ^imaginary * real < 0
+ ifTrue: [(self + (self * self + 1) sqrt) ln]
+ ifFalse: [(self - (self * self + 1) sqrt) ln]"
+
+ ^self i arcSin i negated
+
+ "
+ (10 + 4i) sinh arcSinh -> (-10.0-0.8584073464102067614i)
+ "
+!
- ^ self class real: real sign imaginary: imaginary sign
+arcTan
+ "Answer the arc tangent of the receiver.
+ This is the inverse function of tan."
+
+ | r2 |
+
+ r2 := self squaredNorm.
+ ^ self class
+ real: (real * 2 arcTan: 1 - r2) / 2
+ imaginary: ((r2 + (imaginary * 2) + 1) / (r2 - (imaginary * 2) + 1)) ln / 4
+
+ "
+ (1+4i) tan arcTan -> (1.000000000000000043+4.000000000000004772i)
+ (10+4i) tan arcTan -> (0.5752220392306203171+4.00000000000000004i)
+ "
+!
+
+arcTan:denominator
+ "Answer the four quadrants arc tangent of receiver over denominator."
+
+ |res|
+
+ denominator isZero ifTrue:[
+ self isZero ifTrue:[
+ "shouldn't it be an error ? ^DomainError signal: '0 arcTan: 0'"
+ ^ self class real:0 imaginary:0
+ ] ifFalse:[
+ ^ self class real:Float pi / (real copySignTo:2) imaginary:0
+ ]
+ ].
+
+ res := (self / denominator) arcTan.
+ denominator real < 0 ifTrue:[
+ res := res + Float pi
+ ].
+ res real > Float pi ifTrue:[
+ res := res - (Float pi * 2)
+ ].
+ ^ res
+!
+
+arcTanh
+ "Answer the receiver's area hyperbolic tangent.
+ That is the inverse function of tanh."
+
+ "Some other possible implementation:
+
+ ^((1 + self) / (1 - self)) ln / 2"
+
+ ^ self i arcTan i negated
+
+ "
+ (10+4i) tanh arcTanh
+ "
+!
+
+cos
+ "Answer the receiver's cosine."
+
+ "/ ^ self i cosh
+ ^ self class
+ real:(real cos * imaginary cosh)
+ imaginary:(-1 * real sin * imaginary sinh)
+
+ "
+ (10+4i) cos -> (-22.91356068209214107+14.84629106966035727i)
+ (10+4i) i cosh (-22.91356068209214107+14.84629106966035727i)
+ "
+!
+
+cosh
+ "Answer the receiver's hyperbolic cosine.
+ Hyperbolic cosine is defined by same power series expansion as for real numbers,
+ that is in term of exponential:
+ ^ (self exp + self negated exp) / 2.
+ This implementation avoids creating intermediate objects."
+
+ ^self class
+ real: real cosh * imaginary cos
+ imaginary: real sinh * imaginary sin
+
+ "
+ (10+4i) cosh -> (-7198.729443310666079-8334.842120982836036i)
+ "
+!
+
+sin
+ "Answer the receiver's sine."
+
+ "/ alternative:
+ "/ ^ self i sinh i negated
+
+ ^ self class
+ real:(real sin * imaginary cosh)
+ imaginary:(real cos * imaginary sinh)
+
+ "
+ (10+4i) sin -> (-14.85625516387525498-22.89819255096375875i)
+ "
+!
+
+sinh
+ "Answer the receiver's hyperbolic sine.
+ Hyperbolic sine is defined by same power series expansion as for real numbers,
+ that is in term of exponential:
+ ^ (self exp - self negated exp) / 2.
+ This implementation avoids creating intermediate objects."
+
+ ^ self class
+ real: real sinh * imaginary cos
+ imaginary: real cosh * imaginary sin
+
+ "
+ 10 sinh -> 11013.23287470339338
+ (10+4i) sinh (-7198.729413635291604-8334.842155341616165i)
+ "
+!
+
+tan
+ "Answer the receivers tangent."
+
+ ^ self sin / self cos
+
+ "
+ (10+4i) tan -> (0.0006123503000121616919+0.9997260574022127583i)
+ (10+4i) sin -> (-14.85625516387525498-22.89819255096375875i)
+ (10+4i) cos -> (-22.91356068209214107+14.84629106966035727i)
+ "
+!
+
+tanh
+ "Answer the receiver's hyperbolic tangent."
+
+ "Some possible implementations are:
+ ^self sinh / self cosh
+
+ | tr ti |
+ tr := real tanh.
+ ti := imaginary tan i.
+ ^(tr + ti) / (tr * ti + 1)
+ "
+
+ ^ self i tan i negated
+
+ "
+ (10+4i) tanh
+ "
! !
!Complex methodsFor:'truncation & rounding'!