Fixed #221. The problem with complex sqrt

src/Math-Complex/PMComplex.class.st

@@ -403,6 +403,21 @@ PMComplex >> asComplex [
 	^self
 ]
 
+{ #category : #comparing }
+PMComplex >> closeTo: aNumber precision: aPrecision [
+	"A complex number self is close to a complex number other with precision epsilon if both
+	self real is close to other real with precision epsilon and
+	self imaginary is close to other imaginary with precision epsilon.
+
+	In case aNumber is real and not complex, we convert it to a complex number
+	other = aNumber + 0i"
+	| other |
+	other := aNumber asComplex.
+
+	^ (real closeTo: other real precision: aPrecision) and: [ 
+		imaginary closeTo: other imaginary precision: aPrecision ].
+]
+
 { #category : #arithmetic }
 PMComplex >> complexConjugate [
 
@@ -706,15 +721,18 @@ PMComplex >> sinh [
 
 { #category : #'mathematical functions' }
 PMComplex >> sqrt [
+	"Return the square root of the receiver with a positive imaginary part.
+	Implementation based on: https://en.wikipedia.org/wiki/Complex_number#Square_root"
 
-	"Return the square root of the receiver with a positive imaginary part."
-
-	| u v |
-	(imaginary = 0 and: [ real >= 0 ]) ifTrue: [ 
-		^ self class real: real sqrt imaginary: 0 ].
-	v := (self abs - real / 2) sqrt.
-	u := imaginary / 2 / v.
-	^ self class real: u imaginary: v
+	| newReal newImaginary imaginarySign |
+
+	imaginarySign := imaginary sign = 0
+		ifTrue: [ 1 ] ifFalse: [ imaginary sign ].
+
+	newReal := ((self abs + real) / 2) sqrt.
+	newImaginary := imaginarySign * ((self abs - real) / 2) sqrt.
+
+	^ newReal + newImaginary i
 ]
 
 { #category : #'mathematical functions' }

src/Math-Tests-Complex/PMComplexTest.class.st

@@ -178,6 +178,39 @@ PMComplexTest >> testBug1 [
 	self assert: (0.5 * (2 + 0 i) ln) exp equals: (0.5 * 2 ln) exp
 ]
 
+{ #category : #'testing - close to' }
+PMComplexTest >> testCloseTo [
+
+	self assert: 2 + 3.000000000000001i closeTo: 2 + 3i.
+	self assert: 2.000000000000001 + 3i closeTo: 2 + 3i.
+	self assert: 2.000000000000001 + 3.000000000000001i closeTo: 2 + 3i.
+]
+
+{ #category : #'testing - close to' }
+PMComplexTest >> testCloseToReal [
+
+	self assert: 2 + 0.000000000000001i closeTo: 2.
+	self assert: 2.000000000000001 + 0.000000000000001i closeTo: 2.
+
+	self assert: 2 + 0i closeTo: 2.000000000000001.
+	self deny: 2 + 3i closeTo: 2.000000000000001
+]
+
+{ #category : #'testing - close to' }
+PMComplexTest >> testCloseToRealWithPrecision [
+
+	self assert: 2 + 0.001i closeTo: 2 precision: 0.01.
+	self assert: 2.001 + 0.001i closeTo: 2 precision: 0.01.
+]
+
+{ #category : #'testing - close to' }
+PMComplexTest >> testCloseToWithPrecision [
+
+	self assert: 2 + 3.001i closeTo: 2 + 3i precision: 0.01.
+	self assert: 2.001 + 3i closeTo: 2 + 3i precision: 0.01.
+	self assert: 2.001 + 3.001i closeTo: 2 + 3i precision: 0.01.
+]
+
 { #category : #tests }
 PMComplexTest >> testComplexAsComplex [
 	| ineg |
@@ -511,6 +544,21 @@ PMComplexTest >> testNew [
 	self assert: c imaginary equals: 0
 ]
 
+{ #category : #'testing - close to' }
+PMComplexTest >> testNotCloseToRealWithPrecision [
+
+	self deny: 2 + 0.001i closeTo: 2 precision: 0.000001.
+	self deny: 2.001 + 0.001i closeTo: 2 precision: 0.000001.
+]
+
+{ #category : #'testing - close to' }
+PMComplexTest >> testNotCloseToWithPrecision [
+
+	self deny: 2 + 3.001i closeTo: 2 + 3i precision: 0.000001.
+	self deny: 2.001 + 3i closeTo: 2 + 3i precision: 0.000001.
+	self deny: 2.001 + 3.001i closeTo: 2 + 3i precision: 0.000001.
+]
+
 { #category : #tests }
 PMComplexTest >> testNumberAsComplex [
 	| minusOne |
@@ -685,7 +733,7 @@ PMComplexTest >> testSquareRootOfNegativePureImaginaryNumberIsAComplexNumberWith
 
 	squareRoot := pureImaginaryNumber sqrt.
 
-	expected := 2 sqrt negated + 2 sqrt i.
+	expected := 2 sqrt - 2 sqrt i.
 	self assert: squareRoot real closeTo: expected real.
 	self assert: squareRoot imaginary closeTo: expected imaginary
 
@@ -720,6 +768,44 @@ PMComplexTest >> testSquareRootOfPositiveRealNumberIsAComplexNumberWithOnlyAReal
 	self assert: squareRoot equals: expected
 ]
 
+{ #category : #tests }
+PMComplexTest >> testSquareRootOfVeryLargeRealAndVerySmallImaginary [
+	"This may cause problems because of the loss of precision.
+	Very large and very small floating point numbers have different units of least precision.
+
+	verySmall ~= 0.
+	veryLarge + verySmall = veryLarge."
+
+	| veryLarge verySmall complexNumber expected |
+
+	veryLarge := 1e20.
+	verySmall := 1e-20.
+
+	complexNumber := veryLarge + verySmall i.
+	expected := 1e10 + 1e-31 i.
+
+	self assert: complexNumber sqrt closeTo: expected
+]
+
+{ #category : #tests }
+PMComplexTest >> testSquareRootOfVerySmallRealAndVeryLargeImaginary [
+	"This may cause problems because of the loss of precision.
+	Very large and very small floating point numbers have different units of least precision.
+
+	verySmall ~= 0.
+	veryLarge + verySmall = veryLarge."
+
+	| veryLarge verySmall complexNumber expected |
+
+	veryLarge := 1e20.
+	verySmall := 1e-20.
+
+	complexNumber := verySmall + veryLarge i.
+	expected := 7071067811.865476 + 7071067811.865475 i.
+
+	self assert: complexNumber sqrt closeTo: expected
+]
+
 { #category : #tests }
 PMComplexTest >> testSquareRootOfZeroIsZero [
 	"comment stating purpose of instance-side method"