/[cmucl]/src/code/irrat-dd.lisp
ViewVC logotype

Diff of /src/code/irrat-dd.lisp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.9 by rtoy, Wed May 23 16:48:50 2007 UTC revision 1.10 by rtoy, Thu May 24 19:13:17 2007 UTC
# Line 884  Line 884 
884    (declare (type double-double-float x y)    (declare (type double-double-float x y)
885             (optimize (speed 3) (space 0)))             (optimize (speed 3) (space 0)))
886    (let ((code 0)    (let ((code 0)
887          (w 0w0))          (neg-x (minusp (float-sign x)))
888      (declare (type (integer 0 3) code)          (neg-y (minusp (float-sign y))))
889               (type double-double-float w))      (declare (type (integer 0 3) code))
890      (when (minusp x)      (when neg-x
891        (setf code 2))        (setf code 2))
892      (when (minusp y)      (when neg-y
893        (setf code (logior code 1)))        (setf code (logior code 1)))
894      (when (zerop x)      ;; Handle the case where x or y is zero, taking into account the
895        (unless (zerop (logand code 1))      ;; sign of the zero.
896          (return-from dd-%atan2 (- dd-pi/2)))      (cond ((zerop x)
897        (when (zerop y)             ;; x = 0.
898          (return-from dd-%atan2 0w0))             (cond ((zerop y)
899        (return-from dd-%atan2 dd-pi/2))                    ;; x = 0, y = 0
900      (when (zerop y)                    ;;
901        (return-from dd-%atan2                    ;; CLHS Figure 12-15 shows what the results should
902          (if (zerop (logand code 2))                    ;; be.
903              0w0                    (if neg-x
904              dd-pi)))                        (if neg-y
905      (setf w (ecase code                            (- dd-pi)
906                (0 0w0)                            dd-pi)
907                (1 0w0)                        (if neg-y
908                (2 dd-pi)                            -0w0
909                (3 (- dd-pi))))                            0w0)))
910                     (t
911                      ;; x = 0, y /= 0
912                      (if neg-y
913                          (- dd-pi/2)
914                          dd-pi/2))))
915              ((zerop y)
916               ;; y = 0, x /= 0
917               (cond (neg-x
918                      (if neg-y
919                          (- dd-pi)
920                          dd-pi))
921                     (t
922                      (if neg-y
923                          -0w0
924                          0w0))))
925              (t
926               (let ((w (ecase code
927                          (0
928                           ;; x > 0, y > 0
929                           0w0)
930                          (1
931                           ;; x > 0, y < 0
932                           -0w0)
933                          (2
934                           ;; x < 0, y > 0
935                           dd-pi)
936                          (3
937                           ;; x < 0, y < 0
938                           (- dd-pi)))))
939    
940      (+ w (dd-%atan (/ y x)))))               (+ w (dd-%atan (/ y x))))))))
941    
942    
943  ;;  ;;
# Line 1723  Z may be any number, but the result is a Line 1752  Z may be any number, but the result is a
1752                  (complex (* -0.5w0 (dd-%log1p (/ 2 (- (abs z) 1))))                  (complex (* -0.5w0 (dd-%log1p (/ 2 (- (abs z) 1))))
1753                           (- dd-pi/2)))))                           (- dd-pi/2)))))
1754          (t          (t
1755             (flet ((careful-mul (a b)
1756                      ;; Carefully multiply a and b, taking care to handle
1757                      ;; signed zeroes.  Only need to handle the case of b
1758                      ;; being zero.
1759                      (if (zerop b)
1760                          (if (minusp (* (float-sign a) (float-sign b)))
1761                              -0w0
1762                              0w0)
1763                          (* a b))))
1764           (let* ( ;; Constants           (let* ( ;; Constants
1765                  (theta (/ (sqrt most-positive-double-float) 4.0w0))                  (theta (/ (sqrt most-positive-double-float) 4.0w0))
1766                  (rho (/ 4.0w0 (sqrt most-positive-double-float)))                  (rho (/ 4.0w0 (sqrt most-positive-double-float)))
# Line 1730  Z may be any number, but the result is a Line 1768  Z may be any number, but the result is a
1768                  (rp (float (realpart z) 1.0w0))                  (rp (float (realpart z) 1.0w0))
1769                  (beta (float-sign rp 1.0w0))                  (beta (float-sign rp 1.0w0))
1770                  (x (* beta rp))                  (x (* beta rp))
1771                  (y (* beta (- (float (imagpart z) 1.0w0))))                  (y (careful-mul beta (- (float (imagpart z) 1.0w0))))
1772                  (eta 0.0w0)                  (eta 0.0w0)
1773                  (nu 0.0w0))                  (nu 0.0w0))
1774             ;; Shouldn't need this declare.             ;; Shouldn't need this declare.
# Line 1769  Z may be any number, but the result is a Line 1807  Z may be any number, but the result is a
1807                                                   (+ (square (- 1.0d0 x))                                                   (+ (square (- 1.0d0 x))
1808                                                      (square t1))))))                                                      (square t1))))))
1809                        (setf nu (* 0.5d0                        (setf nu (* 0.5d0
1810                                    (dd-%atan2 (* 2.0d0 y)                                    (dd-%atan2 (careful-mul 2.0d0 y)
1811                                               (- (* (- 1.0d0 x)                                               (- (* (- 1.0d0 x)
1812                                                     (+ 1.0d0 x))                                                     (+ 1.0d0 x))
1813                                                  (square t1))))))))                                                  (square t1))))))))
1814               (complex (* beta eta)               (complex (* beta eta)
1815                        (- (* beta nu))))))))                        (- (* beta nu)))))))))
1816    
1817  (defun dd-complex-tanh (z)  (defun dd-complex-tanh (z)
1818    "Compute tanh z = sinh z / cosh z"    "Compute tanh z = sinh z / cosh z"

Legend:
Removed from v.1.9  
changed lines
  Added in v.1.10

  ViewVC Help
Powered by ViewVC 1.1.5