complex arithmetic and mathematical functions

nicolas cellier ncellier at ifrance.com
Sun Feb 5 22:18:37 UTC 2006


Hello, i propose some extensions to Complex class.

Rationale:

First, a convenient creation message:

 (a i: b) is same as (a + b i)
 except it will not involve any operation, coercion etc... (Try operating on 
long collections of Complex)
 It produces code much shorter than (Complex real: a imaginary: b) found 
everywhere inside Complex class.
 Think of (a at b) and (Point x: a y: b).

Then:
 (aComplex i) has same meaning as (aNumber i) : it is multiply by (1 i).
 This one does exist in VW AT-Numerics for example and seems obvious.

Then, some functions arcCos arcSin arcTan argCosh argSinh argTanh arcTan: tanh

I also modified cos, sin, cosh and sinh only for having things like
 (2 + 0 i) sin imaginary isZero
Most definitions were taken from netlib SLATEC FORTRAN library.

I also propose these extensions in Number: sinh cosh tanh argSinh argCosh 
argTanh (some exist in Complex but not in Number...).
Further work: put some optional primitive in Float...

I also propose adding a new exception DomainError. This can enable trapping 
(-1 sqrt) and answering a complex result for example. (Further work to 
connect this exception if accepted)

Then a #absSecure like there is a #divideFastAndSecureBy: will avoid 
overflow/underflow in intermediate results.
Shouldn't these two become the default implementation ?

Then a #conjugated is missing.
Should we implement a (conjugated ^self) in Number ???

Last, raisedTo: and raisedToInteger: have to be implemented again in 
Complex... This will be the case for a lot of arithmetic classes. (Quaternion 
Matrix Polynomial...)
I see two better solutions
 1) make an ArithmeticValue common superclass like in VW
 2) use Traits to factor common protocol for commutative rings, fields, etc...
What are your plans about it ?

Once again, i am too long, bye.
-------------- next part --------------
'From Squeak3.9alpha of 4 July 2005 [latest update: #6719] on 5 February 2006 at 11:08:19 pm'!
ArithmeticError subclass: #DomainError
	instanceVariableNames: ''
	classVariableNames: ''
	poolDictionaries: ''
	category: 'Exceptions-Kernel'!

!DomainError commentStamp: 'nice 2/5/2006 21:45' prior: 0!
I am a kind of error that should be raised when a mathematical function is used outside its domain of definition.

For example, taking the arc sine of a number greater than one may raise me.

I then make extending of arc sine function definition over the complex field possible.
One just have to trap me, test the originator message and take appropriate action.!


!Complex methodsFor: 'arithmetic' stamp: 'nice 2/5/2006 22:26'!
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]
		ifFalse: [(real / scale i: imaginary / scale) squaredNorm sqrt * scale]! !

!Complex methodsFor: 'arithmetic' stamp: 'nice 2/5/2006 03:24'!
conjugated
	"Return the complex conjugate of this complex number."

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

!Complex methodsFor: 'arithmetic' stamp: 'nice 2/5/2006 03:25'!
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! !

!Complex methodsFor: 'arithmetic' stamp: 'nice 2/5/2006 22:24'!
squaredNorm
	"Answer the square of receiver norm."

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

!Complex methodsFor: 'mathematical functions' stamp: 'nice 2/5/2006 22:20'!
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: [self floatClass pi]
					ifFalse: [0].
				y := real abs argCosh.
				^x i: y]
			ifFalse: [^real arcCos]].
	tmp := self squaredNorm - 1 / 2.
	delta := tmp squared + imaginary squared.
	sh2y := tmp + delta sqrt.
	shy := sh2y sqrt.
	ch2y := 1 + sh2y.
	chy := ch2y sqrt.
	y := shy argSinh * imaginary sign negated.
	x := (real / chy) arcCos.
	^x i: y! !

!Complex methodsFor: 'mathematical functions' stamp: 'nice 2/5/2006 22:17'!
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 := self floatClass pi / 2 * real sign.
					y := real abs argCosh * real sign.
					^x i: y]
				ifFalse: [^real arcSin]].
	tmp := (self squaredNorm - 1) / 2.
	delta := tmp squared + imaginary squared.
	sh2y := tmp + delta sqrt.
	shy := sh2y sqrt.
	ch2y := 1 + sh2y.
	chy := ch2y sqrt.
	y := shy argSinh.
	x := (real / chy) arcSin.
	^imaginary > 0 ifTrue: [x i: y] ifFalse: [x i: y negated]! !

!Complex methodsFor: 'mathematical functions' stamp: 'nice 2/5/2006 22:18'!
arcTan
	"Answer the arc tangent of the receiver.
	This is the inverse function of tan."

	| r2 |
	r2 := self squaredNorm.
	^(real * 2 arcTan: 1 - r2) / 2
		i: ((r2 + (imaginary * 2) + 1) / (r2 - (imaginary * 2) + 1)) ln / 4! !

!Complex methodsFor: 'mathematical functions' stamp: 'nice 2/5/2006 22:18'!
arcTan: denominator 
	"Answer the  four quadrants arc tangent of receiver over denominator.
	This is equivalent to classical C/FORTRAN atan2(self,denominator)."

	^denominator isZero 
		ifTrue: 
			[self isZero 
				ifTrue: 
					["shouldn't it be an error ? ^DomainError signal: '0 arcTan: 0'"
					^0 i: 0]
				ifFalse: 
					[real < 0 
						ifTrue: [self floatClass pi / -2]
						ifFalse: [self floatClass pi / 2]]]
		ifFalse: 
			[| res |
			res := (self / denominator) arcTan.
			denominator real < 0 ifTrue: [res := res + self floatClass pi].
			res real > self floatClass pi 
				ifTrue: [res := res - (self floatClass pi * 2)].
			res]! !

!Complex methodsFor: 'mathematical functions' stamp: 'nice 2/5/2006 22:21'!
argCosh
	"Answer receiver's argument of hyperbolic cosine.
	That is the inverse function of cosh"

	"Some possible implementation:

	^imaginary > 0 
		ifTrue: [(self + (self * self - 1) sqrt) ln]
		ifFalse: [(self + (self * self - 1) sqrt) ln negated]"

	^self arcCos i negated! !

!Complex methodsFor: 'mathematical functions' stamp: 'nice 2/5/2006 22:21'!
argSinh
	"Answer receiver's argument of 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! !

!Complex methodsFor: 'mathematical functions' stamp: 'nice 2/5/2006 22:21'!
argTanh
	"Answer receiver's argument of hyperbolic tangent.
	That is the inverse function of tanh"

	"Some other possible implementation:

	^((1 + self) / (1 - self)) ln / 2"

	^self i arcTan i negated! !

!Complex methodsFor: 'mathematical functions' stamp: 'nice 2/5/2006 22:57'!
cos
	"Answer receiver's cosine."

	"Possible implementations:

	| iself |
	iself := 1 i * self.
	^ (iself exp + iself negated exp) / 2

	^self i exp - self i negated exp / 2i"

	^real cos * imaginary cosh i: real sin negated * imaginary sinh! !

!Complex methodsFor: 'mathematical functions' stamp: 'nice 2/5/2006 22:58'!
cosh
	"Answer receiver's hyperbolic cosine."

	"Possible implementations:

	^ (self exp + self negated exp) / 2
	
	^real cosh * imaginary cos i: real sinh * imaginary sin

	| ex emx |
	ex := real exp.
	emx := ex reciprocal.
	^(ex + emx) / 2 * imaginary cos i: (ex - emx) / 2 * imaginary sin"

	^self i cos! !

!Complex methodsFor: 'mathematical functions' stamp: 'nice 2/5/2006 22:22'!
raisedToInteger: operand 
	"Answer the receiver raised to the power operand, an Integer."

	"implementation note: this code is copied from Number.
	This suggest that both Number and Complex should have an
	ArithmeticValue common superclass like in Visualworks.
	Or maybe should it be a Traits (a property of fields ?)"

	| count result |
	#Numeric.
	"Changed 200/01/19 For ANSI <number> support."
	operand isInteger ifFalse: [^ ArithmeticError signal: 'parameter is not an Integer'"<- Chg"].
	operand = 0 ifTrue: [^ self class one].
	operand = 1 ifTrue: [^ self].
	operand < 0 ifTrue: [^ (self raisedToInteger: operand negated) reciprocal].
	count := 1.
	[(count := count + count) < operand] whileTrue.
	result := self class one.
	[count > 0]
		whileTrue: 
			[result := result * result.
			(operand bitAnd: count)
				= 0 ifFalse: [result := result * self].
			count := count bitShift: -1].
	^ result! !

!Complex methodsFor: 'mathematical functions' stamp: 'nice 2/5/2006 22:23'!
raisedTo: aNumber 
	"Answer the receiver raised to aNumber."

	aNumber isInteger ifTrue:
		["Do the special case of integer power"
		^ self raisedToInteger: aNumber].
	
	aNumber = 0 ifTrue: [^ 1].		"Special case of exponent=0 what if (0 raisedTo: 0) ? Should be an error"
	(self = 0) | (aNumber = 1) ifTrue:
		[^ self].						"Special case of exponent=1"
	^ (aNumber * self ln) exp		"Otherwise use logarithms"! !

!Complex methodsFor: 'mathematical functions' stamp: 'nice 2/5/2006 22:58'!
sin
	"Answer receiver's sine."

	"Possible implementations:

	| iself |
	iself := 1 i * self.
	^ (iself exp - iself negated exp) / 2 i

	^self i exp - self i negated exp / 2i"

	^real sin * imaginary cosh i: real cos * imaginary sinh! !

!Complex methodsFor: 'mathematical functions' stamp: 'nice 2/5/2006 22:58'!
sinh
	"Answer receiver's hyperbolic sine."

	"Possible implementations:

	^ (self exp - self negated exp) / 2

	^real sinh * imaginary cos i: real cosh * imaginary sin

	| ex emx |
	ex := real exp.
	emx := ex reciprocal.
	^(ex - emx) / 2 * imaginary cos i: (ex + emx) / 2 * imaginary sin"

	^self i sin i negated! !

!Complex methodsFor: 'mathematical functions' stamp: 'nice 2/5/2006 22:24'!
sqrt
	"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! !

!Complex methodsFor: 'mathematical functions' stamp: 'nice 2/5/2006 22:23'!
tanh
	"Answer receiver's hyperbolic tangent."

	"Some possible implementation are:

	^self sinh / self cosh

	| tr ti |
	tr := real tanh.
	ti := imaginary tan i.
	^(tr + ti) / (tr * ti + 1)"

	^self i tan i negated! !

!Complex methodsFor: 'private' stamp: 'nice 2/5/2006 03:24'!
floatClass
	"Answer the class suitable for doing floating point operations.
	In default Squeak, this is Float.
	In an image with single and double IEEE 754 floating point numbers,
	this would depend on the class of real and imaginary parts"

	^Float! !


!ComplexTest methodsFor: 'tests' stamp: 'nice 2/5/2006 22:50'!
testAbsSecure
	"self run: #testAbsSecure"
	"self debug: #testAbsSecure"
	
	self assert: (3  + 4 i) absSecure = 5.
	self should: [(3.0e200 + 4.0e200 i) absSecure = 5.0e200]
		description: 'An overflow should be avoided'.
	self 
		should: [(3.0e200 + 4.0e200 i) abs = Float infinity]
		description: 'An overflow can be avoided but is not'.
	self 
		should: 
			[| big |
			big := (2.0d raisedTo: 53) - 1 timesTwoPower: 1023 - 52.
			(big + big i) absSecure > big]
		raise: Error
		description: 'An overflow cannot be avoided'.
	self should: [(3.0e-200 + 4.0e-200 i) absSecure = 5.0e-200]
		description: 'An underflow should be avoided'.
	self should: [(3.0e-200 + 4.0e-200 i) abs = 0]
		description: 'An underflow can be avoided but is not'! !

!ComplexTest methodsFor: 'tests' stamp: 'nice 2/5/2006 22:43'!
testFunctions
	"self debug: #testFunctions"

	"Notice that current implementation is not accurate to machine epsilon...
	It is a naive implementation, an IEEE implementation should do better!!"

	| eps |
	eps := 1.0e-14.
	self assert: ((2 i: 0) sin - 2 sin) abs < eps.
	self assert: (1 / 2 i: 0) arcSin = (1 / 2) arcSin.
	self assert: ((2+ 3 i) arcSin sin - (2+ 3 i)) abs < eps.
	self assert: ((2 - 3 i) arcSin sin - (2 - 3 i)) abs < eps.
	self assert: ((-2+ 3 i) arcSin sin - (-2+ 3 i)) abs < eps.
	self assert: ((-2 - 3 i) arcSin sin - (-2 - 3 i)) abs < eps.

	self assert: ((2 i: 0) cos - 2 cos) abs < eps.
	self assert: (1 / 2 i: 0) arcCos = (1 / 2) arcCos.
	self assert: ((2+ 3 i) arcCos cos - (2+ 3 i)) abs < eps.
	self assert: ((2 - 3 i) arcCos cos - (2 - 3 i)) abs < eps.
	self assert: ((-2+ 3 i) arcCos cos - (-2+ 3 i)) abs < eps.
	self assert: ((-2 - 3 i) arcCos cos - (-2 - 3 i)) abs < eps.

	self assert: ((2 i: 0) tan - 2 tan) abs < eps.
	self assert: ((1 / 2 i: 0) arcTan - (1 / 2) arcTan) abs < eps.
	self assert: ((2+ 3 i) arcTan tan - (2+ 3 i)) abs < eps.
	self assert: ((2 - 3 i) arcTan tan - (2 - 3 i)) abs < eps.
	self assert: ((-2+ 3 i) arcTan tan - (-2+ 3 i)) abs < eps.
	self assert: ((-2 - 3 i) arcTan tan - (-2 - 3 i)) abs < eps.

	self assert: (2 i: 0) exp = 2 exp.
	self assert: (1 / 2 i: 0) ln = (1 / 2) ln.
	self assert: ((2+ 3 i) ln exp - (2+ 3 i)) abs < eps.
	self assert: ((2 - 3 i) ln exp - (2 - 3 i)) abs < eps.
	self assert: ((-2+ 3 i) ln exp - (-2+ 3 i)) abs < eps.
	self assert: ((-2 - 3 i) ln exp - (-2 - 3 i)) abs < eps.

	eps := 1.0e-14.

	self assert: ((2 + 0 i) sinh - 2 sinh) abs < eps.
	self assert: (1 / 2 + 0 i) argSinh = (1 / 2) argSinh.
	self assert: ((2+ 3 i) argSinh sinh - (2+ 3 i)) abs < eps.
	self assert: ((2 - 3 i) argSinh sinh - (2 - 3 i)) abs < eps.
	self assert: ((-2+ 3 i) argSinh sinh - (-2+ 3 i)) abs < eps.
	self assert: ((-2 - 3 i) argSinh sinh - (-2 - 3 i)) abs < eps.

	self assert: ((2 + 0 i) cosh - 2 cosh) abs < eps.
	self assert: ((2 + 0 i) argCosh - 2 argCosh) abs < eps.
	self assert: ((2+ 3 i) argCosh cosh - (2+ 3 i)) abs < eps.
	self assert: ((2 - 3 i) argCosh cosh - (2 - 3 i)) abs < eps.
	self assert: ((-2+ 3 i) argCosh cosh - (-2+ 3 i)) abs < eps.
	self assert: ((-2 - 3 i) argCosh cosh - (-2 - 3 i)) abs < eps.

	self assert: ((2 + 0 i) tanh - 2 tanh) abs < eps.
	self assert: ((1 / 2 + 0 i) argTanh - (1 / 2) argTanh) abs < eps.
	self assert: ((2+ 3 i) argTanh tanh - (2+ 3 i)) abs < eps.
	self assert: ((2 - 3 i) argTanh tanh - (2 - 3 i)) abs < eps.
	self assert: ((-2+ 3 i) argTanh tanh - (-2+ 3 i)) abs < eps.
	self assert: ((-2 - 3 i) argTanh tanh - (-2 - 3 i)) abs < eps.! !


!Number methodsFor: 'mathematical functions' stamp: 'nice 2/5/2006 22:10'!
argCosh
	"Answer receiver's argument of hyperbolic cosine.
	That is the inverse function of cosh"

	^self asFloat argCosh! !

!Number methodsFor: 'mathematical functions' stamp: 'nice 2/5/2006 22:10'!
argSinh
	"Answer receiver's argument of hyperbolic sine.
	That is the inverse function of sinh"

	^self asFloat argSinh! !

!Number methodsFor: 'mathematical functions' stamp: 'nice 2/5/2006 22:10'!
argTanh
	"Answer receiver's argument of hyperbolic tangent.
	That is the inverse function of tanh"

	^self asFloat argTanh! !

!Number methodsFor: 'mathematical functions' stamp: 'nice 2/5/2006 22:12'!
cosh
	"Answer receivers hyperbolic cosine."
	
	^self asFloat cosh! !

!Number methodsFor: 'mathematical functions' stamp: 'nice 2/5/2006 22:05'!
sinc
	"Answer receivers cardinal sine."
	
	^ self isZero
		ifTrue: [self class one]
		ifFalse: [self sin / self]! !

!Number methodsFor: 'mathematical functions' stamp: 'nice 2/5/2006 22:12'!
sinh
	"Answer receivers hyperbolic sine"
	
	^self asFloat sinh! !

!Number methodsFor: 'mathematical functions' stamp: 'nice 2/5/2006 21:54'!
tanh
	"Answer receivers hyperbolic tangent"
	
	^self asFloat tanh! !

!Number methodsFor: 'converting' stamp: 'nice 2/5/2006 03:31'!
i: aNumber
	"Form a complex number with
		receiver as realPart
		aNumber as imaginaryPart
	this is the same as (self + aNumber i) but a little bit more efficient."

	aNumber isNumber ifFalse: [self error: 'Badly formed complex number'].
	^Complex real: self imaginary: aNumber! !


!Float methodsFor: 'mathematical functions' stamp: 'nice 2/5/2006 22:09'!
argCosh
	"Answer receiver's argument of hyperbolic cosine.
	That is the inverse function of cosh"

	self < 1 
		ifTrue: 
			[^DomainError signal: 'Receiver must be greater or equal to 1'].
	^self + 1 = self 
		ifTrue: [self abs ln + 2 ln]
		ifFalse: [((self squared - 1) sqrt + self) ln]! !

!Float methodsFor: 'mathematical functions' stamp: 'nice 2/5/2006 22:09'!
argSinh
	"Answer receiver's argument of hyperbolic sine.
	That is the inverse function of sinh"

	^self + 1 = self 
		ifTrue: [self abs ln + 2 ln]
		ifFalse: [((self squared + 1) sqrt + self) ln]! !

!Float methodsFor: 'mathematical functions' stamp: 'nice 2/5/2006 22:09'!
argTanh
	"Answer receiver's argument of hyperbolic tangent.
	That is the inverse function of tanh"

	self abs >= 1 
		ifTrue: 
			[^DomainError signal: 'Receiver must be strictly between 1.0 and -1.0'].
	^((1 + self) / (1 - self)) ln / 2! !

!Float methodsFor: 'mathematical functions' stamp: 'nice 2/5/2006 22:11'!
cosh
	"Answer receivers hyperbolic cosine."
	
	| ex |
	ex := self exp.
	^(ex + ex reciprocal) / 2! !

!Float methodsFor: 'mathematical functions' stamp: 'nice 2/5/2006 22:11'!
sinh
	"Answer receivers hyperbolic sine"
	
	| ex |
	ex := self exp.
	^(ex - ex reciprocal) / 2! !

!Float methodsFor: 'mathematical functions' stamp: 'nice 2/5/2006 03:30'!
tanh
	"Answer hyperbolic tangent of receiver.
	Trivial implementation is:
		^self sinh/self cosh
	This implementation takes care not to overflow."

	| ex emx |
	self > 20.0 ifTrue: [^1.0].
	self < -20.0 ifTrue: [^-1.0].
	ex := self exp.
	emx := ex reciprocal.
	^(ex - emx) / (ex + emx)! !


!NumberTest methodsFor: 'tests' stamp: 'nice 2/5/2006 22:07'!
testFunctions
	"self debug: #testFunctions"
	
	| eps |
	eps := 1.0e-14.
	self assert: 0 sinc = 1.
	self assert: (Float pi / 2) sinc = (2 / Float pi).
	self assert: (Float pi) sinc abs < eps.
	
	self assert: 0 cosh = 1.
	self assert: 0 sinh = 0.
	self assert: 0 tanh = 0.
	self assert: (2 cosh argCosh - 2) abs < eps.
	self assert: (2 sinh argSinh - 2) abs < eps.
	self assert: (2 tanh argTanh - 2) abs < eps.! !



More information about the Squeak-dev mailing list