[squeak-dev] The Trunk: Kernel-nice.1459.mcz

commits at source.squeak.org commits at source.squeak.org
Sat Apr 23 19:41:25 UTC 2022


Nicolas Cellier uploaded a new version of Kernel to project The Trunk:
http://source.squeak.org/trunk/Kernel-nice.1459.mcz

==================== Summary ====================

Name: Kernel-nice.1459
Author: nice
Time: 23 April 2022, 9:40:29.771092 pm
UUID: 5aa845dd-8c95-4a0f-9b68-afdac84b11e7
Ancestors: Kernel-mt.1458

Fix some avoidable floating point overflow/underflow/divisionByZero in Complex>>#sqrt

This was reported at: https://github.com/PolyMathOrg/PolyMath/issues/221

=============== Diff against Kernel-mt.1458 ===============

Item was changed:
  ----- Method: Complex>>sqrt (in category 'mathematical functions') -----
  sqrt
+ 	"Return the square root of the receiver with a positive imaginary part.
+ 	Implementation notes:
+ 	the formulation used ensure a protection against floating point overflow/underflow.
+ 	it also result in a reasonable precision (around 3 ulp).
+ 	It is inspired by following reference, except that it uses pre-scaling rather than eception handling:
+ 	Implementing Complex Elementary Function Using Exception Handling
+ 	ACM Transactions on Mathematical Software - October 1994
+ 	Ping Tang and 3 other authors"
- 	"Return the square root of the receiver with a positive imaginary part."
  
+ 	| x y r s t scale |
+ 	real isZero
+ 		ifTrue:
+ 			[t := imaginary abs sqrt.
+ 			^self class real: t imaginary: (imaginary copySignTo: t)].
+ 	scale := real abs max: imaginary abs.
+ 	x := real / scale.
+ 	y := imaginary / scale.
+ 	r := (x squared + y squared) sqrt.
+ 	s := scale sqrt.
+ 	t := (r + x abs * 2) sqrt.
+ 	^real > 0
+ 		ifTrue: [self class real: t * s / 2 imaginary: y * s / t]
+ 		ifFalse: [self class real: y abs * s / t imaginary: (y copySignTo: t * s / 2)]!
- 	| 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!



More information about the Squeak-dev mailing list