Ignore:
Timestamp:
Mar 18, 2011, 9:19:07 PM (9 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:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/source/level-0/l0-numbers.lisp

    r14673 r14681  
    786786
    787787(defun /-2 (x y)
    788   (number-case x
    789     (double-float (number-case y
    790                     (double-float (/ (the double-float x) (the double-float y)))
    791                     (short-float (with-stack-double-floats ((dy y))
    792                                    (/ (the double-float x) (the double-float dy))))
    793                     (rational (dfloat-rat / x y))
    794                     (complex (let* ((ry (%realpart y))
    795                                     (iy (%imagpart y))
    796                                     (dn (+ (* ry ry) (* iy iy))))
    797                                (canonical-complex (/ (* x ry) dn) (/ (- (* x iy)) dn))))))
    798     (short-float (number-case y
    799                    (short-float (/ (the short-float x) (the short-float y)))
    800                    (double-float (with-stack-double-floats ((dx x))
    801                                    (/ (the double-float dx) (the double-float y))))
    802                    (rational (sfloat-rat / x y))
    803                    (complex (let* ((ry (%realpart y))
    804                                     (iy (%imagpart y))
    805                                     (dn (+ (* ry ry) (* iy iy))))
    806                                (canonical-complex (/ (* x ry) dn) (/ (- (* x iy)) dn))))))                   
    807     (integer (number-case y
    808                (double-float (rat-dfloat / x y))
    809                (short-float (rat-sfloat / x y))
    810                (integer (integer-/-integer x y))
    811                (complex (let* ((ry (%realpart y))
    812                                (iy (%imagpart y))
    813                                (dn (+ (* ry ry) (* iy iy))))
    814                           (canonical-complex (/ (* x ry) dn) (/ (- (* x iy)) dn))))
    815                (ratio
    816                 (if (eql 0 x)
    817                   0
    818                   (let* ((ny (%numerator y))
    819                          (dy (%denominator y))
    820                          (gcd (gcd x ny)))
    821                     (build-ratio (* (maybe-truncate x gcd) dy)
    822                                  (maybe-truncate ny gcd)))))))
    823     (complex (number-case y
    824                (complex (let* ((rx (%realpart x))
    825                                (ix (%imagpart x))
    826                                (ry (%realpart y))
    827                                (iy (%imagpart y))
    828                                (dn (+ (* ry ry) (* iy iy))))
    829                           (canonical-complex (/ (+ (* rx ry) (* ix iy)) dn)
    830                                              (/ (- (* ix ry) (* rx iy)) dn))))
    831                ((rational float)
    832                 (canonical-complex (/ (%realpart x) y) (/ (%imagpart x) y)))))
    833     (ratio (number-case y
    834              (double-float (rat-dfloat / x y))
    835              (short-float (rat-sfloat / x y))
    836              (integer
    837               (when (eql y 0)
    838                 (divide-by-zero-error '/ x y))
    839               (let* ((nx (%numerator x)) (gcd (gcd nx y)))
    840                 (build-ratio (maybe-truncate nx gcd)
    841                              (* (maybe-truncate y gcd) (%denominator x)))))
    842              (complex (let* ((ry (%realpart y))
    843                              (iy (%imagpart y))
    844                              (dn (+ (* ry ry) (* iy iy))))
    845                         (canonical-complex (/ (* x ry) dn) (/ (- (* x iy)) dn))))
    846              (ratio
    847               (let* ((nx (%numerator x))
    848                      (dx (%denominator x))
    849                      (ny (%numerator y))
    850                      (dy (%denominator y))
    851                      (g1 (gcd nx ny))
    852                      (g2 (gcd dx dy)))
    853                 (build-ratio (* (maybe-truncate nx g1) (maybe-truncate dy g2))
    854                              (* (maybe-truncate dx g2) (maybe-truncate ny g1)))))))))
     788  (macrolet ((real-complex-/ (x y)
     789               (let ((ry (gensym))
     790                     (iy (gensym))
     791                     (r (gensym))
     792                     (dn (gensym)))
     793                 `(let* ((,ry (%realpart ,y))
     794                         (,iy (%imagpart ,y)))
     795                    (if (> (abs ,ry) (abs ,iy))
     796                      (let* ((,r (/ ,iy ,ry))
     797                             (,dn (* ,ry (+ 1 (* ,r ,r)))))
     798                        (canonical-complex (/ ,x ,dn)
     799                                           (/ (- (* ,x ,r)) ,dn)))
     800                      (let* ((,r (/ ,ry ,iy))
     801                             (,dn (* ,iy (+ 1 (* ,r ,r)))))
     802                        (canonical-complex (/ (* ,x ,r) ,dn)
     803                                           (/ (- ,x) ,dn))))))))
     804    (number-case x
     805      (double-float (number-case y
     806                      (double-float (/ (the double-float x) (the double-float y)))
     807                      (short-float (with-stack-double-floats ((dy y))
     808                                     (/ (the double-float x)
     809                                        (the double-float dy))))
     810                      (rational (dfloat-rat / x y))
     811                      (complex (real-complex-/ x y))))
     812      (short-float (number-case y
     813                     (short-float (/ (the short-float x) (the short-float y)))
     814                     (double-float (with-stack-double-floats ((dx x))
     815                                     (/ (the double-float dx)
     816                                        (the double-float y))))
     817                     (rational (sfloat-rat / x y))
     818                     (complex (real-complex-/ x y))))
     819      (integer (number-case y
     820                 (double-float (rat-dfloat / x y))
     821                 (short-float (rat-sfloat / x y))
     822                 (integer (integer-/-integer x y))
     823                 (complex (real-complex-/ x y))
     824                 (ratio
     825                  (if (eql 0 x)
     826                    0
     827                    (let* ((ny (%numerator y))
     828                           (dy (%denominator y))
     829                           (gcd (gcd x ny)))
     830                      (build-ratio (* (maybe-truncate x gcd) dy)
     831                                   (maybe-truncate ny gcd)))))))
     832      (complex (number-case y
     833                 (complex (let* ((rx (%realpart x))
     834                                 (ix (%imagpart x))
     835                                 (ry (%realpart y))
     836                                 (iy (%imagpart y)))
     837                            (if (> (abs ry) (abs iy))
     838                              (let* ((r (/ iy ry))
     839                                     (dn (+ ry (* r iy))))
     840                                (canonical-complex (/ (+ rx (* ix r)) dn)
     841                                                   (/ (- ix (* rx r)) dn)))
     842                              (let* ((r (/ ry iy))
     843                                     (dn (+ iy (* r ry))))
     844                                (canonical-complex (/ (+ (* rx r) ix) dn)
     845                                                   (/ (- (* ix r) rx) dn))))))
     846                 ((rational float)
     847                  (canonical-complex (/ (%realpart x) y) (/ (%imagpart x) y)))))
     848      (ratio (number-case y
     849               (double-float (rat-dfloat / x y))
     850               (short-float (rat-sfloat / x y))
     851               (integer
     852                (when (eql y 0)
     853                  (divide-by-zero-error '/ x y))
     854                (let* ((nx (%numerator x)) (gcd (gcd nx y)))
     855                  (build-ratio (maybe-truncate nx gcd)
     856                               (* (maybe-truncate y gcd) (%denominator x)))))
     857               (complex (real-complex-/ x y))
     858               (ratio
     859                (let* ((nx (%numerator x))
     860                       (dx (%denominator x))
     861                       (ny (%numerator y))
     862                       (dy (%denominator y))
     863                       (g1 (gcd nx ny))
     864                       (g2 (gcd dx dy)))
     865                  (build-ratio (* (maybe-truncate nx g1)
     866                                  (maybe-truncate dy g2))
     867                               (* (maybe-truncate dx g2)
     868                                  (maybe-truncate ny g1))))))))))
    855869
    856870
Note: See TracChangeset for help on using the changeset viewer.