123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129 |
- diff -r -U1 ccl.orig/lib/format.lisp ccl/lib/format.lisp
- --- ccl.orig/lib/format.lisp 2015-11-07 02:10:10.000000000 +0600
- +++ ccl/lib/format.lisp 2015-11-20 22:51:51.736191995 +0600
- @@ -1296,5 +1296,2 @@
-
- -
- -
- -
- ;;; Given a non-negative floating point number, SCALE-EXPONENT returns a
- @@ -1305,41 +1302,74 @@
-
- -
- -(defconstant long-log10-of-2 0.30103d0)
- -
- -#|
- -(defun scale-exponent (x)
- - (if (floatp x )
- - (scale-expt-aux (abs x) 0.0d0 1.0d0 1.0d1 1.0d-1 long-log10-of-2)
- - (report-bad-arg x 'float)))
- -
- -#|this is the slisp code that was in the place of the error call above.
- - before floatp was put in place of shortfloatp.
- - ;(scale-expt-aux x (%sp-l-float 0) (%sp-l-float 1) %long-float-ten
- - ; %long-float-one-tenth long-log10-of-2)))
- -|#
- -
- -; this dies with floating point overflow (?) if fed least-positive-double-float
- -
- -(defun scale-expt-aux (x zero one ten one-tenth log10-of-2)
- - (let ((exponent (nth-value 1 (decode-float x))))
- - (if (= x zero)
- - (values zero 1)
- - (let* ((e (round (* exponent log10-of-2)))
- - (x (if (minusp e) ;For the end ranges.
- - (* x ten (expt ten (- -1 e)))
- - (/ x ten (expt ten (1- e))))))
- - (do ((d ten (* d ten))
- - (y x (/ x d))
- - (e e (1+ e)))
- - ((< y one)
- - (do ((m ten (* m ten))
- - (z y (* z m))
- - (e e (1- e)))
- - ((>= z one-tenth) (values x e)))))))))
- -|#
- -
- -(defun scale-exponent (n)
- - (let ((exp (nth-value 1 (decode-float n))))
- - (values (round (* exp long-log10-of-2)))))
- -
- +(defconstant single-float-min-e
- + (nth-value 1 (decode-float least-positive-single-float)))
- +(defconstant double-float-min-e
- + (nth-value 1 (decode-float least-positive-double-float)))
- +
- +;;; Adapted from CMUCL.
- +
- +;; This is a modified version of the scale computation from Burger and
- +;; Dybvig's paper "Printing floating-point quickly and accurately."
- +;; We only want the exponent, so most things not needed for the
- +;; computation of the exponent have been removed. We also implemented
- +;; the floating-point log approximation given in Burger and Dybvig.
- +;; This is very noticeably faster for large and small numbers. It is
- +;; slower for intermediate sized numbers.
- +(defun accurate-scale-exponent (v)
- + (declare (type float v))
- + (if (zerop v)
- + 1
- + (let ((float-radix 2) ; b
- + (float-digits (float-digits v)) ; p
- + (min-e
- + (etypecase v
- + (single-float single-float-min-e)
- + (double-float double-float-min-e))))
- + (multiple-value-bind (f e)
- + (integer-decode-float v)
- + (let ( ;; FIXME: these even tests assume normal IEEE rounding
- + ;; mode. I wonder if we should cater for non-normal?
- + (high-ok (evenp f)))
- + ;; We only want the exponent here.
- + (labels ((flog (x)
- + (declare (type (float (0.0)) x))
- + (let ((xd (etypecase x
- + (single-float
- + (float x 1d0))
- + (double-float
- + x))))
- + (ceiling (- (the (double-float -400d0 400d0)
- + (log xd 10d0))
- + 1d-10))))
- + (fixup (r s m+ k)
- + (if (if high-ok
- + (>= (+ r m+) s)
- + (> (+ r m+) s))
- + (+ k 1)
- + k))
- + (scale (r s m+)
- + (let* ((est (flog v))
- + (scale (the integer (10-to-e (abs est)))))
- + (if (>= est 0)
- + (fixup r (* s scale) m+ est)
- + (fixup (* r scale) s (* m+ scale) est)))))
- + (let (r s m+)
- + (if (>= e 0)
- + (let* ((be (expt float-radix e))
- + (be1 (* be float-radix)))
- + (if (/= f (expt float-radix (1- float-digits)))
- + (setf r (* f be 2)
- + s 2
- + m+ be)
- + (setf r (* f be1 2)
- + s (* float-radix 2)
- + m+ be1)))
- + (if (or (= e min-e)
- + (/= f (expt float-radix (1- float-digits))))
- + (setf r (* f 2)
- + s (* (expt float-radix (- e)) 2)
- + m+ 1)
- + (setf r (* f float-radix 2)
- + s (* (expt float-radix (- 1 e)) 2)
- + m+ float-radix)))
- + (scale r s m+))))))))
-
- @@ -1922,3 +1952,3 @@
- (format-error "incompatible values for k and d")))
- - (when (not exp) (setq exp (scale-exponent number)))
- + (when (not exp) (setq exp (accurate-scale-exponent (abs number))))
- AGAIN
|