610 
;;;; 
;;;; 
611 
;;;; Internal functions: 
;;;; Internal functions: 
612 
;;;; square coercetocomplextype cssqs complexlogscaled 
;;;; square coercetocomplextype cssqs complexlogscaled 

;;;; suppressdividebyzero 

613 
;;;; 
;;;; 
614 
;;;; 
;;;; 
615 
;;;; Please send any bug reports, comments, or improvements to Raymond 
;;;; Please send any bug reports, comments, or improvements to Raymond 
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* ((fpmodes (getfloatingpointmodes)) 
(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 (doublefloat x y) 
(declare (doublefloat x y) 
698 
(type (doublefloat 0d0) rho) 
(type (doublefloat 0d0) rho) 
699 
(fixnum k)) 
(fixnum k)) 
700 
(unwindprotect 
;; 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 
(withfloattrapsmasked (:underflow :overflow) 
705 
;; see if any signal would have been raised. 
(setf rho (+ (square x) (square y))) 
706 

(cond ((and (or (floattrappingnanp rho) 
707 
(setfloatingpointmodes :traps nil :accruedexceptions nil) 
(floatinfinityp rho)) 
708 

(or (floatinfinityp (abs x)) 
709 
(setf rho (+ (square x) (square y))) 
(floatinfinityp (abs y)))) 
710 
(cond ((and (or (floattrappingnanp rho) 
(setf rho #.ext:doublefloatpositiveinfinity)) 
711 
(floatinfinityp rho)) 
((let ((threshold #.(/ leastpositivedoublefloat 
712 
(or (floatinfinityp (abs x)) 
doublefloatepsilon)) 
713 
(floatinfinityp (abs y)))) 
(traps (ldb vm::floatexceptionsbyte 
714 
(setf rho #.ext:doublefloatpositiveinfinity)) 
(vm:floatingpointmodes)))) 
715 
((let* ((exceptions (getf (getfloatingpointmodes) 
;; Overflow raised or (underflow raised and rho < 
716 
:accruedexceptions)) 
;; lambda/eps) 
717 
(threshold #.(/ leastpositivedoublefloat 
(or (not (zerop (logand vm:floatoverflowtrapbit traps))) 
718 
doublefloatepsilon))) 
(and (not (zerop (logand vm:floatunderflowtrapbit 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 #'setfloatingpointmodes fpmodes)))) 

724 


725 
(defun complexsqrt (z) 
(defun complexsqrt (z) 
726 
"Principle square root of Z 
"Principle square root of Z 
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 dividebyzero errors in the body. After 


;; execution of the body, the floating point modes are returned as 


;; they were before entry. 





(defmacro suppressdividebyzero (&body body) 


(let ((fpmodes (gensym))) 


`(let ((,fpmodes (getfloatingpointmodes))) 


(unwindprotect 


(progn 


;; Remove the dividebyzero trap but leave everything 


;; else undisturbed. 


(setfloatingpointmodes 


:traps (remove :dividebyzero (getf ,fpmodes :traps))) 


,@body) 


(apply #'setfloatingpointmodes ,fpmodes))))) 




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. 
928 
(declare (number z)) 
(declare (number z)) 
929 
(let ((sqrt1+z (complexsqrt (+ 1 z))) 
(let ((sqrt1+z (complexsqrt (+ 1 z))) 
930 
(sqrt1z (complexsqrt ( 1 z)))) 
(sqrt1z (complexsqrt ( 1 z)))) 
931 
(suppressdividebyzero 
(withfloattrapsmasked (:dividebyzero) 
932 
(complex (* 2 (atan (/ (realpart sqrt1z) 
(complex (* 2 (atan (/ (realpart sqrt1z) 
933 
(realpart sqrt1+z)))) 
(realpart sqrt1+z)))) 
934 
(asinh (imagpart (* (conjugate sqrt1+z) 
(asinh (imagpart (* (conjugate sqrt1+z) 
935 
sqrt1z))))))) 
sqrt1z))))))) 
936 


937 
(defun complexacosh (z) 
(defun complexacosh (z) 
938 
"Compute acosh z = 2 * log(sqrt((z+1)/2) + sqrt((z1)/2)) 
"Compute acosh z = 2 * log(sqrt((z+1)/2) + sqrt((z1)/2)) 
941 
(declare (number z)) 
(declare (number z)) 
942 
(let ((sqrtz1 (complexsqrt ( z 1))) 
(let ((sqrtz1 (complexsqrt ( z 1))) 
943 
(sqrtz+1 (complexsqrt (+ z 1)))) 
(sqrtz+1 (complexsqrt (+ z 1)))) 
944 
(suppressdividebyzero 
(withfloattrapsmasked (:dividebyzero) 
945 
(complex (asinh (realpart (* (conjugate sqrtz1) 
(complex (asinh (realpart (* (conjugate sqrtz1) 
946 
sqrtz+1))) 
sqrtz+1))) 
947 
(* 2 (atan (/ (imagpart sqrtz1) 
(* 2 (atan (/ (imagpart sqrtz1) 
948 
(realpart sqrtz+1)))))))) 
(realpart sqrtz+1)))))))) 
949 


950 


951 
(defun complexasin (z) 
(defun complexasin (z) 
955 
(declare (number z)) 
(declare (number z)) 
956 
(let ((sqrt1z (complexsqrt ( 1 z))) 
(let ((sqrt1z (complexsqrt ( 1 z))) 
957 
(sqrt1+z (complexsqrt (+ 1 z)))) 
(sqrt1+z (complexsqrt (+ 1 z)))) 
958 
(suppressdividebyzero 
(withfloattrapsmasked (:dividebyzero) 
959 
(complex (atan (/ (realpart z) 
(complex (atan (/ (realpart z) 
960 
(realpart (* sqrt1z sqrt1+z)))) 
(realpart (* sqrt1z sqrt1+z)))) 
961 
(asinh (imagpart (* (conjugate sqrt1z) 
(asinh (imagpart (* (conjugate sqrt1z) 
962 
sqrt1+z))))))) 
sqrt1+z))))))) 
963 


964 
(defun complexasinh (z) 
(defun complexasinh (z) 
965 
"Compute asinh z = log(z + sqrt(1 + z*z)) 
"Compute asinh z = log(z + sqrt(1 + z*z)) 