From: Christophe Rhodes <crhodes@us...>  20040517 09:44:56

Update of /cvsroot/sbcl/sbcl/src/code In directory sc8prcvs1.sourceforge.net:/tmp/cvsserv5116/src/code Modified Files: numbers.lisp irrat.lisp sxhash.lisp Log Message: 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 nonx86 platform ever passed it and (b) the IMAGPART changes have knockon 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. Index: numbers.lisp =================================================================== RCS file: /cvsroot/sbcl/sbcl/src/code/numbers.lisp,v retrieving revision 1.34 retrieving revision 1.35 diff u d r1.34 r1.35  numbers.lisp 3 May 2004 23:01:28 0000 1.34 +++ numbers.lisp 17 May 2004 09:44:45 0000 1.35 @@ 265,7 +265,7 @@ ((complex rational) (sb!kernel:%imagpart number)) (float  (float 0 number)) + (* 0 number)) (t 0))) Index: irrat.lisp =================================================================== RCS file: /cvsroot/sbcl/sbcl/src/code/irrat.lisp,v retrieving revision 1.23 retrieving revision 1.24 diff u d r1.23 r1.24  irrat.lisp 9 Mar 2004 12:08:40 0000 1.23 +++ irrat.lisp 17 May 2004 09:44:45 0000 1.24 @@ 335,7 +335,7 @@ (coerce (%sqrt (coerce number 'doublefloat)) 'singlefloat))) (((foreach singlefloat doublefloat)) (if (minusp number)  (complexsqrt number) + (complexsqrt (complex number)) (coerce (%sqrt (coerce number 'doublefloat)) '(dispatchtype number)))) ((complex) @@ 431,7 +431,7 @@ (((foreach singlefloat doublefloat)) (if (or (> number (coerce 1 '(dispatchtype number))) (< number (coerce 1 '(dispatchtype number))))  (complexasin number) + (complexasin (complex number)) (coerce (%asin (coerce number 'doublefloat)) '(dispatchtype number)))) ((complex) @@ 448,7 +448,7 @@ (((foreach singlefloat doublefloat)) (if (or (> number (coerce 1 '(dispatchtype number))) (< number (coerce 1 '(dispatchtype number))))  (complexacos number) + (complexacos (complex number)) (coerce (%acos (coerce number 'doublefloat)) '(dispatchtype number)))) ((complex) @@ 539,7 +539,7 @@ (coerce (%acosh (coerce number 'doublefloat)) 'singlefloat))) (((foreach singlefloat doublefloat)) (if (< number (coerce 1 '(dispatchtype number)))  (complexacosh number) + (complexacosh (complex number)) (coerce (%acosh (coerce number 'doublefloat)) '(dispatchtype number)))) ((complex) @@ 557,7 +557,7 @@ (((foreach singlefloat doublefloat)) (if (or (> number (coerce 1 '(dispatchtype number))) (< number (coerce 1 '(dispatchtype number))))  (complexatanh number) + (complexatanh (complex number)) (coerce (%atanh (coerce number 'doublefloat)) '(dispatchtype number)))) ((complex) @@ 640,6 +640,19 @@ ;;; should be used instead? (KLUDGED 20040308 CSR, by replacing the ;;; special variable references with (probably equally slow) ;;; constructors) +;;; +;;; FIXME: As of 200405, 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, 20040516 (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 complexsqrt (z)  (declare (number z)) + ;; KLUDGE: Here and below, we can't just declare Z to be of type + ;; COMPLEX, because onearg 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, 20040516 + (declare (type (or complex rational) z)) (multiplevaluebind (rho k) (cssqs z) (declare (type (or (member 0d0) (doublefloat 0d0)) rho) @@ 799,7 +818,7 @@ ;;; ;;; This is for use with J /= 0 only when z is huge. (defun complexlogscaled (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 complexlog (z)  (declare (number z)) + (declare (type (or rational complex) z)) (complexlogscaled 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(1z))/2. (defun complexatanh (z)  (declare (number z)) + (declare (type (or rational complex) z)) (let* (;; constants (theta (/ (sqrt mostpositivedoublefloat) 4.0d0)) (rho (/ 4.0d0 (sqrt mostpositivedoublefloat))) @@ 900,7 +919,7 @@ ;;; Compute tanh z = sinh z / cosh z. (defun complextanh (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 ((sqrt1+z (complexsqrt (+ 1 z))) (sqrt1z (complexsqrt ( 1 z)))) (withfloattrapsmasked (:dividebyzero) @@ 972,7 +991,7 @@ ;;; ;;; Z may be any NUMBER, but the result is always a COMPLEX. (defun complexacosh (z)  (declare (number z)) + (declare (type (or rational complex) z)) (let ((sqrtz1 (complexsqrt ( z 1))) (sqrtz+1 (complexsqrt (+ z 1)))) (withfloattrapsmasked (:dividebyzero) @@ 985,7 +1004,7 @@ ;;; ;;; Z may be any NUMBER, but the result is always a COMPLEX. (defun complexasin (z)  (declare (number z)) + (declare (type (or rational complex) z)) (let ((sqrt1z (complexsqrt ( 1 z))) (sqrt1+z (complexsqrt (+ 1 z)))) (withfloattrapsmasked (:dividebyzero) @@ 998,7 +1017,7 @@ ;;; ;;; Z may be any number, but the result is always a complex. (defun complexasinh (z)  (declare (number z)) + (declare (type (or rational complex) z)) ;; asinh z = i * asin (i*z) (let* ((iz (complex ( (imagpart z)) (realpart z))) (result (complexasin iz))) @@ 1009,7 +1028,7 @@ ;;; ;;; Z may be any number, but the result is always a complex. (defun complexatan (z)  (declare (number z)) + (declare (type (or rational complex) z)) ;; atan z = i * atanh (i*z) (let* ((iz (complex ( (imagpart z)) (realpart z))) (result (complexatanh iz))) @@ 1020,7 +1039,7 @@ ;;; ;;; Z may be any number, but the result is always a complex. (defun complextan (z)  (declare (number z)) + (declare (type (or rational complex) z)) ;; tan z = i * tanh(i*z) (let* ((iz (complex ( (imagpart z)) (realpart z))) (result (complextanh iz))) Index: sxhash.lisp =================================================================== RCS file: /cvsroot/sbcl/sbcl/src/code/sxhash.lisp,v retrieving revision 1.5 retrieving revision 1.6 diff u d r1.5 r1.6  sxhash.lisp 15 Sep 2003 09:21:38 0000 1.5 +++ sxhash.lisp 17 May 2004 09:44:45 0000 1.6 @@ 17,14 +17,15 @@ ;;; SXHASH of FLOAT values is defined directly in terms of DEFTRANSFORM in ;;; order to avoid boxing. (deftransform sxhash ((x) (singlefloat))  '(let ((bits (singlefloatbits x))) + '(let* ((val (+ 0.0f0 x)) + (bits (singlefloatbits val))) (logxor 66194023 (sxhash (the fixnum (logand mostpositivefixnum (logxor bits (ash bits 7)))))))) (deftransform sxhash ((x) (doublefloat))  '(let* ((val x) + '(let* ((val (+ 0.0d0 x)) (hi (doublefloathighbits val)) (lo (doublefloatlowbits val)) (hilo (logxor hi lo))) 