Skip to content
Snippets Groups Projects
Commit c275fc68 authored by toy's avatar toy
Browse files

o Fix the declaration bug in complex-log-scaled

o Remove the old special function routines
o Add logb-finite to help optimize the use of logb
o Remove some unneeded declarations since the compiler is smarter now
  than when this was originally written.
o Add inhibit-warnings to coerce-to-complex-type since their
  unavoidable.
o The cores of some routines are compiled with speed 3 and space 0 to
  get some maybe-inline routines inlined.
parent 2e346a66
No related branches found
No related tags found
No related merge requests found
...@@ -5,7 +5,7 @@ ...@@ -5,7 +5,7 @@
;;; Carnegie Mellon University, and has been placed in the public domain. ;;; Carnegie Mellon University, and has been placed in the public domain.
;;; ;;;
(ext:file-comment (ext:file-comment
"$Header: /Volumes/share2/src/cmucl/cvs2git/cvsroot/src/code/irrat.lisp,v 1.31 2001/03/04 20:12:37 pw Exp $") "$Header: /Volumes/share2/src/cmucl/cvs2git/cvsroot/src/code/irrat.lisp,v 1.32 2001/04/16 16:13:56 toy Exp $")
;;; ;;;
;;; ********************************************************************** ;;; **********************************************************************
;;; ;;;
...@@ -524,17 +524,6 @@ ...@@ -524,17 +524,6 @@
((complex) ((complex)
(complex-atan y))))) (complex-atan y)))))
;; It seems that everyone has a C version of sinh, cosh, and
;; tanh. Let's use these for reals because the original
;; implementations based on the definitions lose big in round-off
;; error. These bad definitions also mean that sin and cos for
;; complex numbers can also lose big.
#+nil
(defun sinh (number)
"Return the hyperbolic sine of NUMBER."
(/ (- (exp number) (exp (- number))) 2))
(defun sinh (number) (defun sinh (number)
"Return the hyperbolic sine of NUMBER." "Return the hyperbolic sine of NUMBER."
(number-dispatch ((number number)) (number-dispatch ((number number))
...@@ -545,11 +534,6 @@ ...@@ -545,11 +534,6 @@
(complex (* (sinh x) (cos y)) (complex (* (sinh x) (cos y))
(* (cosh x) (sin y))))))) (* (cosh x) (sin y)))))))
#+nil
(defun cosh (number)
"Return the hyperbolic cosine of NUMBER."
(/ (+ (exp number) (exp (- number))) 2))
(defun cosh (number) (defun cosh (number)
"Return the hyperbolic cosine of NUMBER." "Return the hyperbolic cosine of NUMBER."
(number-dispatch ((number number)) (number-dispatch ((number number))
...@@ -607,8 +591,9 @@ ...@@ -607,8 +591,9 @@
((complex) ((complex)
(complex-atanh number)))) (complex-atanh number))))
;;; HP-UX does not supply a C version of log1p, so ;;; HP-UX does not supply a C version of log1p, so use the definition.
;;; use the definition. ;;; We really need to fix this. The definition really loses big-time
;;; in roundoff as x gets small.
#+hpux #+hpux
(declaim (inline %log1p)) (declaim (inline %log1p))
...@@ -619,60 +604,6 @@ ...@@ -619,60 +604,6 @@
(the double-float (log (the (double-float 0d0) (+ number 1d0))))) (the double-float (log (the (double-float 0d0) (+ number 1d0)))))
#+old-elfun
(progn
;;; Here are the old definitions of the special functions, for
;;; complex-valued arguments. Some of these functions suffer from
;;; severe round-off error or unnecessary overflow.
(declaim (inline mult-by-i))
(defun mult-by-i (number)
(complex (- (imagpart number))
(realpart number)))
(defun complex-sqrt (x)
(exp (/ (log x) 2)))
(defun complex-log (x)
(complex (log (abs x))
(phase x)))
(defun complex-atanh (number)
(/ (- (log (1+ number)) (log (- 1 number))) 2))
(defun complex-tanh (number)
(/ (- (exp number) (exp (- number)))
(+ (exp number) (exp (- number)))))
(defun complex-acos (number)
(* -2 (mult-by-i (log (+ (sqrt (/ (1+ number) 2))
(mult-by-i (sqrt (/ (- 1 number) 2))))))))
(defun complex-acosh (number)
(* 2 (log (+ (sqrt (/ (1+ number) 2)) (sqrt (/ (1- number) 2))))))
(defun complex-asin (number)
(- (mult-by-i (log (+ (mult-by-i number) (sqrt (- 1 (* number number))))))))
(defun complex-asinh (number)
(log (+ number (sqrt (1+ (* number number))))))
(defun complex-atan (y)
(let ((im (imagpart y))
(re (realpart y)))
(/ (- (log (complex (- 1 im) re))
(log (complex (+ 1 im) (- re))))
(complex 0 2))))
(defun complex-tan (number)
(let* ((num (sin number))
(denom (cos number)))
(if (zerop denom) (error "~S undefined tangent." number)
(/ num denom))))
)
#-old-specfun
(progn
;;;; ;;;;
;;;; This is a set of routines that implement many elementary ;;;; This is a set of routines that implement many elementary
;;;; transcendental functions as specified by ANSI Common Lisp. The ;;;; transcendental functions as specified by ANSI Common Lisp. The
...@@ -723,11 +654,10 @@ ...@@ -723,11 +654,10 @@
;;;; State of the Art in Numerical Analysis", pp. 165-211, Clarendon ;;;; State of the Art in Numerical Analysis", pp. 165-211, Clarendon
;;;; Press, 1987 ;;;; Press, 1987
;;;; ;;;;
(declaim (inline square)) (declaim (inline square))
(declaim (ftype (function (double-float) (double-float 0d0)) square))
(defun square (x) (defun square (x)
(declare (double-float x) (declare (double-float x))
(values (double-float 0d0)))
(* x x)) (* x x))
;; If you have these functions in libm, perhaps they should be used ;; If you have these functions in libm, perhaps they should be used
...@@ -742,8 +672,20 @@ underlying floating-point format" ...@@ -742,8 +672,20 @@ underlying floating-point format"
(type double-float-exponent n)) (type double-float-exponent n))
(scale-float x n)) (scale-float x n))
(declaim (inline logb-finite))
(defun logb-finite (x)
"Same as logb but X is not infinity and non-zero and not a NaN, so
that we can always return an integer"
(declare (type double-float x))
(multiple-value-bind (signif expon sign)
(decode-float x)
(declare (ignore signif sign))
;; decode-float is almost right, except that the exponent
;; is off by one
(1- expon)))
(defun logb (x) (defun logb (x)
"Compute an integer N such that 1 <= |2^N * x| < 2. "Compute an integer N such that 1 <= |2^(-N) * x| < 2.
For the special cases, the following values are used: For the special cases, the following values are used:
x logb x logb
...@@ -758,17 +700,13 @@ For the special cases, the following values are used: ...@@ -758,17 +700,13 @@ For the special cases, the following values are used:
#.ext:double-float-positive-infinity) #.ext:double-float-positive-infinity)
((zerop x) ((zerop x)
;; The answer is negative infinity, but we are supposed to ;; The answer is negative infinity, but we are supposed to
;; signal divide-by-zero. ;; signal divide-by-zero, so do the actual division
;; (error 'division-by-zero :operation 'logb :operands (list x))
(/ -1.0d0 x) (/ -1.0d0 x)
) )
(t (t
(multiple-value-bind (signif expon sign) (logb-finite x))))
(decode-float x)
(declare (ignore signif sign))
;; decode-float is almost right, except that the exponent
;; is off by one
(1- expon)))))
;; This function is used to create a complex number of the appropriate ;; This function is used to create a complex number of the appropriate
;; type. ;; type.
...@@ -779,12 +717,13 @@ For the special cases, the following values are used: ...@@ -779,12 +717,13 @@ For the special cases, the following values are used:
it has the same type as Z. If Z has type (complex rational), the X it has the same type as Z. If Z has type (complex rational), the X
and Y are coerced to single-float." and Y are coerced to single-float."
(declare (double-float x y) (declare (double-float x y)
(number z)) (number z)
(if (subtypep (type-of (realpart z)) 'double-float) (optimize (extensions:inhibit-warnings 3)))
(if (typep (realpart z) 'double-float)
(complex x y) (complex x y)
;; Convert anything that's not a double-float to a single-float. ;; Convert anything that's not a double-float to a single-float.
(complex (float x 1.0) (complex (float x 1f0)
(float y 1.0)))) (float y 1f0))))
(defun cssqs (z) (defun cssqs (z)
;; Compute |(x+i*y)/2^k|^2 scaled to avoid over/underflow. The ;; Compute |(x+i*y)/2^k|^2 scaled to avoid over/underflow. The
...@@ -792,36 +731,37 @@ and Y are coerced to single-float." ...@@ -792,36 +731,37 @@ and Y are coerced to single-float."
;; Save all FP flags ;; Save all FP flags
(let ((x (float (realpart z) 1d0)) (let ((x (float (realpart z) 1d0))
(y (float (imagpart z) 1d0)) (y (float (imagpart z) 1d0)))
(k 0)
(rho 0d0))
(declare (double-float x y)
(type (double-float 0d0) rho)
(fixnum k))
;; Would this be better handled using an exception handler to ;; Would this be better handled using an exception handler to
;; catch the overflow or underflow signal? For now, we turn all ;; catch the overflow or underflow signal? For now, we turn all
;; traps off and look at the accrued exceptions to see if any ;; traps off and look at the accrued exceptions to see if any
;; signal would have been raised. ;; signal would have been raised.
(with-float-traps-masked (:underflow :overflow) (with-float-traps-masked (:underflow :overflow)
(setf rho (+ (square x) (square y))) (let ((rho (+ (square x) (square y))))
(cond ((and (or (float-nan-p rho) (declare (optimize (speed 3) (space 0)))
(float-infinity-p rho)) (cond ((and (or (float-nan-p rho)
(or (float-infinity-p (abs x)) (float-infinity-p rho))
(float-infinity-p (abs y)))) (or (float-infinity-p (abs x))
(setf rho #.ext:double-float-positive-infinity)) (float-infinity-p (abs y))))
((let ((threshold #.(/ least-positive-double-float (values ext:double-float-positive-infinity 0))
double-float-epsilon)) ((let ((threshold #.(/ least-positive-double-float
(traps (ldb vm::float-sticky-bits double-float-epsilon))
(vm:floating-point-modes)))) (traps (ldb vm::float-sticky-bits
;; Overflow raised or (underflow raised and rho < (vm:floating-point-modes))))
;; lambda/eps) ;; Overflow raised or (underflow raised and rho <
(or (not (zerop (logand vm:float-overflow-trap-bit traps))) ;; lambda/eps)
(and (not (zerop (logand vm:float-underflow-trap-bit traps))) (or (not (zerop (logand vm:float-overflow-trap-bit traps)))
(< rho threshold)))) (and (not (zerop (logand vm:float-underflow-trap-bit traps)))
(setf k (logb (max (abs x) (abs y)))) (< rho threshold))))
(setf rho (+ (square (scalb x (- k))) ;; If we're here, neither x nor y are infinity and at
(square (scalb y (- k)))))))) ;; least one is non-zero.. Thus logb returns a nice
(values rho k))) ;; integer.
(let ((k (- (logb-finite (max (abs x) (abs y))))))
(values (+ (square (scalb x k))
(square (scalb y k)))
(- k))))
(t
(values rho 0)))))))
(defun complex-sqrt (z) (defun complex-sqrt (z)
"Principle square root of Z "Principle square root of Z
...@@ -830,36 +770,40 @@ Z may be any number, but the result is always a complex." ...@@ -830,36 +770,40 @@ Z may be any number, but the result is always a complex."
(declare (number z)) (declare (number z))
(multiple-value-bind (rho k) (multiple-value-bind (rho k)
(cssqs z) (cssqs z)
(declare (type (double-float 0d0) rho) (declare (type (or (member 0d0) (double-float 0d0)) rho)
(fixnum k)) (type fixnum k))
(let ((x (float (realpart z) 1.0d0)) (let ((x (float (realpart z) 1.0d0))
(y (float (imagpart z) 1.0d0)) (y (float (imagpart z) 1.0d0))
(eta 0d0) (eta 0d0)
(nu 0d0)) (nu 0d0))
(declare (double-float x y eta nu)) (declare (double-float x y eta nu))
(if (not (float-nan-p x)) (locally
(setf rho (+ (scalb (abs x) (- k)) (sqrt rho)))) ;; space 0 to get maybe-inline functions inlined.
(declare (optimize (speed 3) (space 0)))
(if (not (float-nan-p x))
(setf rho (+ (scalb (abs x) (- k)) (sqrt rho))))
(cond ((oddp k)
(setf k (ash k -1)))
(t
(setf k (1- (ash k -1)))
(setf rho (+ rho rho))))
(setf rho (scalb (sqrt rho) k))
(setf eta rho)
(setf nu y)
(when (/= rho 0d0)
(when (not (float-infinity-p (abs nu)))
(setf nu (/ (/ nu rho) 2d0)))
(when (< x 0d0)
(setf eta (abs nu))
(setf nu (float-sign y rho))))
(coerce-to-complex-type eta nu z)))))
(cond ((oddp k)
(setf k (ash k -1)))
(t
(setf k (1- (ash k -1)))
(setf rho (+ rho rho))))
(setf rho (scalb (sqrt rho) k))
(setf eta rho)
(setf nu y)
(when (/= rho 0d0)
(when (not (float-infinity-p (abs nu)))
(setf nu (/ (/ nu rho) 2d0)))
(when (< x 0d0)
(setf eta (abs nu))
(setf nu (float-sign y rho))))
(coerce-to-complex-type eta nu z))))
(defun complex-log-scaled (z j) (defun complex-log-scaled (z j)
"Compute log(2^j*z). "Compute log(2^j*z).
...@@ -879,23 +823,21 @@ This is for use with J /= 0 only when |z| is huge." ...@@ -879,23 +823,21 @@ This is for use with J /= 0 only when |z| is huge."
(y (float (imagpart z) 1.0d0))) (y (float (imagpart z) 1.0d0)))
(multiple-value-bind (rho k) (multiple-value-bind (rho k)
(cssqs z) (cssqs z)
(declare (type (double-float 0d0) rho) (declare (optimize (speed 3)))
(fixnum k))
(let ((beta (max (abs x) (abs y))) (let ((beta (max (abs x) (abs y)))
(theta (min (abs x) (abs y)))) (theta (min (abs x) (abs y))))
(declare (type (double-float 0d0) beta theta)) (coerce-to-complex-type (if (and (zerop k)
(if (and (zerop k) (< t0 beta)
(< t0 beta) (or (<= beta t1)
(or (<= beta t1) (< rho t2)))
(< rho t2))) (/ (%log1p (+ (* (- beta 1.0d0)
(setf rho (/ (%log1p (+ (* (- beta 1.0d0) (+ beta 1.0d0))
(+ beta 1.0d0)) (* theta theta)))
(* theta theta))) 2d0)
2d0)) (+ (/ (log rho) 2d0)
(setf rho (+ (/ (log rho) 2d0) (* (+ k j) ln2)))
(* (+ k j) ln2)))) (atan y x)
(setf theta (atan y x)) z)))))
(coerce-to-complex-type rho theta z)))))
(defun complex-log (z) (defun complex-log (z)
"Log of Z = log |Z| + i * arg Z "Log of Z = log |Z| + i * arg Z
...@@ -913,92 +855,88 @@ Z may be any number, but the result is always a complex." ...@@ -913,92 +855,88 @@ Z may be any number, but the result is always a complex."
"Compute atanh z = (log(1+z) - log(1-z))/2" "Compute atanh z = (log(1+z) - log(1-z))/2"
(declare (number z)) (declare (number z))
(let* (;; Constants (let* (;; Constants
(theta #.(/ (sqrt most-positive-double-float) 4.0d0)) (theta (/ (sqrt most-positive-double-float) 4.0d0))
(rho #.(/ 4.0d0 (sqrt most-positive-double-float))) (rho (/ 4.0d0 (sqrt most-positive-double-float)))
(half-pi #.(/ pi 2.0d0)) (half-pi (/ pi 2.0d0))
(rp (float (realpart z) 1.0d0)) (rp (float (realpart z) 1.0d0))
(beta (float-sign rp 1.0d0)) (beta (float-sign rp 1.0d0))
(x (* beta rp)) (x (* beta rp))
(y (* beta (- (float (imagpart z) 1.0d0)))) (y (* beta (- (float (imagpart z) 1.0d0))))
(eta 0.0d0) (eta 0.0d0)
(nu 0.0d0)) (nu 0.0d0))
(declare (double-float theta rho half-pi rp beta y eta nu) ;; Shouldn't need this declare.
(type (double-float 0d0) x)) (declare (double-float x y))
(cond ((or (> x theta) (locally
(> (abs y) theta)) (declare (optimize (speed 3)))
;; To avoid overflow... (cond ((or (> x theta)
(setf eta (float-sign y half-pi)) (> (abs y) theta))
;; nu is real part of 1/(x + iy). This is x/(x^2+y^2), ;; To avoid overflow...
;; which can cause overflow. Arrange this computation so (setf eta (float-sign y half-pi))
;; that it won't overflow. ;; nu is real part of 1/(x + iy). This is x/(x^2+y^2),
(setf nu (let* ((x-bigger (> x (abs y))) ;; which can cause overflow. Arrange this computation so
(r (if x-bigger (/ y x) (/ x y))) ;; that it won't overflow.
(d (+ 1.0d0 (* r r)))) (setf nu (let* ((x-bigger (> x (abs y)))
(declare (double-float r d)) (r (if x-bigger (/ y x) (/ x y)))
(if x-bigger (d (+ 1.0d0 (* r r))))
(/ (/ x) d) (if x-bigger
(/ (/ r y) d))))) (/ (/ x) d)
((= x 1.0d0) (/ (/ r y) d)))))
;; Should this be changed so that if y is zero, eta is set ((= x 1.0d0)
;; to +infinity instead of approx 176? In any case ;; Should this be changed so that if y is zero, eta is set
;; tanh(176) is 1.0d0 within working precision. ;; to +infinity instead of approx 176? In any case
(let ((t1 (+ 4d0 (square y))) ;; tanh(176) is 1.0d0 within working precision.
(t2 (+ (abs y) rho))) (let ((t1 (+ 4d0 (square y)))
(declare (type (double-float 0d0) t1 t2)) (t2 (+ (abs y) rho)))
#+nil (setf eta (log (/ (sqrt (sqrt t1)))
(setf eta (log (/ (sqrt (sqrt t1))) (sqrt t2)))
(sqrt t2))) (setf nu (* 0.5d0
(setf eta (* 0.5d0 (log (the (double-float 0.0d0) (float-sign y
(/ (sqrt t1) t2))))) (+ half-pi (atan (* 0.5d0 t2))))))))
(setf nu (* 0.5d0 (t
(float-sign y (let ((t1 (+ (abs y) rho)))
(+ half-pi (atan (* 0.5d0 t2)))))))) ;; Normal case using log1p(x) = log(1 + x)
(t (setf eta (* 0.25d0
(let ((t1 (+ (abs y) rho))) (%log1p (/ (* 4.0d0 x)
(declare (double-float t1)) (+ (square (- 1.0d0 x))
;; Normal case using log1p(x) = log(1 + x) (square t1))))))
(setf eta (* 0.25d0 (setf nu (* 0.5d0
(%log1p (/ (* 4.0d0 x) (atan (* 2.0d0 y)
(+ (square (- 1.0d0 x)) (- (* (- 1.0d0 x)
(square t1)))))) (+ 1.0d0 x))
(setf nu (* 0.5d0 (square t1))))))))
(atan (* 2.0d0 y) (coerce-to-complex-type (* beta eta)
(- (* (- 1.0d0 x) (- (* beta nu))
(+ 1.0d0 x)) z))))
(square t1))))))))
(coerce-to-complex-type (* beta eta)
(- (* beta nu))
z)))
(defun complex-tanh (z) (defun complex-tanh (z)
"Compute tanh z = sinh z / cosh z" "Compute tanh z = sinh z / cosh z"
(declare (number z)) (declare (number z))
(let ((x (float (realpart z) 1.0d0)) (let ((x (float (realpart z) 1.0d0))
(y (float (imagpart z) 1.0d0))) (y (float (imagpart z) 1.0d0)))
(declare (double-float x y)) (locally
(cond ((> (abs x) ;; space 0 to get maybe-inline functions inlined
#-(or linux hpux) #.(/ (%asinh most-positive-double-float) 4d0) (declare (optimize (speed 3) (space 0)))
;; This is more accurate under linux. (cond ((> (abs x)
#+(or linux hpux) #.(/ (+ (%log 2.0d0) #-(or linux hpux) #.(/ (%asinh most-positive-double-float) 4d0)
(%log most-positive-double-float)) 4d0)) ;; This is more accurate under linux.
(complex (float-sign x) #+(or linux hpux) #.(/ (+ (%log 2.0d0)
(float-sign y 0.0d0))) (%log most-positive-double-float)) 4d0))
(t (coerce-to-complex-type (float-sign x)
(let* ((tv (%tan y)) (float-sign y) z))
(beta (+ 1.0d0 (* tv tv))) (t
(s (sinh x)) (let* ((tv (%tan y))
(rho (sqrt (+ 1.0d0 (* s s))))) (beta (+ 1.0d0 (* tv tv)))
(declare (double-float tv s) (s (sinh x))
(type (double-float 0.0d0) beta rho)) (rho (sqrt (+ 1.0d0 (* s s)))))
(if (float-infinity-p (abs tv)) (if (float-infinity-p (abs tv))
(coerce-to-complex-type (/ rho s) (coerce-to-complex-type (/ rho s)
(/ tv) (/ tv)
z) z)
(let ((den (+ 1.0d0 (* beta s s)))) (let ((den (+ 1.0d0 (* beta s s))))
(coerce-to-complex-type (/ (* beta rho s) (coerce-to-complex-type (/ (* beta rho s)
den) den)
(/ tv den) (/ tv den)
z)))))))) z)))))))))
;; Kahan says we should only compute the parts needed. Thus, the ;; Kahan says we should only compute the parts needed. Thus, the
;; realpart's below should only compute the real part, not the whole ;; realpart's below should only compute the real part, not the whole
...@@ -1094,4 +1032,4 @@ Z may be any number, but the result is always a complex." ...@@ -1094,4 +1032,4 @@ Z may be any number, but the result is always a complex."
(result (complex-tanh iz))) (result (complex-tanh iz)))
(complex (imagpart result) (complex (imagpart result)
(- (realpart result))))) (- (realpart result)))))
)
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment