Object subclass: #DhbLagrangeInterpolator
    instanceVariableNames: 'pointCollection '
    classVariableNames: ''
    poolDictionaries: ''!

DhbLagrangeInterpolator subclass: #DhbNevilleInterpolator
    instanceVariableNames: 'leftErrors rightErrors '
    classVariableNames: ''
    poolDictionaries: ''!

DhbNevilleInterpolator subclass: #DhbBulirschStoerInterpolator
    instanceVariableNames: ''
    classVariableNames: ''
    poolDictionaries: ''!

DhbLagrangeInterpolator subclass: #DhbNewtonInterpolator
    instanceVariableNames: 'coefficients '
    classVariableNames: ''
    poolDictionaries: ''!

DhbNewtonInterpolator subclass: #DhbSplineInterpolator
    instanceVariableNames: 'startPointDerivative endPointDerivative '
    classVariableNames: ''
    poolDictionaries: ''!

Object subclass: #DhbPolynomial
    instanceVariableNames: 'coefficients '
    classVariableNames: ''
    poolDictionaries: ''!

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



!DhbBulirschStoerInterpolator privateMethods !

computeDifference: aNumber at: anInteger1 order: anInteger2
		"Private -
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 28/5/99 "
	| diff ratio |
	ratio := ( ( self xPointAt: anInteger1) - aNumber) * ( rightErrors at: anInteger1)
							/ ( (  self xPointAt: ( anInteger1 + anInteger2)) - aNumber).
	diff := ( ( leftErrors at: ( anInteger1 + 1)) - ( rightErrors at: anInteger1))
							/ ( ratio - ( leftErrors at: ( anInteger1 + 1))).
	rightErrors at: anInteger1 put: ( leftErrors at: ( anInteger1 + 1)) * diff.	
	leftErrors at: anInteger1 put: ratio * diff.! !

!DhbLagrangeInterpolator class publicMethods !

new
		"Create a new instance of the receiver without points. Points must be added with add:
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 11/2/99 "
	^super new initialize!

points: aCollectionOfPoints
		"Create a new instance of the receiver with given points.
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 11/2/99 "
	^self new initialize: aCollectionOfPoints! !

!DhbLagrangeInterpolator publicMethods !

add: aPoint
		"Add a point to the collection of points.
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 11/2/99 "
	^pointCollection add: aPoint!

size
		"(c) Copyrights Didier BESSET, 2000, all rights reserved.
		 Initial code: 3/12/00 "
	^pointCollection size!

value: aNumber
		"Compute the value of the Lagrange interpolation polynomial on the receiver's points at aNumber.
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 11/2/99 "
	| norm dx products answer size |
	norm := 1.
	size := pointCollection size.
	products := Array new: size.
	products atAllPut: 1.
	1 to: size
		do: [ :n |
			  dx := aNumber - ( self xPointAt: n).
			  dx = 0
				ifTrue: [ ^( self yPointAt: n)].
			  norm := norm * dx.
			  1 to: size
				do: [ :m |
					  m = n
						ifFalse:[ products at: m put: ( (( self xPointAt: m) - ( self xPointAt: n)) * ( products at: m))].
					].
			].
	answer := 0.
	1 to: size do:
		[ :n | answer := ( self yPointAt: n) / ( ( products at: n) * ( aNumber - ( self xPointAt: n))) + answer].
	^norm * answer
! !

!DhbLagrangeInterpolator privateMethods !

defaultSamplePoints
		"Private -
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 27/5/99 "
	^OrderedCollection new!

initialize
		"Private - Create an empty point collection for the receiver.
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 11/2/99 "
	^self initialize: self defaultSamplePoints!

initialize: aCollectionOfPoints
		"Private - Defines the collection of points for the receiver.
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 11/2/99 "
	pointCollection := aCollectionOfPoints.
	^self!

xPointAt: anInteger
		"Private -
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 19/4/99 "
	^( pointCollection at: anInteger) x!

yPointAt: anInteger
		"Private -
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 19/4/99 "
	^( pointCollection at: anInteger) y! !

!DhbNevilleInterpolator publicMethods !

value: aNumber
		"Compute the value of the Lagrange interpolation polynomial on the receiver's points at aNumber.
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 7/4/99 "
	^(self valueAndError: aNumber) first!

valueAndError: aNumber
		"Compute and return the interpolated value of the interpolation Lagranage polynomial
		 and its estimated error.
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 7/4/99 "
	| size nearestIndex answer error |
	nearestIndex := self initializeDifferences: aNumber.
	nearestIndex < 0
		ifTrue: [ ^Array with: ( self yPointAt: nearestIndex negated) with: 0].
	answer := leftErrors at: nearestIndex.
	nearestIndex := nearestIndex - 1.
	size := pointCollection size.
	1 to: ( size - 1) do:
		[ :m |
		  1 to: ( size - m) do:
			[ :n | self computeDifference: aNumber at: n order: m].
		  size - m > ( 2 * nearestIndex)
				ifTrue: [ error := leftErrors at: ( nearestIndex + 1) ]
				ifFalse:[ error := rightErrors at: ( nearestIndex).
							  nearestIndex := nearestIndex - 1.
							].
		  answer := answer + error.
		].
	^Array with: answer with: error abs! !

!DhbNevilleInterpolator privateMethods !

computeDifference: aNumber at: anInteger1 order: anInteger2
		"Private -
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 28/5/99 "
	| leftDist rightDist ratio |
	leftDist := ( self xPointAt: anInteger1) - aNumber.
	rightDist := (  self xPointAt: ( anInteger1 + anInteger2)) - aNumber.
	ratio := ( ( leftErrors at: ( anInteger1 + 1)) - ( rightErrors at: anInteger1)) / ( leftDist - rightDist).
	leftErrors at: anInteger1 put: ratio * leftDist.
	rightErrors at: anInteger1 put: ratio * rightDist.!

defaultSamplePoints
		"Private -
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 27/5/99 "
	^SortedCollection sortBlock: [ :a :b | a x < b x]!

initializeDifferences: aNumber
		"Private -
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 28/5/99 "
	| size nearestIndex dist minDist |
	size := pointCollection size.
	leftErrors size = size
		ifFalse:[ leftErrors := Array new: size.
				  rightErrors := Array new: size.
				].
	minDist := ( ( self xPointAt: 1) - aNumber) abs.
	nearestIndex := 1.
	leftErrors at: 1 put: ( self yPointAt: 1).
	rightErrors at: 1 put: leftErrors first.
	2 to: size do:
		[ :n |
		  dist := ( ( self xPointAt: n) - aNumber) abs.
		  dist < minDist
			ifTrue: [ dist = 0
						ifTrue: [ ^n negated].
					  nearestIndex := n.
					  minDist := dist.
					].
		 leftErrors at: n put: ( self yPointAt: n).
		 rightErrors at: n put: ( leftErrors at: n).
		].
	^nearestIndex
! !

!DhbNewtonInterpolator publicMethods !

add: aPoint
		"Add a point to the collection of points.
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 11/2/99 "
	self resetCoefficients.
	^super add: aPoint!

value: aNumber
		"Compute the value of the Lagrange interpolation polynomial on the receiver's points at aNumber.
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 11/2/99 "
	| answer size |
	coefficients isNil
		ifTrue: [ self computeCoefficients].
	size := coefficients size.
	answer := coefficients at: size.
	(size - 1) to: 1 by: -1
		do: [ :n | answer := answer * ( aNumber - (self xPointAt:  n)) + ( coefficients at: n)].
	^answer
! !

!DhbNewtonInterpolator privateMethods !

computeCoefficients
		"Private - Computes the coefficients for the receiver.
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 11/2/99 "
	| size k1 kn|
	size := pointCollection size.
	coefficients := ( 1 to: size) collect: [ :n | self yPointAt: n].
	1 to: (size - 1)
		do: [ :n |
			  size to: ( n + 1)  by: -1
				do: [ :k |
					  k1 := k - 1.
					  kn := k - n.
					  coefficients at: k put: ( (( coefficients at: k) - ( coefficients at: k1)) 
											/ ((self xPointAt: k) - (self xPointAt: kn))).
					].
			].
!

resetCoefficients
		"Private - Reset the coefficients of the receiver to force a new computation.
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 11/2/99 "
	coefficients := nil.! !

!DhbPolynomial class publicMethods !

coefficients: anArray
		" Creates a new instance with given coefficients
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 4/1/99 "
	^self new initialize: anArray reverse! !

!DhbPolynomial publicMethods !

* aNumberOrPolynomial
		"
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 17/4/99 "
	^aNumberOrPolynomial timesPolynomial: self!

+ aNumberOrPolynomial
		"(c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 19/4/99 "
	^aNumberOrPolynomial addPolynomial: self!

- aNumberOrPolynomial
		"(c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 19/4/99 "
	^aNumberOrPolynomial subtractToPolynomial: self!

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

addNumber: aNumber
		"(c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 19/4/99 "
	| newCoefficients |
	newCoefficients := coefficients reverse.
	newCoefficients at: 1 put: newCoefficients first + aNumber.
	^self class new: newCoefficients!

addPolynomial: aPolynomial
		"(c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 19/4/99 "
	^self class new: ( ( 0 to: (self degree max: aPolynomial degree)) collect: [ :n | ( aPolynomial at: n) + ( self at: n)])!

at: anInteger
		"Answers the coefficient of order anInteger.
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 17/4/99 "
	^anInteger < coefficients size
		ifTrue: [ coefficients at: ( coefficients size - anInteger)]
		ifFalse:[ 0]!

coefficients
		"(c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 21/5/99 "
	^coefficients deepCopy!

degree
		"(c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 17/4/99 "
	^coefficients size - 1!

derivative
		"Answer a new polynomial, derivative of the receiver.
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 4/1/99 "
	| n |
	n := coefficients size.
	^self class new: ( ( coefficients collect: [ :each | n := n - 1. each * n]) reverse copyFrom: 2 to: coefficients size)!

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

dividingPolynomialWithRemainder: aPolynomial
		"(c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 17/4/99 "
	| remainderCoefficients quotientCoefficients n m norm quotientDegree |
	n := self degree.
	m := aPolynomial degree.
	quotientDegree := m - n.
	quotientDegree < 0
		ifTrue: [ ^Array with: ( self class new: #(0)) with: aPolynomial].
	quotientCoefficients := Array new: quotientDegree + 1.
	remainderCoefficients := ( 0 to: m) collect: [ :k | aPolynomial at: k].
	norm := 1 / coefficients first.
	quotientDegree to: 0 by: -1
		do: [ :k | | x |
			  x := ( remainderCoefficients at: n + k + 1) * norm.
			  quotientCoefficients at: (quotientDegree + 1 - k) put: x.
			  (n + k - 1) to: k by: -1
				do: [ :j | 
				remainderCoefficients at: j + 1 put: 
							( ( remainderCoefficients at: j + 1) - ( x * (self at: j - k)))
				].
			].
	^Array with: ( self class new: quotientCoefficients reverse)
		   with: ( self class new: ( remainderCoefficients copyFrom: 1 to: n))!

integral
		"Answer a new polynomial, integral of the receiver with value 0 at x=0.
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 4/1/99 "
	^self integral: 0!

integral: aValue
		"Answer a new polynomial, integral of the receiver with given value at x=0.
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 4/1/99 "
	| n |
	n := coefficients size + 1.
	^self class new: ( ( coefficients collect: [ :each | n := n - 1. each / n]) copyWith: aValue) reverse!

printOn: aStream
		"Append to aStream a written representation of the receiver.
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 4/1/99 "
	| n firstNonZeroCoefficientPrinted  |
	n := 0.
	firstNonZeroCoefficientPrinted := false.
	coefficients reverse do:
		[ :each |
		  each = 0
			ifFalse:[ firstNonZeroCoefficientPrinted
							ifTrue: [ aStream space.
										 each < 0
											ifFalse:[ aStream nextPut: $+].
										 aStream space.
										]
							ifFalse:[ firstNonZeroCoefficientPrinted := true].
						  ( each = 1 and: [ n > 0])
							ifFalse:[ each printOn: aStream].
						  n > 0
							ifTrue: [ aStream nextPutAll: ' X'.
										 n > 1
											ifTrue: [ aStream nextPut: $^.
										 				 n printOn: aStream.
														].
										].
						].
		 n := n + 1.
		].!

subtractToPolynomial: aPolynomial
		"(c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 19/4/99 "
	^self class new: ( ( 0 to: (self degree max: aPolynomial degree)) collect: [ :n | ( aPolynomial at: n) - ( self at: n)])!

timesNumber: aNumber
		"(c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 17/4/99 "
	^self class new: ( coefficients reverse collect: [ :each | each * aNumber])!

timesPolynomial: aPolynomial
		"(c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 17/4/99 "
	| productCoefficients degree|
	degree := aPolynomial degree + self degree.
	productCoefficients := ( degree to: 0 by: -1)
			collect:[ :n | | sum |
					  sum := 0.
					  0 to: (degree - n)
						do: [ :k | sum := (self at: k) * (aPolynomial at: ( degree - n - k)) + sum].
					  sum
					].
	^self class new: productCoefficients!

value: aNumber
		"Answer the value of the polynomial for the specified variable value.
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 4/1/99 "
	^coefficients inject: 0 into: [ :sum :each | sum * aNumber + each]! !

!DhbPolynomial privateMethods !

initialize: anArray
		"Private - Initialize the coefficients of the receiver.
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 4/1/99 "
	coefficients := anArray.
	^self! !

!DhbSplineInterpolator publicMethods !

endPointDerivative: aNumber
		"Defines the end point derivatives.
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 7/4/99 "
	endPointDerivative := aNumber.
	self resetCoefficients.!

resetEndPointDerivatives
		"Set the end point derivatives to undefined.
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 7/4/99 "
	self setEndPointDerivatives: ( Array new: 2).!

setEndPointDerivatives: anArray
		"Defines the end point derivatives.
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 7/4/99 "
	startPointDerivative := anArray at: 1.
	endPointDerivative := anArray at: 2.
	self resetCoefficients.!

startPointDerivative: aNumber
		"Defines the end point derivatives.
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 7/4/99 "
	startPointDerivative := aNumber.
	self resetCoefficients.!

value: aNumber 
	"Computes the value of a cubic spline interpolation over the points of the receiver.
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 7/4/99 "

	| answer n1 n2 n step a b |
	coefficients isNil ifTrue: [self computeSecondDerivatives].
	n2 := pointCollection size.
	n1 := 1.
	[n2 - n1 > 1] whileTrue: 
			[n := (n1 + n2) // 2.
			(self xPointAt:  n) > aNumber ifTrue: [n2 := n] ifFalse: [n1 := n]].
	step := (self xPointAt: n2) - (self xPointAt: n1).
	a := ((self xPointAt: n2) - aNumber) / step.
	b := (aNumber - (self xPointAt: n1)) / step.
	^a * (self yPointAt: n1) + (b * (self yPointAt: n2)) 
		+ ((a * (a squared - 1) * (coefficients at: n1) 
				+ (b * (b squared - 1) * (coefficients at: n2))) * step squared 
				/ 6)! !

!DhbSplineInterpolator privateMethods !

computeSecondDerivatives
	"Private -
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 7/4/99 "

	| size u w s dx inv2dx |
	size := pointCollection size.
	coefficients := Array new: size.
	u := Array new: size - 1.
	startPointDerivative isNil 
		ifTrue: 
			[coefficients at: 1 put: 0.
			u at: 1 put: 0]
		ifFalse: 
			[coefficients at: 1 put: -1 / 2.
			s := 1 / (( self xPointAt: 2) x - ( self xPointAt: 1) x).
			u at: 1
				put: 3 * s 
						* (s * (( self yPointAt: size) - ( self yPointAt: size - 1)) 
								- startPointDerivative)].
	2 to: size - 1
		do: 
			[:n | 
			dx := (self xPointAt: n) - (self xPointAt: ( n - 1)).
			inv2dx := 1 / (( self xPointAt: n + 1) - (self xPointAt: n - 1)).
			s := dx * inv2dx.
			w := 1 / (s * (coefficients at: n - 1) + 2).
			coefficients at: n put: (s - 1) * w.
			u at: n
				put: (((( self yPointAt: n + 1) - ( self yPointAt: n)) 
						/ (( self xPointAt: n + 1) - ( self xPointAt: n)) 
							- ((( self yPointAt: n) - ( self yPointAt: n - 1)) / dx)) * 6 
						* inv2dx - ((u at: n - 1) * s)) 
						* w].
	endPointDerivative isNil 
		ifTrue: [coefficients at: size put: 0]
		ifFalse: 
			[w := 1 / 2.
			s := 1 / ((self xPointAt:  size) - (self xPointAt: ( size - 1))).
			u at: 1
				put: 3 * s * (endPointDerivative 
								- (s * (self yPointAt: size) - (self yPointAt: size - 1))).
			coefficients at: size
				put: s - (w * (u at: size - 1) / ((coefficients at: size - 1) * w + 1))].
	size - 1 to: 1
		by: -1
		do: 
			[:n | 
			coefficients at: n
				put: (coefficients at: n) * (coefficients at: n + 1) + (u at: n)]!

defaultSamplePoints
		"Private -
		 (c) Copyrights Didier BESSET, 1999, all rights reserved.
		 Initial code: 27/5/99 "
	^SortedCollection sortBlock: [ :a :b | a x < b x]! !

!Number publicMethods !

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

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

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

DhbLagrangeInterpolator initializeAfterLoad!

DhbNevilleInterpolator initializeAfterLoad!

DhbBulirschStoerInterpolator initializeAfterLoad!

DhbNewtonInterpolator initializeAfterLoad!

DhbSplineInterpolator initializeAfterLoad!

DhbPolynomial initializeAfterLoad!

DhbFunctionEvaluation initializeAfterLoad!