+\f
+;;;; not-OLD-SPECFUN stuff
+;;;;
+;;;; (This was conditional on #-OLD-SPECFUN in the CMU CL sources,
+;;;; but OLD-SPECFUN was mentioned nowhere else, so it seems to be
+;;;; the standard special function system.)
+;;;;
+;;;; This is a set of routines that implement many elementary
+;;;; transcendental functions as specified by ANSI Common Lisp. The
+;;;; implementation is based on Kahan's paper.
+;;;;
+;;;; I believe I have accurately implemented the routines and are
+;;;; correct, but you may want to check for your self.
+;;;;
+;;;; These functions are written for CMU Lisp and take advantage of
+;;;; some of the features available there. It may be possible,
+;;;; however, to port this to other Lisps.
+;;;;
+;;;; Some functions are significantly more accurate than the original
+;;;; definitions in CMU Lisp. In fact, some functions in CMU Lisp
+;;;; give the wrong answer like (acos #c(-2.0 0.0)), where the true
+;;;; answer is pi + i*log(2-sqrt(3)).
+;;;;
+;;;; All of the implemented functions will take any number for an
+;;;; input, but the result will always be a either a complex
+;;;; single-float or a complex double-float.
+;;;;
+;;;; general functions:
+;;;; complex-sqrt
+;;;; complex-log
+;;;; complex-atanh
+;;;; complex-tanh
+;;;; complex-acos
+;;;; complex-acosh
+;;;; complex-asin
+;;;; complex-asinh
+;;;; complex-atan
+;;;; complex-tan
+;;;;
+;;;; utility functions:
+;;;; scalb logb
+;;;;
+;;;; internal functions:
+;;;; square coerce-to-complex-type cssqs complex-log-scaled
+;;;;
+;;;; references:
+;;;; Kahan, W. "Branch Cuts for Complex Elementary Functions, or Much
+;;;; Ado About Nothing's Sign Bit" in Iserles and Powell (eds.) "The
+;;;; State of the Art in Numerical Analysis", pp. 165-211, Clarendon
+;;;; Press, 1987
+;;;;
+;;;; The original CMU CL code requested:
+;;;; Please send any bug reports, comments, or improvements to
+;;;; Raymond Toy at <email address deleted during 2002 spam avalanche>.
+
+;;; FIXME: In SBCL, the floating point infinity constants like
+;;; SB!EXT:DOUBLE-FLOAT-POSITIVE-INFINITY aren't available as
+;;; constants at cross-compile time, because the cross-compilation
+;;; host might not have support for floating point infinities. Thus,
+;;; they're effectively implemented as special variable references,
+;;; and the code below which uses them might be unnecessarily
+;;; inefficient. Perhaps some sort of MAKE-LOAD-TIME-VALUE hackery
+;;; should be used instead? (KLUDGED 2004-03-08 CSR, by replacing the
+;;; special variable references with (probably equally slow)
+;;; constructors)
+;;;
+;;; FIXME: As of 2004-05, when PFD noted that IMAGPART and COMPLEX
+;;; differ in their interpretations of the real line, IMAGPART was
+;;; patch, which without a certain amount of effort would have altered
+;;; all the branch cut treatment. Clients of these COMPLEX- routines
+;;; were patched to use explicit COMPLEX, rather than implicitly
+;;; passing in real numbers for treatment with IMAGPART, and these
+;;; COMPLEX- functions altered to require arguments of type COMPLEX;
+;;; however, someone needs to go back to Kahan for the definitive
+;;; answer for treatment of negative real floating point numbers and
+;;; branch cuts. If adjustment is needed, it is probably the removal
+;;; of explicit calls to COMPLEX in the clients of irrational
+;;; functions. -- a slightly bitter CSR, 2004-05-16
+
+(declaim (inline square))
+(defun square (x)
+ (declare (double-float x))
+ (* x x))
+
+;;; original CMU CL comment, apparently re. SCALB and LOGB and
+;;; perhaps CSSQS:
+;;; If you have these functions in libm, perhaps they should be used
+;;; instead of these Lisp versions. These versions are probably good
+;;; enough, especially since they are portable.
+
+;;; Compute 2^N * X without computing 2^N first. (Use properties of
+;;; the underlying floating-point format.)
+(declaim (inline scalb))
+(defun scalb (x n)
+ (declare (type double-float x)
+ (type double-float-exponent n))
+ (scale-float x n))
+
+;;; This is like LOGB, but X is not infinity and non-zero and not a
+;;; NaN, so we can always return an integer.
+(declaim (inline logb-finite))
+(defun logb-finite (x)
+ (declare (type double-float x))
+ (multiple-value-bind (signif exponent sign)
+ (decode-float x)
+ (declare (ignore signif sign))
+ ;; DECODE-FLOAT is almost right, except that the exponent is off
+ ;; by one.
+ (1- exponent)))
+
+;;; Compute an integer N such that 1 <= |2^N * x| < 2.
+;;; For the special cases, the following values are used:
+;;; x logb
+;;; NaN NaN
+;;; +/- infinity +infinity
+;;; 0 -infinity
+(defun logb (x)
+ (declare (type double-float x))
+ (cond ((float-nan-p x)
+ x)
+ ((float-infinity-p x)
+ ;; DOUBLE-FLOAT-POSITIVE-INFINITY
+ (double-from-bits 0 (1+ sb!vm:double-float-normal-exponent-max) 0))
+ ((zerop x)
+ ;; The answer is negative infinity, but we are supposed to
+ ;; signal divide-by-zero, so do the actual division
+ (/ -1.0d0 x)
+ )
+ (t
+ (logb-finite x))))
+
+;;; This function is used to create a complex number of the
+;;; appropriate type:
+;;; Create complex number with real part X and imaginary part Y
+;;; such that has the same type as Z. If Z has type (complex
+;;; rational), the X and Y are coerced to single-float.
+#!+long-float (eval-when (:compile-toplevel :load-toplevel :execute)
+ (error "needs work for long float support"))
+(declaim (inline coerce-to-complex-type))
+(defun coerce-to-complex-type (x y z)
+ (declare (double-float x y)
+ (number z))
+ (if (typep (realpart z) 'double-float)
+ (complex x y)
+ ;; Convert anything that's not already a DOUBLE-FLOAT (because
+ ;; the initial argument was a (COMPLEX DOUBLE-FLOAT) and we
+ ;; haven't done anything to lose precision) to a SINGLE-FLOAT.
+ (complex (float x 1f0)
+ (float y 1f0))))
+
+;;; Compute |(x+i*y)/2^k|^2 scaled to avoid over/underflow. The
+;;; result is r + i*k, where k is an integer.
+#!+long-float (eval-when (:compile-toplevel :load-toplevel :execute)
+ (error "needs work for long float support"))
+(defun cssqs (z)
+ (let ((x (float (realpart z) 1d0))
+ (y (float (imagpart z) 1d0)))
+ ;; Would this be better handled using an exception handler to
+ ;; catch the overflow or underflow signal? For now, we turn all
+ ;; traps off and look at the accrued exceptions to see if any
+ ;; signal would have been raised.
+ (with-float-traps-masked (:underflow :overflow)
+ (let ((rho (+ (square x) (square y))))
+ (declare (optimize (speed 3) (space 0)))
+ (cond ((and (or (float-nan-p rho)
+ (float-infinity-p rho))
+ (or (float-infinity-p (abs x))
+ (float-infinity-p (abs y))))
+ ;; DOUBLE-FLOAT-POSITIVE-INFINITY
+ (values
+ (double-from-bits 0 (1+ sb!vm:double-float-normal-exponent-max) 0)
+ 0))
+ ((let ((threshold #.(/ least-positive-double-float
+ double-float-epsilon))
+ (traps (ldb sb!vm::float-sticky-bits
+ (sb!vm:floating-point-modes))))
+ ;; Overflow raised or (underflow raised and rho <
+ ;; lambda/eps)
+ (or (not (zerop (logand sb!vm:float-overflow-trap-bit traps)))
+ (and (not (zerop (logand sb!vm:float-underflow-trap-bit
+ traps)))
+ (< rho threshold))))
+ ;; If we're here, neither x nor y are infinity and at
+ ;; least one is non-zero.. Thus logb returns a nice
+ ;; integer.
+ (let ((k (- (logb-finite (max (abs x) (abs y))))))
+ (values (+ (square (scalb x k))
+ (square (scalb y k)))
+ (- k))))
+ (t
+ (values rho 0)))))))
+
+;;; principal square root of Z
+;;;
+;;; Z may be RATIONAL or COMPLEX; the result is always a COMPLEX.
+(defun complex-sqrt (z)
+ ;; KLUDGE: Here and below, we can't just declare Z to be of type
+ ;; COMPLEX, because one-arg COMPLEX on rationals returns a rational.
+ ;; Since there isn't a rational negative zero, this is OK from the
+ ;; point of view of getting the right answer in the face of branch
+ ;; cuts, but declarations of the form (OR RATIONAL COMPLEX) are
+ ;; still ugly. -- CSR, 2004-05-16
+ (declare (type (or complex rational) z))
+ (multiple-value-bind (rho k)
+ (cssqs z)
+ (declare (type (or (member 0d0) (double-float 0d0)) rho)
+ (type fixnum k))
+ (let ((x (float (realpart z) 1.0d0))
+ (y (float (imagpart z) 1.0d0))
+ (eta 0d0)
+ (nu 0d0))
+ (declare (double-float x y eta nu))
+
+ (locally
+ ;; 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)))))