"{ 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

 any significant part of it which is separately transmitted or stored:

	* You must ensure that every copy includes this notice, and that

	  source and author(s) of the material are acknowledged.

	* These conditions must be imposed on anyone who receives a copy.

	* The material shall not be used for commercial gain without the prior

	  written consent of the author(s).

 For more information about the Manchester Goodies Library (from which

 this file was distributed) send e-mail:

	To: goodies-lib@cs.man.ac.uk

	Subject: help

 This is an additional goody-class, which is NOT covered by the

 ST/X license. It has been packaged with the ST/X distribution to

 make your live easier instead. NO WARRANTY.

"

"{ Package: 'stx:libbasic' }"

"{ NameSpace: Smalltalk }"

ArithmeticValue subclass:#Complex

	instanceVariableNames:'real imaginary'

	classVariableNames:'ComplexOne ComplexZero'

	poolDictionaries:''

	category:'Magnitude-Numbers'

!

!Complex class methodsFor:'documentation'!

copyright

"

 This is a Manchester Goodie.  It is distributed freely on condition

 that you observe these conditions in respect of the whole Goodie, and on

 any significant part of it which is separately transmitted or stored:

	* You must ensure that every copy includes this notice, and that

	  source and author(s) of the material are acknowledged.

	* These conditions must be imposed on anyone who receives a copy.

	* The material shall not be used for commercial gain without the prior

	  written consent of the author(s).

 For more information about the Manchester Goodies Library (from which

 this file was distributed) send e-mail:

	To: goodies-lib@cs.man.ac.uk

	Subject: help

 This is an additional goody-class, which is NOT covered by the

 ST/X license. It has been packaged with the ST/X distribution to

 make your live easier instead. NO WARRANTY.

"

!

documentation

"

    This class implements complex numbers.

    A complex number has real and imaginary parts which must be manipulated simultaneously

    in any numeric processing.

    Complex numbers can be used in many of the same places that regular numbers

    can be used with one major exception of comparisons, since complex numbers cannot

    be directly compared for size

    (except through lengths of vectors (see absolute value)).

    [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)

    [Constructors:]                         

        5 i

        6 + 7 i.

        5.6 - 8 i.

        Complex real: 10 imaginary: 5.

        Complex abs: 5 arg: (Float pi / 4)

        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.

        additions (trigonometric) and fixes by cg.

"

!

examples

"

    -25 sqrt                                  -> error

    Complex trapImaginary:[ -25 sqrt ]        -> 5i (0.0+5.0i)

    Complex trapImaginary:[ -25 integerSqrt ] -> 5i (0+5i)

    (Complex trapImaginary:[ -5397346292805549782720214077673687804022210808238353958670041357153884304 integerSqrt ])

        squared

    1 + 3i

    Number i + 1

    1 + Number i

    1i * 1i

    Number i * Number i

    (5 % 7) real

    (5 % 7) imaginary

    (5 % 7) = 5

    (5 % 0) = 5

    (5.0 % 0) = 5

    (1 % 0) + (2 % 0)

    (1 % 0) + (0 % 2)

    (1 % 0) + (2 % 3)

    (1 % 0) * (2 % 0)

    (1 % 0) * (0 % 2)

    (1 % 0) * (2 % 3)

    (1 % 2) + 2

    (1 % 2) * 2

    2 + (1 % 2)

    2 * (1 % 2)

    (Number i raisedTo:-3) -> Number i

    (Number i raisedTo:-2) -> -1

    (Number i raisedTo:-1) -> Number i negated

    (Number i raisedTo:0) -> 1

    (Number i raisedTo:1) -> Number i

    (Number i raisedTo:2) -> -1

    (Number i raisedTo:3) -> Number i negated

    (Number i raisedTo:4) -> 1

    (Number i raisedTo:6) -> -1

    3 raisedTo:Number i

    3 i raisedTo:Number i

    5i * 5i     -> -25

    5i squared  -> -25 

"

! !

!Complex class methodsFor:'instance creation'!

abs:aNumber1 arg:aNumber2

    ^ self 

        real:(aNumber1 * aNumber2 cos)

        imaginary:(aNumber1 * aNumber2 sin).

    "

     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

    "

     Complex fromReal:1.0 

    "

    "Modified (comment): / 12-06-2017 / 20:42:56 / cg"

!

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

    "

     Complex imaginary:1.0 

     (0.0 % 1.0)

    "

    "Modified (comment): / 12-06-2017 / 20:44:51 / cg"

!

real:aNumber

    "Create a new complex number from the given real number."

    ^ self basicNew setReal:aNumber setImaginary:0

    "

     Complex real:1.0 

    "

    "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."

    v = 0 ifTrue: [^ u].

    ^ self basicNew setReal:u setImaginary:v

    "

     Complex real:1.0 imaginary:2.0

     (1.0 % 2.0)

    "

    "Modified (comment): / 12-06-2017 / 20:44:32 / cg"

! !

!Complex class methodsFor:'coercing & converting'!

coerce: aNumber

    "convert the argument aNumber into an instance of the receiver's class and return it."

    ^ aNumber asComplex

! !

!Complex class methodsFor:'constants access'!

unity

    "Answer the value which allows, for any given arithmetic value, the following to be true:

        aNumber * aNumber class unity = aNumber

     This must be true regardless of how a given subclass chooses to define #*"

    ComplexOne isNil ifTrue:[

        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

     This must be true regardless of how a given subclass chooses to define #+"

    ComplexZero isNil ifTrue:[

        ComplexZero := self basicNew setReal:0 setImaginary:0

    ].

    ^ ComplexZero

! !

!Complex methodsFor:'accessing'!

imaginary

    "Return the imaginary part of the complex number."

    ^ imaginary

!

imaginaryPart

    "Return the imaginary part of the complex number.

     An alias for imaginary (for compatibility with other complex implementations)"

    ^ imaginary

!

real

    "Return the real part of the complex number."

    ^ real

!

realPart

    "Return the real part of the complex number.

     An alias for real (for compatibility with other complex implementations)"

    ^ real

! !

!Complex methodsFor:'arithmetic'!

* 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:[

"/        u := aNumber real.

"/        v := aNumber imaginary.

"/        r := (real * u) - (imaginary * v).

"/        i  := (real * v) + (imaginary * u).

"/        i = 0 ifTrue:[ ^ r ].

"/        ^ Complex real:r imaginary:i

"/    ].

    ^ aNumber productFromComplex:self.

    "Modified: / 8.7.1998 / 12:12:37 / cg"

!

+ 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: [

"/        r := aNumber real + real.

"/        i := aNumber imaginary + imaginary.

"/        i = 0 ifTrue:[ ^ r ].

"/        ^ Complex real:r imaginary:i

"/    ].

    ^ aNumber sumFromComplex:self.

    "Modified: / 8.7.1998 / 12:15:42 / cg"

!

- 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: [

"/        r := real - aNumber real.

"/        i := imaginary - aNumber imaginary.

"/        i = 0 ifTrue:[ ^ r ].

"/        ^ Complex real:r imaginary:i.

"/    ].

    ^ aNumber differenceFromComplex:self.

    "Modified: / 8.7.1998 / 12:15:38 / cg"

!

/ 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:[

"/        u := aNumber real.

"/        v := aNumber imaginary.

"/        denom := u * u + (v * v).

"/        r := u * real + (v * imaginary) / denom.

"/        i := u * imaginary - (v * real) / denom.

"/        i = 0 ifTrue:[ ^ r ].

"/        ^ Complex real:r imaginary:i

"/    ].

    ^ aNumber quotientFromComplex:self.

    "Modified: / 8.7.1998 / 12:15:34 / cg"

!

abs

    "Return the magnitude (or absolute value) of the complex number

     that's the distance from zero (the origin in the complex plane)."

    ^ (real * real + (imaginary * imaginary)) sqrt

    "

     (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

    "

!

conjugated

    "Return the complex conjugate of this complex number

     (i.e. with imaginary part negated)."

    ^ self class

        real: real

        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.