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

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

    NOTE (from the original author): 
        Although Complex seems 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.

    NOTE (from porter): 
        moved to Number hierarchy.
        Makes live of users much easier (isNumber, reading, etc)
        
    [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.
        ^ 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 |

    absReal := real abs.
    absImag := imaginary abs.

    absReal >= absImag ifTrue: [
	multiplicand := absReal.
	quotient := imaginary / real
    ] ifFalse: [
	multiplicand := absImag.
	quotient := real / imaginary
    ].
    ^ multiplicand * ((1 + (quotient * quotient)) sqrt)
!

negated
    "return a new complex with both real and imaginary parts negated"

    ^ self class
        real: real negated
        imaginary: imaginary negated
! !

!Complex methodsFor:'coercing & 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
!

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'!

< aNumber
    "raises an error - complex numbers are not well ordered"

    ^ Number
	raise: #unorderedSignal
	receiver: self
	selector: #<
	arg: aNumber
	errorString: 'Complex numbers are not well ordered'

    "
     1 < (2 % 2)
     (2 % 2) < 1
    "
!

= anObject
    "return true, if the argument represents the same numeric value
     as the receiver, false otherwise."

    anObject class == self class ifTrue:[
        ^ (real = anObject real) and:[ (imaginary = anObject imaginary)]
    ].
    ^ anObject equalFromComplex:self

    "
     (Complex real:1.0 imaginary:2.0) = (Complex real:1.0 imaginary:2.0)
     (Complex real:1.0 imaginary:0) = 1.0
    "

    "Modified (comment): / 12-06-2017 / 20:43:41 / cg"
    "Modified: / 26-05-2019 / 10:04:33 / Claus Gittinger"
!

hash
    "Hash is implemented because = is implemented."

    ^ (real hash) bitXor:(imaginary hash bitShift:16)

    "
     (1+0i) hash
     (1+1i) hash
    "
! !

!Complex methodsFor:'double dispatching'!

differenceFromComplex:aComplex
    "Return the difference of the argument, aComplex and the receiver."

    | r i |

    r := aComplex real - real.
    i := aComplex imaginary - imaginary.
    i = 0 ifTrue:[ ^ r ].
    ^ self class real:r imaginary:i.
!

differenceFromFixedPoint: aFixedPoint
    "Return the difference of the argument, aFixedPoint and the receiver."

    ^ aFixedPoint asComplex - self
!

differenceFromFloat:aFloat
    "Return the difference of the argument, aFloat and the receiver."

    "/ ^ aFloat asComplex - self

    | r |

    r := aFloat - real.
    imaginary = 0 ifTrue:[ ^ r ].
    ^ self class real:r imaginary:imaginary negated

    "
     (1 % 1) - 1.0
     1.0 - (1 % 1)
    "
!

differenceFromFraction: aFraction
    "Return the difference of the argument, aFraction and the receiver."

    ^ aFraction asComplex - self
!

differenceFromInteger: anInteger
    "Return the difference of the argument, anInteger and the receiver."

    ^ anInteger asComplex - self
!

equalFromComplex:aComplex
    "return true if aComplex represents the same number as myself"

    ^ (aComplex real = real) and:[aComplex imaginary = imaginary]

    "Modified (comment): / 12-06-2017 / 20:41:04 / cg"
!

equalFromFloat:aFloat
    "return true if aFloat represents the same number as myself"

    imaginary = 0 ifFalse:[^ false].
    ^ real = aFloat

    "Modified (comment): / 12-06-2017 / 20:41:00 / cg"
!

productFromComplex:aComplex
    "Return the product of the receiver and the argument, aComplex."

    | u v r i |

    u := aComplex real.
    v := aComplex imaginary.
    r := (real * u) - (imaginary * v).
    i  := (real * v) + (imaginary * u).
    i = 0 ifTrue:[ ^ r ].
    ^ self class real:r imaginary:i
!

productFromFixedPoint: aFixedPoint
    "Return the product of the receiver and the argument, aFixedPoint."

    ^ aFixedPoint asComplex * self
!

productFromFloat: aFloat
    "Return the product of the receiver and the argument, aFloat."

    "/  ^ aFloat asComplex * self

    | u r i |

    u := aFloat.
    r := (real * aFloat).
    i  := (imaginary * aFloat).
    i = 0 ifTrue:[ ^ r ].
    ^ self class real:r imaginary:i

    "
     (1 % 1) * 2.0
     (1 % 1) * 0.0
     2.0 * (1 % 1)
    "
!

productFromFraction: aFraction
    "Return the product of the receiver and the argument, aFraction."

    ^ aFraction asComplex * self
!

productFromInteger: anInteger
    "sent when an integer does not know how to multiply the receiver, a complex.
     Return the product of the receiver and the argument, anInteger."

    ^ anInteger asComplex * self

    "Modified (comment): / 12-06-2017 / 20:30:20 / cg"
!

quotientFromComplex:aComplex
    "Return the quotient of the argument, aComplex and the receiver."

    | denom nr ni r i |

    nr := aComplex real.
    ni := aComplex imaginary.
    denom := real * real + (imaginary * imaginary).
    r := (real * nr + (imaginary * ni)) / denom.
    i := (real * ni - (imaginary * nr)) / denom.
    i = 0 ifTrue:[ ^ r ].
    ^ self class real:r imaginary:i

"/ is the stuff below better ?
"/    "Implement complex division (a + ib) / (c + id).
"/     Due to double dispatch, in this routine
"/        self = (c + id) and aComplex = (a + ib)."
"/
"/    | quotient denominator |
"/
"/    self realPart abs >= (self imaginaryPart abs)
"/        ifTrue: [
"/            quotient := self imaginaryPart / self realPart.
"/            denominator := self realPart + (self imaginaryPart * quotient).
"/            ^ Complex
"/                real: (aComplex realPart + (aComplex imaginaryPart * quotient)) / denominator
"/                imaginary: (aComplex imaginaryPart - (aComplex realPart * quotient)) / denominator ]
"/        ifFalse: [
"/            quotient := self realPart / self imaginaryPart.
"/            denominator := (self realPart * quotient) + self imaginaryPart.
"/            ^ Complex
"/                real: ((aComplex realPart * quotient) + aComplex imaginaryPart) / denominator
"/                imaginary: ((aComplex imaginaryPart * quotient) - aComplex realPart) / denominator ]
!

quotientFromFixedPoint:aFixedPoint
    "Return the quotient of the argument, aFixedPoint and the receiver."

    ^ aFixedPoint asComplex / self
!

quotientFromFloat:aFloat
    "Return the quotient of the argument, aFloat and the receiver."

    ^ aFloat asComplex / self
!

quotientFromFraction:aFraction
    "Return the quotient of the argument, aFraction and the receiver."

    ^ aFraction asComplex / self
!

quotientFromInteger:anInteger
    "Return the quotient of the argument, anInteger and the receiver."

    ^ anInteger asComplex / self
!

raisedFromFloat:aNumber
    ^ self raisedFromNumber:aNumber
    
    "
     2 raisedTo:(2 + 2i)
     2 ** (2 + 2i)

     2.0 raisedTo:(2 + 2i)
     2.0 ** (2 + 2i)
    "

    "Created: / 26-01-2019 / 10:56:06 / Claus Gittinger"
!

raisedFromNumber:aNumber
    "see http://www.math.toronto.edu/mathnet/questionCorner/complexexp.html"

    "/ a ^ (b+i*c) = a^b * ( cos(c*ln(a)) + i*sin(c*ln(a)) )

    |cLNa a_b|

    cLNa := imaginary * aNumber ln.
    a_b := aNumber raisedTo:real.
    
    ^ Complex
        real:(a_b * cLNa cos)
        imaginary:(a_b * cLNa sin)

    "
     2 raisedTo:(2 + 2i)
     2 ** (2 + 2i)
    "

    "Created: / 01-07-2017 / 20:33:47 / cg"
    "Modified (comment): / 03-07-2017 / 14:08:45 / cg"
!

sumFromComplex:aComplex
    "Return the sum of the receiver and the argument, aComplex."

    | r i |

    r := aComplex real + real.
    i := aComplex imaginary + imaginary.
    i = 0 ifTrue:[ ^ r ].
    ^ self class real:r imaginary:i
!

sumFromFixedPoint: aFixedPoint
    "Return the sum of the receiver and the argument, aFixedPoint."

    ^ aFixedPoint asComplex + self
!

sumFromFloat: aFloat
    "Return the sum of the receiver and the argument, aFloat."

    "/ ^ aFloat asComplex + self

    | r |

    r := aFloat + real.
    imaginary = 0 ifTrue:[ ^ r ].
    ^ self class real:r imaginary:imaginary

    "
     (1 % 1) + 1.0
     1.0 + (1 % 1)
    "
!

sumFromFraction: aFraction
    "Return the sum of the receiver and the argument, aFraction."

    ^ aFraction asComplex + self
!

sumFromInteger: anInteger
    "Return the sum of the receiver and the argument, anInteger."

    ^ anInteger asComplex + self
! !

!Complex methodsFor:'mathematical functions'!

angle
    "Return the radian angle for this Complex number."

    real < 0 ifTrue: [
	imaginary < 0 ifTrue: [
	    ^ (imaginary / real) arcTan - Float pi
	].
	^ Float pi + (imaginary / real) arcTan
    ].
    ^ (imaginary / real) arcTan

    "
     (1 % 1) angle radiansToDegrees
    "
!

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

    ^ (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
    "Return the square root of the receiver"

    | w quotient absReal absImag |

    ((real = 0) and: [ imaginary = 0 ]) ifTrue: [
        ^ self class zero
    ].
    absReal := real abs.
    absImag := imaginary abs.

    absReal >= absImag ifTrue:[
        quotient := imaginary / real.
        w := (absReal sqrt) * (((1 + (1 + (quotient * quotient)) sqrt) / 2) sqrt)
    ] ifFalse: [
        quotient := real / imaginary.
        w := (absImag sqrt) * (((quotient abs + (1 + (quotient * quotient)) sqrt) / 2) sqrt)
    ].

    real >= 0 ifTrue:[
        ^ self class real: w imaginary: (imaginary / (2 * w))
    ].
    imaginary >= 0 ifTrue: [
        ^ self class real: absImag / (2 * w) imaginary: w
    ].
    ^ self class real: absImag / (2 * w) imaginary: -1 * w
!

sqrt_bad
    "Return the square root of the receiver"

    | u v |

    (imaginary = 0 and: [real >= 0]) ifTrue: [^ real sqrt].
    v := ((self abs - real) / 2) sqrt.
    u := imaginary / 2 / v.
    ^ self class real: u imaginary: v

    "
     -4 asComplex sqrt
     4 asComplex sqrt
    "
    "
     -4 asComplex sqrt squared
    "
!

squaredNorm
    "Answer the square of receiver's norm."

    ^real * real + (imaginary * imaginary)
! !

!Complex methodsFor:'printing & storing'!

displayOn: aGCOrStream
    "/ what a kludge - Dolphin and Squeak mean: printOn: a stream;
    "/ old ST80 means: draw-yourself on a GC.
    (aGCOrStream isStream) ifFalse:[
        ^ super displayOn:aGCOrStream
    ].

    aGCOrStream nextPut: $(.
    self realPart printOn: aGCOrStream.
    self imaginaryPart >= 0
        ifTrue: [ aGCOrStream nextPut: $+ ]
        ifFalse: [ aGCOrStream nextPut: $- ].
    self imaginaryPart abs printOn: aGCOrStream.
    aGCOrStream nextPutAll: 'i)'

    "
     Complex real:1 imaginary:1
     (Complex real:1 imaginary:1) printString
     (Complex real:1 imaginary:1) displayString
    "

    "Modified (format): / 22-02-2017 / 17:01:02 / cg"
    "Modified (comment): / 26-06-2018 / 21:00:28 / Claus Gittinger"
!

printOn: aStream
    aStream nextPut: $(.
    real storeOn: aStream.
    aStream nextPutAll: '%'.
    imaginary storeOn: aStream.
    aStream nextPut: $).
!

printString
    ^ '(' , real printString, '%', imaginary printString, ')'
!

storeOn: aStream
    self printOn:aStream
! !

!Complex methodsFor:'private'!

setReal: u setImaginary: v
    real := u.
    imaginary := v.
! !

!Complex methodsFor:'testing'!

isComplex
    "Answer whether the receiver has an imaginary part
     (i.e. if it is a complex number). Always true here."

    ^ true
!

isReal
    "Return true if this Complex number has a zero imaginary part."

    ^ imaginary = 0
!

isZero
    "Answer whether 'self = self class zero'.
     We can't use #= because #= is defined in terms of #isZero"

    ^real isZero and: [imaginary isZero]
!

sign
    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)
    "
!

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 receiver's 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'!

ceiling
    "blocked: complex numbers have no ceiling"

    ^ self shouldNotImplement
!

floor
    "blocked: complex numbers have no floor"

    ^ self shouldNotImplement
! !

!Complex class methodsFor:'documentation'!

version
    ^ '$Header$'
!

version_CVS
    ^ '$Header$'
! !