(coerce (%sqrt (coerce number 'double-float)) 'single-float)))
(((foreach single-float double-float))
(if (minusp number)
- (complex-sqrt number)
+ (complex-sqrt (complex number))
(coerce (%sqrt (coerce number 'double-float))
'(dispatch-type number))))
((complex)
(((foreach single-float double-float))
(if (or (> number (coerce 1 '(dispatch-type number)))
(< number (coerce -1 '(dispatch-type number))))
- (complex-asin number)
+ (complex-asin (complex number))
(coerce (%asin (coerce number 'double-float))
'(dispatch-type number))))
((complex)
(((foreach single-float double-float))
(if (or (> number (coerce 1 '(dispatch-type number)))
(< number (coerce -1 '(dispatch-type number))))
- (complex-acos number)
+ (complex-acos (complex number))
(coerce (%acos (coerce number 'double-float))
'(dispatch-type number))))
((complex)
(coerce (%acosh (coerce number 'double-float)) 'single-float)))
(((foreach single-float double-float))
(if (< number (coerce 1 '(dispatch-type number)))
- (complex-acosh number)
+ (complex-acosh (complex number))
(coerce (%acosh (coerce number 'double-float))
'(dispatch-type number))))
((complex)
(((foreach single-float double-float))
(if (or (> number (coerce 1 '(dispatch-type number)))
(< number (coerce -1 '(dispatch-type number))))
- (complex-atanh number)
+ (complex-atanh (complex number))
(coerce (%atanh (coerce number 'double-float))
'(dispatch-type number))))
((complex)
;;; 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)
;;; principal square root of Z
;;;
-;;; Z may be any NUMBER, but the result is always a COMPLEX.
+;;; Z may be RATIONAL or COMPLEX; the result is always a COMPLEX.
(defun complex-sqrt (z)
- (declare (number 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)
;;;
;;; This is for use with J /= 0 only when |z| is huge.
(defun complex-log-scaled (z j)
- (declare (number z)
+ (declare (type (or rational complex) z)
(fixnum j))
;; The constants t0, t1, t2 should be evaluated to machine
;; precision. In addition, Kahan says the accuracy of log1p
;;;
;;; Z may be any number, but the result is always a complex.
(defun complex-log (z)
- (declare (number z))
+ (declare (type (or rational complex) z))
(complex-log-scaled z 0))
;;; KLUDGE: Let us note the following "strange" behavior. atanh 1.0d0
;;; i*y is never 0 since we have positive and negative zeroes. -- rtoy
;;; Compute atanh z = (log(1+z) - log(1-z))/2.
(defun complex-atanh (z)
- (declare (number z))
+ (declare (type (or rational complex) z))
(let* (;; constants
(theta (/ (sqrt most-positive-double-float) 4.0d0))
(rho (/ 4.0d0 (sqrt most-positive-double-float)))
;;; Compute tanh z = sinh z / cosh z.
(defun complex-tanh (z)
- (declare (number z))
+ (declare (type (or rational complex) z))
(let ((x (float (realpart z) 1.0d0))
(y (float (imagpart z) 1.0d0)))
(locally
;;
;; and these two expressions are equal if and only if arg conj z =
;; -arg z, which is clearly true for all z.
- (declare (number z))
+ (declare (type (or rational complex) z))
(let ((sqrt-1+z (complex-sqrt (+ 1 z)))
(sqrt-1-z (complex-sqrt (- 1 z))))
(with-float-traps-masked (:divide-by-zero)
;;;
;;; Z may be any NUMBER, but the result is always a COMPLEX.
(defun complex-acosh (z)
- (declare (number z))
+ (declare (type (or rational complex) z))
(let ((sqrt-z-1 (complex-sqrt (- z 1)))
(sqrt-z+1 (complex-sqrt (+ z 1))))
(with-float-traps-masked (:divide-by-zero)
;;;
;;; Z may be any NUMBER, but the result is always a COMPLEX.
(defun complex-asin (z)
- (declare (number z))
+ (declare (type (or rational complex) z))
(let ((sqrt-1-z (complex-sqrt (- 1 z)))
(sqrt-1+z (complex-sqrt (+ 1 z))))
(with-float-traps-masked (:divide-by-zero)
;;;
;;; Z may be any number, but the result is always a complex.
(defun complex-asinh (z)
- (declare (number z))
+ (declare (type (or rational complex) z))
;; asinh z = -i * asin (i*z)
(let* ((iz (complex (- (imagpart z)) (realpart z)))
(result (complex-asin iz)))
;;;
;;; Z may be any number, but the result is always a complex.
(defun complex-atan (z)
- (declare (number z))
+ (declare (type (or rational complex) z))
;; atan z = -i * atanh (i*z)
(let* ((iz (complex (- (imagpart z)) (realpart z)))
(result (complex-atanh iz)))
;;;
;;; Z may be any number, but the result is always a complex.
(defun complex-tan (z)
- (declare (number z))
+ (declare (type (or rational complex) z))
;; tan z = -i * tanh(i*z)
(let* ((iz (complex (- (imagpart z)) (realpart z)))
(result (complex-tanh iz)))
+++ /dev/null
-;;;; tests of irrational floating point functions
-
-;;;; This software is part of the SBCL system. See the README file for
-;;;; more information.
-;;;;
-;;;; While most of SBCL is derived from the CMU CL system, the test
-;;;; files (like this one) were written from scratch after the fork
-;;;; from CMU CL.
-;;;;
-;;;; This software is in the public domain and is provided with
-;;;; absolutely no warranty. See the COPYING and CREDITS files for
-;;;; more information.
-
-(in-package :cl-user)
-\f
-;;;; old bugs
-
-;;; This used to fail with
-;;; The value -0.44579905382680446d0 is not of type (DOUBLE-FLOAT 0.0d0).
-;;; MNA's port of Raymond Toy work on CMU CL fixed this in sbcl-0.6.12.53.
-(assert (equal (log #c(0.4 0.5)) #C(-0.44579905 0.8960554)))
-\f
-;;;; other tests
-
-;;; expt
-(assert (equal (expt #c(0 1) 2) -1))
-(assert (equal (prin1-to-string (expt 2 #c(0 1))) "#C(0.7692389 0.63896126)"))
-
-;;; log
-(assert (equal (prin1-to-string (log -3 10)) "#C(0.47712126 1.3643764)"))
-(assert (= (log 3 0) 0))
-
-;;; sqrt, isqrt
-(assert (= (sqrt 9) 3.0))
-(assert (= (sqrt -9.0) #c(0.0 3.0)))
-(assert (= (isqrt 9) 3))
-(assert (= (isqrt 26) 5))
-
-
-;;; sin, sinh, asin, asinh
-(assert (equal (prin1-to-string (sin (* 8 (/ pi 2)))) "-4.898425415289509d-16"))
-(assert (equal (prin1-to-string (sin (expt 10 3))) "0.82687956"))
-(assert (= (sinh 0) 0.0))
-(assert (equal (prin1-to-string (sinh #c(5.0 -9.6)))
- "#C(-73.06699 12.936809)"))
-(assert (= (sin (* #c(0 1) 5)) (* #c(0 1) (sinh 5))))
-(assert (= (sinh (* #c(0 1) 5)) (* #c(0 1) (sin 5))))
-(assert (equal (prin1-to-string (asin -1)) "-1.5707964"))
-(assert (= (asin 0) 0.0))
-(assert (= (asin 2) #c(1.5707964 -1.3169578)))
-(assert (equal (prin1-to-string (asinh 0.5)) "0.4812118"))
-(assert (equal (prin1-to-string (asinh 3/7)) "0.41643077"))
-
-;;; cos, cosh, acos, acosh
-(assert (= (cos 0) 1.0))
-(assert (equal (prin1-to-string (cos (/ pi 2))) "6.123031769111886d-17"))
-(assert (= (cosh 0) 1.0))
-(assert (equal (prin1-to-string (cosh 1)) "1.5430807"))
-(assert (= (cos (* #c(0 1) 5)) (cosh 5)))
-(assert (= (cosh (* #c(0 1) 5)) (cos 5)))
-(assert (equal (prin1-to-string (acos 0)) "1.5707964"))
-(assert (equal (prin1-to-string (acos -1)) "3.1415927"))
-(assert (equal (prin1-to-string (acos 2)) "#C(0.0 1.3169578)"))
-(assert (= (acos 1.00001) #c(0.0 0.0044751678)))
-(assert (= (acosh 0) #c(0 1.5707964)))
-(assert (= (acosh 1) 0))
-(assert (= (acosh -1) #c(0 3.1415927)))
-
-;;; tan, tanh
-(assert (equal (prin1-to-string (tan 1)) "1.5574077"))
-(assert (equal (prin1-to-string (tan (/ pi 2))) "1.6331778728383844d+16"))
-(assert (equal (prin1-to-string (tanh 0.00753)) "0.0075298576"))
-(assert (= (tanh 50) 1.0))
-(assert (= (tan (* #c(0 1) 5)) (* #c(0 1) (tanh 5))))
-(assert (= (atan 1) 0.7853982))
-(assert (equal (prin1-to-string (atanh 0.5) ) "0.54930615"))
-(assert (equal (prin1-to-string (atanh 3/7)) "0.45814538"))
projects / sbcl.git / commitdiff