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

Diff of /src/code/irrat.lisp

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

revision 1.20 by dtc, Tue Aug 26 18:00:24 1997 UTC revision 1.21 by dtc, Sat Aug 30 18:21:36 1997 UTC
# Line 610  Line 610 
610  ;;;;  ;;;;
611  ;;;; Internal functions:  ;;;; Internal functions:
612  ;;;;    square coerce-to-complex-type cssqs complex-log-scaled  ;;;;    square coerce-to-complex-type cssqs complex-log-scaled
 ;;;;    suppress-divide-by-zero  
613  ;;;;  ;;;;
614  ;;;;  ;;;;
615  ;;;; Please send any bug reports, comments, or improvements to Raymond  ;;;; Please send any bug reports, comments, or improvements to Raymond
# Line 687  and Y are coerced to single-float." Line 686  and Y are coerced to single-float."
686                 (float y 1.0))))                 (float y 1.0))))
687    
688  (defun cssqs (z)  (defun cssqs (z)
689    ;; Compute |x+i*y|^2/2^k carefully.  The result is r + i*k, where k    ;; Compute |(x+i*y)/2^k|^2 scaled to avoid over/underflow. The
690    ;; is an integer.    ;; result is r + i*k, where k is an integer.
691    
692    ;; Save all FP flags    ;; Save all FP flags
693    (let* ((fp-modes (get-floating-point-modes))    (let ((x (float (realpart z) 1d0))
694           (x (float (realpart z) 1d0))          (y (float (imagpart z) 1d0))
695           (y (float (imagpart z) 1d0))          (k 0)
696           (k 0)          (rho 0d0))
          (rho 0d0)  
          )  
697      (declare (double-float x y)      (declare (double-float x y)
698               (type (double-float 0d0) rho)               (type (double-float 0d0) rho)
699               (fixnum k))               (fixnum k))
700      (unwind-protect      ;; Would this be better handled using an exception handler to
701           (progn      ;; catch the overflow or underflow signal?  For now, we turn all
702             ;; Would this be better handled using an exception handler      ;; traps off and look at the accrued exceptions to see if any
703             ;; to catch the overflow or underflow signal?  For now, we      ;; signal would have been raised.
704             ;; turn all traps off and look at the accrued exceptions to      (with-float-traps-masked (:underflow :overflow)
705             ;; see if any signal would have been raised.        (setf rho (+ (square x) (square y)))
706          (cond ((and (or (float-trapping-nan-p rho)
707             (set-floating-point-modes :traps nil :accrued-exceptions nil)                        (float-infinity-p rho))
708                      (or (float-infinity-p (abs x))
709             (setf rho (+ (square x) (square y)))                        (float-infinity-p (abs y))))
710             (cond ((and (or (float-trapping-nan-p rho)               (setf rho #.ext:double-float-positive-infinity))
711                             (float-infinity-p rho))              ((let ((threshold #.(/ least-positive-double-float
712                         (or (float-infinity-p (abs x))                                     double-float-epsilon))
713                             (float-infinity-p (abs y))))                     (traps (ldb vm::float-exceptions-byte
714                    (setf rho #.ext:double-float-positive-infinity))                                 (vm:floating-point-modes))))
715                   ((let* ((exceptions (getf (get-floating-point-modes)                 ;; Overflow raised or (underflow raised and rho <
716                                             :accrued-exceptions))                 ;; lambda/eps)
717                           (threshold #.(/ least-positive-double-float                 (or (not (zerop (logand vm:float-overflow-trap-bit traps)))
718                                           double-float-epsilon)))                     (and (not (zerop (logand vm:float-underflow-trap-bit traps)))
719                      ;; overflow raised or (underflow raised and rho <                          (< rho threshold))))
720                      ;; lambda/eps)               (setf k (logb (max (abs x) (abs y))))
721                      (or (member :overflow exceptions)               (setf rho (+ (square (scalb x (- k)))
722                          (and (member :underflow exceptions)                            (square (scalb y (- k))))))))
723                               (< rho threshold))))      (values rho k)))
                   (setf k (logb (max (abs x) (abs y))))  
                   (setf rho (+ (square (scalb x (- k)))  
                                (square (scalb y (- k)))))))  
            (values rho k))  
       ;; restore over/underflow flags  
       (apply #'set-floating-point-modes fp-modes))))  
724    
725  (defun complex-sqrt (z)  (defun complex-sqrt (z)
726    "Principle square root of Z    "Principle square root of Z
# Line 916  Z may be any number, but the result is a Line 907  Z may be any number, but the result is a
907  ;; However, we take a pragmatic approach and just use the whole  ;; However, we take a pragmatic approach and just use the whole
908  ;; expression.  ;; expression.
909    
 ;; This macro suppresses any divide-by-zero errors in the body.  After  
 ;; execution of the body, the floating point modes are returned as  
 ;; they were before entry.  
   
 (defmacro suppress-divide-by-zero (&body body)  
   (let ((fp-modes (gensym)))  
     `(let ((,fp-modes (get-floating-point-modes)))  
       (unwind-protect  
            (progn  
              ;; Remove the divide-by-zero trap but leave everything  
              ;; else undisturbed.  
              (set-floating-point-modes  
               :traps (remove :divide-by-zero (getf ,fp-modes :traps)))  
              ,@body)  
         (apply #'set-floating-point-modes ,fp-modes)))))  
   
910  ;; NOTE: The formula given by Kahan is somewhat ambiguous in whether  ;; NOTE: The formula given by Kahan is somewhat ambiguous in whether
911  ;; it's the conjugate of the square root or the square root of the  ;; it's the conjugate of the square root or the square root of the
912  ;; conjugate.  This needs to be checked.  ;; conjugate.  This needs to be checked.
# Line 953  Z may be any number, but the result is a Line 928  Z may be any number, but the result is a
928    (declare (number z))    (declare (number z))
929    (let ((sqrt-1+z (complex-sqrt (+ 1 z)))    (let ((sqrt-1+z (complex-sqrt (+ 1 z)))
930          (sqrt-1-z (complex-sqrt (- 1 z))))          (sqrt-1-z (complex-sqrt (- 1 z))))
931      (suppress-divide-by-zero      (with-float-traps-masked (:divide-by-zero)
932       (complex (* 2 (atan (/ (realpart sqrt-1-z)        (complex (* 2 (atan (/ (realpart sqrt-1-z)
933                              (realpart sqrt-1+z))))                               (realpart sqrt-1+z))))
934                (asinh (imagpart (* (conjugate sqrt-1+z)                 (asinh (imagpart (* (conjugate sqrt-1+z)
935                                    sqrt-1-z)))))))                                     sqrt-1-z)))))))
936    
937  (defun complex-acosh (z)  (defun complex-acosh (z)
938    "Compute acosh z = 2 * log(sqrt((z+1)/2) + sqrt((z-1)/2))    "Compute acosh z = 2 * log(sqrt((z+1)/2) + sqrt((z-1)/2))
# Line 966  Z may be any number, but the result is a Line 941  Z may be any number, but the result is a
941    (declare (number z))    (declare (number z))
942    (let ((sqrt-z-1 (complex-sqrt (- z 1)))    (let ((sqrt-z-1 (complex-sqrt (- z 1)))
943          (sqrt-z+1 (complex-sqrt (+ z 1))))          (sqrt-z+1 (complex-sqrt (+ z 1))))
944      (suppress-divide-by-zero      (with-float-traps-masked (:divide-by-zero)
945       (complex (asinh (realpart (* (conjugate sqrt-z-1)        (complex (asinh (realpart (* (conjugate sqrt-z-1)
946                                    sqrt-z+1)))                                     sqrt-z+1)))
947                (* 2 (atan (/ (imagpart sqrt-z-1)                 (* 2 (atan (/ (imagpart sqrt-z-1)
948                              (realpart sqrt-z+1))))))))                               (realpart sqrt-z+1))))))))
949    
950    
951  (defun complex-asin (z)  (defun complex-asin (z)
# Line 980  Z may be any number, but the result is a Line 955  Z may be any number, but the result is a
955    (declare (number z))    (declare (number z))
956    (let ((sqrt-1-z (complex-sqrt (- 1 z)))    (let ((sqrt-1-z (complex-sqrt (- 1 z)))
957          (sqrt-1+z (complex-sqrt (+ 1 z))))          (sqrt-1+z (complex-sqrt (+ 1 z))))
958      (suppress-divide-by-zero      (with-float-traps-masked (:divide-by-zero)
959       (complex (atan (/ (realpart z)        (complex (atan (/ (realpart z)
960                         (realpart (* sqrt-1-z sqrt-1+z))))                          (realpart (* sqrt-1-z sqrt-1+z))))
961                (asinh (imagpart (* (conjugate sqrt-1-z)                 (asinh (imagpart (* (conjugate sqrt-1-z)
962                                    sqrt-1+z)))))))                                     sqrt-1+z)))))))
963    
964  (defun complex-asinh (z)  (defun complex-asinh (z)
965    "Compute asinh z = log(z + sqrt(1 + z*z))    "Compute asinh z = log(z + sqrt(1 + z*z))

Legend:
Removed from v.1.20  
changed lines
  Added in v.1.21

  ViewVC Help
Powered by ViewVC 1.1.5