Changeset 14673

Timestamp:

Mar 11, 2011, 11:30:08 AM (14 years ago)

Author:

R. Matthew Emerson

Message:

When computing the absolute value of a complex number,
be a little more careful about it.

hypot in libm is even more careful, and one could
certainly make the argument that we should just arrange
to ff-call that.

This new version conses more than what we were
previously doing, but seems to be about the same speed.
If this ends up being a performance issue, writing a
LAP version probably would not be too tough.

Fixes ticket:830.

File:

• trunk/source/level-0/l0-numbers.lisp (modified) (2 diffs)

trunk/source/level-0/l0-numbers.lisp

@@ -915 +915 @@
          (sqrt (+ (* rx rx) (* ix ix))))
         (short-float
-         (%short-float (%hypot (%double-float rx)
-                               (%double-float ix))))
+         (%short-float (%double-float-hypot (%double-float rx)
+                                            (%double-float ix))))
         (double-float
-         (%hypot rx ix)))))))
+         (%double-float-hypot rx ix)))))))
 
 
@@ -2002 +2002 @@
   (/ 1 n))
 
-; x & y must both be double floats
-(defun %hypot (x y)
-  (with-stack-double-floats ((x**2) (y**2))
-    (let ((res**2 x**2))
-      (%double-float*-2! x x x**2)
-      (%double-float*-2! y y y**2)
-      (%double-float+-2! x**2 y**2 res**2)
-      (fsqrt res**2))))
-
-
+;; Compute (sqrt (+ (* x x) (* y y))), but
+;; try to be a little more careful about it.
+;; Both x and y must be double-floats.
+(defun %double-float-hypot (x y)
+  (with-stack-double-floats ((a) (b) (c))
+    (%%double-float-abs! x a)
+    (%%double-float-abs! y b)
+    (when (> a b)
+      (psetq a b b a))
+    (if (= b 0d0)
+      0d0
+      (progn
+        (%double-float/-2! a b c)
+        (* b (fsqrt (+ 1d0 (* c c))))))))
+