0.8.10.27:

Two fixes for ugly specification regarding negative zeros
	... SXHASH is defined to respect similarity.  This is good in general,
		but bad in the presence of negative floating point zeros,
		which are similar to positive floating point zeros, and
		must therefore hash to the same value.  Make it so,
		courtesy of a neat trick (add 0.0) from
	... IMAGPART is specified to return (* 0 <real>) for reals.  This
		is different from what we were doing for negative floating
		point reals.  Make it so, but...
	... adjust the irrational function code so that floating point
		reals are treated identically to #c(<real> 0.0), so that
		we don't get a discontinuity in the real line.
	... adjust the hash.impure.lisp test to cope with the new caveat;
	... delete irrat.pure.lisp, because (a) no non-x86 platform ever
		passed it and (b) the IMAGPART changes have knock-on effects
		that cause x86 not to pass it either.  Replacing it with
		tests based on the IEEE floating point test vectors
		(available from netlib) would be a good thing.
This commit is contained in:
Christophe Rhodes 2004-05-17 09:44:45 +00:00
parent 2f1071f50a
commit b0a51fec91
7 changed files with 53 additions and 100 deletions

5
NEWS
View file

@ -2443,6 +2443,11 @@ changes in sbcl-0.8.11 relative to sbcl-0.8.10:
more likely. (thanks also to Bruno Haible for a similar report
and discussions)
** removed stack cleaning in the cleanup part of UNWIND-PROTECT.
** IMAGPART is specified (infelicitously) to return (* 0 <thing>)
for objects of type REAL. Make it so.
** SXHASH is specified (infelicitously) to respect similarity,
which means that (SXHASH 0.0) must equal (SXHASH -0.0). Make
it so.
planned incompatible changes in 0.8.x:
* (not done yet, but planned:) When the profiling interface settles

View file

@ -335,7 +335,7 @@
(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)
@ -431,7 +431,7 @@
(((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)
@ -448,7 +448,7 @@
(((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)
@ -539,7 +539,7 @@
(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)
@ -557,7 +557,7 @@
(((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)
@ -640,6 +640,19 @@
;;; 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)
@ -756,9 +769,15 @@
;;; 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)
@ -799,7 +818,7 @@
;;;
;;; 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
@ -834,7 +853,7 @@
;;;
;;; 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
@ -843,7 +862,7 @@
;;; 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)))
@ -900,7 +919,7 @@
;;; 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
@ -959,7 +978,7 @@
;;
;; 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)
@ -972,7 +991,7 @@
;;;
;;; 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)
@ -985,7 +1004,7 @@
;;;
;;; 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)
@ -998,7 +1017,7 @@
;;;
;;; 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)))
@ -1009,7 +1028,7 @@
;;;
;;; 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)))
@ -1020,7 +1039,7 @@
;;;
;;; 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)))

View file

@ -265,7 +265,7 @@
((complex rational)
(sb!kernel:%imagpart number))
(float
(float 0 number))
(* 0 number))
(t
0)))

View file

@ -17,14 +17,15 @@
;;; SXHASH of FLOAT values is defined directly in terms of DEFTRANSFORM in
;;; order to avoid boxing.
(deftransform sxhash ((x) (single-float))
'(let ((bits (single-float-bits x)))
'(let* ((val (+ 0.0f0 x))
(bits (single-float-bits val)))
(logxor 66194023
(sxhash (the fixnum
(logand most-positive-fixnum
(logxor bits
(ash bits -7))))))))
(deftransform sxhash ((x) (double-float))
'(let* ((val x)
'(let* ((val (+ 0.0d0 x))
(hi (double-float-high-bits val))
(lo (double-float-low-bits val))
(hilo (logxor hi lo)))

View file

@ -194,8 +194,13 @@
(unless (typep (funcall #'sxhash i) '(and fixnum unsigned-byte))
(error "bad SXHASH behavior for ~S" i))
(dolist (j sxhash-tests)
(unless (eq (t->boolean (equal i j))
(t->boolean (= (sxhash i) (sxhash j))))
(unless (or (eq (t->boolean (equal i j))
(t->boolean (= (sxhash i) (sxhash j))))
(and (typep i 'number)
(typep j 'number)
(= i j)
(subtypep (type-of i) (type-of j))
(subtypep (type-of j) (type-of i))))
;; (If you get a surprising failure here, maybe you were
;; just very unlucky; see the notes above.)
(error "bad SXHASH behavior for ~S ~S" i j))))

View file

@ -1,77 +0,0 @@
;;;; 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)
;;;; 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)))
;;;; 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"))

View file

@ -17,4 +17,4 @@
;;; checkins which aren't released. (And occasionally for internal
;;; versions, especially for internal versions off the main CVS
;;; branch, it gets hairier, e.g. "0.pre7.14.flaky4.13".)
"0.8.10.26"
"0.8.10.27"