Object subclass: #DhbLanczosFormula
    instanceVariableNames: 'coefficients sqrt2Pi '
    classVariableNames: 'UniqueInstance '
    poolDictionaries: ''!

SubApplication subclass: #DhbGammaFunction
    instanceVariableNames: ''
    classVariableNames: ''
    poolDictionaries: ''!



!DhbLanczosFormula class publicMethods !

new
		"Answer a unique instance. Create it if it does not exist.
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 4/1/99 "
	UniqueInstance isNil
		ifTrue: [ UniqueInstance := super new.
					 UniqueInstance initialize.
					].
	^UniqueInstance! !

!DhbLanczosFormula publicMethods !

gamma: aNumber
		"
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 1/3/99 "
	^( self leadingFactor: aNumber) exp * ( self series: aNumber) * sqrt2Pi / aNumber!

logGamma: aNumber
		"
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 1/3/99 "
	^( self leadingFactor: aNumber) + ( ( self series: aNumber) * sqrt2Pi / aNumber) ln! !

!DhbLanczosFormula privateMethods !

initialize
		"Private - Initialize the coefficients of the receiver.
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 1/3/99 "
	sqrt2Pi := ( Float pi * 2) sqrt.
	coefficients := #( 76.18009172947146 -86.50532032941677 24.01409824083091 -1.231739572450155 0.1208650973866179e-2 -0.5395239384953e-5).
	^self!

leadingFactor: aNumber
		"Private - Answers the log of the leading factor in Lanczos' formula.
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 1/3/99 "
	| temp |
	temp := aNumber + 5.5.
	^( temp ln * ( aNumber + 0.5) - temp)!

series: aNumber
		"Private - Answer the value of the series of Lanczos' formula.
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 1/3/99 "
	| term |
	term := aNumber.
	^coefficients inject: 1.000000000190015
						into: [ :sum :each | term := term + 1. each / term + sum]! !

!Integer publicMethods !

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

!Number publicMethods !

gamma
		"Compute the Gamma function for the receiver.
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 11/2/99 "
	^self > 1
		ifTrue: [ ^DhbLanczosFormula new gamma: self]
		ifFalse:[ self < 0
						ifTrue: [ Float pi / ( ( Float pi * self) sin * ( 1 - self) gamma)]
						ifFalse:[ ( DhbLanczosFormula new gamma: ( self + 1)) / self]
					]!

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: [ DhbLanczosFormula new logGamma: self]
		ifFalse:[ self > 0
						ifTrue: [ ( DhbLanczosFormula new logGamma: ( self + 1) ) - self ln ]
						ifFalse: [ ^self error: 'Argument for the log gamma function must be positive']
					]! !

DhbLanczosFormula initializeAfterLoad!

DhbGammaFunction initializeAfterLoad!