[squeakdev] The Trunk: Kernelnice.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/Kernelnice.1459.mcz
==================== Summary ====================
Name: Kernelnice.1459
Author: nice
Time: 23 April 2022, 9:40:29.771092 pm
UUID: 5aa845dd8c954a0f9b68afdac84b11e7
Ancestors: Kernelmt.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 Kernelmt.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 prescaling 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 Squeakdev
mailing list
