"======================================================================
|
|   Numerical methods - Additions to the base classes
|
|
 ======================================================================"


"======================================================================
|
| Written by Didier Besset.
|
| This file is distributed together with GNU Smalltalk.
|
 ======================================================================"

!Smalltalk.Integer methodsFor: 'numerics'!

gamma
    "Compute the Gamma function for the receiver.
     (c) Copyrights Didier BESSET, 1999, all rights reserved.
     Initial code: 11/2/99 "

    self > 0
        ifFalse:[ ^self error: 'Attempt to compute the Gamma function of a non-positive integer'].
    ^( self - 1) factorial!

random
    "Answer a random integer between 0 and the receiver.
     (c) Copyrights Didier BESSET, 1999, all rights reserved.
     Initial code: 15/2/99 "

    ^Dhb.DhbMitchellMooreGenerator new integerValue: self! !

!Smalltalk.Number class methodsFor: 'numerics'!

random
    "Answers a random number between 0 and the receiver
     (c) Copyrights Didier BESSET, 1999, all rights reserved.
     Initial code: 17/2/99 "

    ^Dhb.DhbMitchellMooreGenerator new floatValue * self! !

!Smalltalk.Number methodsFor: 'numerics'!

addPolynomial: aPolynomial
    "(c) Copyrights Didier BESSET, 1999, all rights reserved.
     Initial code: 19/4/99 "

    ^aPolynomial addNumber: self!

asLimitedPrecisionReal
    "Convert the receiver to an instance of
    some subclass of LimitedPrecisionReal.
    This method defines what the default is."

    ^self asFloat!

beta: aNumber
    "Computes the beta function of the receiver and aNumber
     (c) Copyrights Didier BESSET, 1999, all rights reserved.
     Initial code: 1/3/99 "

    ^( self logBeta: aNumber) exp!

dividingPolynomial: aPolynomial
    "(c) Copyrights Didier BESSET, 1999, all rights reserved.
     Initial code: 17/4/99 "

    ^aPolynomial timesNumber: (1 / self)!

equalsTo: aNumber
    " (c) Copyrights Didier BESSET, 1999, all rights reserved.
     Initial code: 21/4/99 "

    ^self relativelyEqualsTo: aNumber upTo: Dhb.DhbFloatingPointMachine new defaultNumericalPrecision
!

errorFunction
    "Answer the error function for the receiver.
     (c) Copyrights Didier BESSET, 1999, all rights reserved.
     Initial code: 11/2/99 "

    ^Dhb.DhbErfApproximation new value: self
!

gamma
    "Compute the Gamma function for the receiver. 
    (c) Copyrights Didier BESSET, 1999, all rights reserved. 
    Initial code: 11/2/99"

    ^self > 1
        ifTrue: [^Dhb.DhbLanczosFormula new gamma: self]
        ifFalse: [self < 0
                ifTrue: [Float pi / ((Float pi * self) sin * (1 - self) gamma)]
                ifFalse: [(Dhb.DhbLanczosFormula new gamma: self + 1)
                        / self]]!

logBeta: aNumber
    "Computes the logarithm of the beta function of the receiver and aNumber
     (c) Copyrights Didier BESSET, 1999, all rights reserved.
     Initial code: 1/3/99 "

    ^self logGamma + aNumber logGamma - ( self + aNumber) logGamma!

logGamma
    "Computes the log of the Gamma function (for positive numbers only)
     (c) Copyrights Didier BESSET, 1999, all rights reserved.
     Initial code: 1/3/99 "

    ^self > 1
        ifTrue: [ Dhb.DhbLanczosFormula new logGamma: self]
        ifFalse:[ self > 0
                        ifTrue: [ ( Dhb.DhbLanczosFormula new logGamma: ( self + 1) ) - self ln ]
                        ifFalse: [ ^self error: 'Argument for the log gamma function must be positive']
        ]!

productWithMatrix: aMatrix
    "Answer a new matrix, product of aMatrix with the receiver.
     (c) Copyrights Didier BESSET, 1999, all rights reserved.
     Initial code: 11/2/99 "

    ^aMatrix class rows: ( aMatrix rowsCollect: [ :each | each * self])!

productWithVector: aVector
    "Answers a new vector product of the receiver with aVector.
     (c) Copyrights Didier BESSET, 1999, all rights reserved.
     Initial code: 11/2/99 "

    ^aVector collect: [ :each | each * self]!

random
    "Answers a random number distributed between 0 and the receiver.
     (c) Copyrights Didier BESSET, 1999, all rights reserved.
     Initial code: 17/2/99 "

    ^self class random * self!

relativelyEqualsTo: aNumber upTo: aSmallNumber
    "(c) Copyrights Didier BESSET, 1999, all rights reserved.
     Initial code: 21/4/99 "

    | norm |
    norm := self abs max: aNumber abs.
    ^norm <= Dhb.DhbFloatingPointMachine new defaultNumericalPrecision
        or: [ (self - aNumber) abs < ( aSmallNumber * norm)]!

subtractToPolynomial: aPolynomial
    "(c) Copyrights Didier BESSET, 1999, all rights reserved.
     Initial code: 19/4/99 "

    ^aPolynomial addNumber: self negated!

timesPolynomial: aPolynomial
    "(c) Copyrights Didier BESSET, 1999, all rights reserved.
     Initial code: 17/4/99 "

    ^aPolynomial timesNumber: self! !

!Smalltalk.Point methodsFor: 'numerics'!

extentFromBottomLeft: aPoint
    "(c) Copyrights Didier BESSET, 1998, all rights reserved
     Initial code: 21/4/98 "

    ^Rectangle origin: (self - (0 @ (aPoint isInteger ifTrue: [ aPoint ]
                                                      ifFalse:[ aPoint y ])))
               extent: aPoint!

extentFromBottomRight: aPoint
    "(c) Copyrights Didier BESSET, 1998, all rights reserved
     Initial code: 21/4/98 "

    ^Rectangle origin: (self - aPoint) extent: aPoint!

extentFromCenter: aPoint
    "(c) Copyrights Didier BESSET, 1998, all rights reserved
     Initial code: 21/4/98 "

    ^Rectangle origin: ( self - ( aPoint // 2)) extent: aPoint!

extentFromTopLeft: aPoint
    "(c) Copyrights Didier BESSET, 1998, all rights reserved
     Initial code: 21/4/98 "

    ^Rectangle origin: self extent: aPoint!

extentFromTopRight: aPoint
    "(c) Copyrights Didier BESSET, 1998, all rights reserved
     Initial code: 21/4/98 "

    ^Rectangle origin: (self - ((aPoint isInteger ifTrue: [ aPoint ]
                                                  ifFalse:[ aPoint x ]) @ 0))
               extent: aPoint! !

!Smalltalk.Rectangle methodsFor: 'numerics'!

positiveRectangle
    "(c) Copyrights Didier BESSET, 1998, all rights reserved
     Initial code: 21/4/98 "

    ^( origin min: corner) corner: ( origin max: corner)! !

!Smalltalk.Collection methodsFor: 'numerics'!

asVector
    ^(Dhb.DhbVector new: self size)
        replaceElementsFrom: 1 to: self size withArray: self startingAt: 1
! !

!Dhb.DhbPolynomial methodsFor: 'numerics'!

generality
    ^nil
! !

!Dhb.DhbVector methodsFor: 'numerics'!

generality
    ^nil
! !

