Changeset 14681 for trunk/source/level-0/l0-numbers.lisp

Timestamp:

Mar 18, 2011, 9:19:07 PM (10 years ago)

Author:

rme

Message:

In /-2, be more careful with division involving a complex operand.

The code was adapted from CMUCL. Fixes ticket:674.

File:

• trunk/source/level-0/l0-numbers.lisp (modified) (1 diff)

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

@@ -786 +786 @@
 
 (defun /-2 (x y)
-  (number-case x
-    (double-float (number-case y
-                    (double-float (/ (the double-float x) (the double-float y)))
-                    (short-float (with-stack-double-floats ((dy y))
-                                   (/ (the double-float x) (the double-float dy))))
-                    (rational (dfloat-rat / x y))
-                    (complex (let* ((ry (%realpart y))
-                                    (iy (%imagpart y))
-                                    (dn (+ (* ry ry) (* iy iy))))
-                               (canonical-complex (/ (* x ry) dn) (/ (- (* x iy)) dn))))))
-    (short-float (number-case y
-                   (short-float (/ (the short-float x) (the short-float y)))
-                   (double-float (with-stack-double-floats ((dx x))
-                                   (/ (the double-float dx) (the double-float y))))
-                   (rational (sfloat-rat / x y))
-                   (complex (let* ((ry (%realpart y))
-                                    (iy (%imagpart y))
-                                    (dn (+ (* ry ry) (* iy iy))))
-                               (canonical-complex (/ (* x ry) dn) (/ (- (* x iy)) dn))))))                   
-    (integer (number-case y
-               (double-float (rat-dfloat / x y))
-               (short-float (rat-sfloat / x y))
-               (integer (integer-/-integer x y))
-               (complex (let* ((ry (%realpart y))
-                               (iy (%imagpart y))
-                               (dn (+ (* ry ry) (* iy iy))))
-                          (canonical-complex (/ (* x ry) dn) (/ (- (* x iy)) dn))))
-               (ratio
-                (if (eql 0 x)
-                  0
-                  (let* ((ny (%numerator y)) 
-                         (dy (%denominator y)) 
-                         (gcd (gcd x ny)))
-                    (build-ratio (* (maybe-truncate x gcd) dy)
-                                 (maybe-truncate ny gcd)))))))
-    (complex (number-case y
-               (complex (let* ((rx (%realpart x))
-                               (ix (%imagpart x))
-                               (ry (%realpart y))
-                               (iy (%imagpart y))
-                               (dn (+ (* ry ry) (* iy iy))))
-                          (canonical-complex (/ (+ (* rx ry) (* ix iy)) dn)
-                                             (/ (- (* ix ry) (* rx iy)) dn))))
-               ((rational float)
-                (canonical-complex (/ (%realpart x) y) (/ (%imagpart x) y)))))
-    (ratio (number-case y
-             (double-float (rat-dfloat / x y))
-             (short-float (rat-sfloat / x y))
-             (integer
-              (when (eql y 0)
-                (divide-by-zero-error '/ x y))
-              (let* ((nx (%numerator x)) (gcd (gcd nx y)))
-                (build-ratio (maybe-truncate nx gcd)
-                             (* (maybe-truncate y gcd) (%denominator x)))))
-             (complex (let* ((ry (%realpart y))
-                             (iy (%imagpart y))
-                             (dn (+ (* ry ry) (* iy iy))))
-                        (canonical-complex (/ (* x ry) dn) (/ (- (* x iy)) dn))))
-             (ratio
-              (let* ((nx (%numerator x))
-                     (dx (%denominator x))
-                     (ny (%numerator y))
-                     (dy (%denominator y))
-                     (g1 (gcd nx ny))
-                     (g2 (gcd dx dy)))
-                (build-ratio (* (maybe-truncate nx g1) (maybe-truncate dy g2))
-                             (* (maybe-truncate dx g2) (maybe-truncate ny g1)))))))))
+  (macrolet ((real-complex-/ (x y)
+               (let ((ry (gensym))
+                     (iy (gensym))
+                     (r (gensym))
+                     (dn (gensym)))
+                 `(let* ((,ry (%realpart ,y))
+                         (,iy (%imagpart ,y)))
+                    (if (> (abs ,ry) (abs ,iy))
+                      (let* ((,r (/ ,iy ,ry))
+                             (,dn (* ,ry (+ 1 (* ,r ,r)))))
+                        (canonical-complex (/ ,x ,dn)
+                                           (/ (- (* ,x ,r)) ,dn)))
+                      (let* ((,r (/ ,ry ,iy))
+                             (,dn (* ,iy (+ 1 (* ,r ,r)))))
+                        (canonical-complex (/ (* ,x ,r) ,dn)
+                                           (/ (- ,x) ,dn))))))))
+    (number-case x
+      (double-float (number-case y
+                      (double-float (/ (the double-float x) (the double-float y)))
+                      (short-float (with-stack-double-floats ((dy y))
+                                     (/ (the double-float x)
+                                        (the double-float dy))))
+                      (rational (dfloat-rat / x y))
+                      (complex (real-complex-/ x y))))
+      (short-float (number-case y
+                     (short-float (/ (the short-float x) (the short-float y)))
+                     (double-float (with-stack-double-floats ((dx x))
+                                     (/ (the double-float dx)
+                                        (the double-float y))))
+                     (rational (sfloat-rat / x y))
+                     (complex (real-complex-/ x y))))
+      (integer (number-case y
+                 (double-float (rat-dfloat / x y))
+                 (short-float (rat-sfloat / x y))
+                 (integer (integer-/-integer x y))
+                 (complex (real-complex-/ x y))
+                 (ratio
+                  (if (eql 0 x)
+                    0
+                    (let* ((ny (%numerator y)) 
+                           (dy (%denominator y)) 
+                           (gcd (gcd x ny)))
+                      (build-ratio (* (maybe-truncate x gcd) dy)
+                                   (maybe-truncate ny gcd)))))))
+      (complex (number-case y
+                 (complex (let* ((rx (%realpart x))
+                                 (ix (%imagpart x))
+                                 (ry (%realpart y))
+                                 (iy (%imagpart y)))
+                            (if (> (abs ry) (abs iy))
+                              (let* ((r (/ iy ry))
+                                     (dn (+ ry (* r iy))))
+                                (canonical-complex (/ (+ rx (* ix r)) dn)
+                                                   (/ (- ix (* rx r)) dn)))
+                              (let* ((r (/ ry iy))
+                                     (dn (+ iy (* r ry))))
+                                (canonical-complex (/ (+ (* rx r) ix) dn)
+                                                   (/ (- (* ix r) rx) dn))))))
+                 ((rational float)
+                  (canonical-complex (/ (%realpart x) y) (/ (%imagpart x) y)))))
+      (ratio (number-case y
+               (double-float (rat-dfloat / x y))
+               (short-float (rat-sfloat / x y))
+               (integer
+                (when (eql y 0)
+                  (divide-by-zero-error '/ x y))
+                (let* ((nx (%numerator x)) (gcd (gcd nx y)))
+                  (build-ratio (maybe-truncate nx gcd)
+                               (* (maybe-truncate y gcd) (%denominator x)))))
+               (complex (real-complex-/ x y))
+               (ratio
+                (let* ((nx (%numerator x))
+                       (dx (%denominator x))
+                       (ny (%numerator y))
+                       (dy (%denominator y))
+                       (g1 (gcd nx ny))
+                       (g2 (gcd dx dy)))
+                  (build-ratio (* (maybe-truncate nx g1)
+                                  (maybe-truncate dy g2))
+                               (* (maybe-truncate dx g2)
+                                  (maybe-truncate ny g1))))))))))