Ignore:
Timestamp:
Nov 17, 2015, 6:21:44 PM (5 years ago)
Author:
rme
Message:

Replace scale-exponent with accurate-scaled-exponent, borrowed from
CMUCL.

Closes ticket:536 and ticket:1186.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/source/lib/format.lisp

    r16338 r16639  
    12951295|#
    12961296     
    1297 
    1298 
    1299 
    13001297;;; Given a non-negative floating point number, SCALE-EXPONENT returns a
    13011298;;; new floating point number Z in the range (0.1, 1.0] and and exponent
     
    13041301;;; JUST do the EXPONENT since thats all we use
    13051302
    1306 
    1307 (defconstant long-log10-of-2 0.30103d0)
    1308 
    1309 #|
    1310 (defun scale-exponent (x)
    1311   (if (floatp x )
    1312       (scale-expt-aux (abs x) 0.0d0 1.0d0 1.0d1 1.0d-1 long-log10-of-2)
    1313       (report-bad-arg x 'float)))
    1314 
    1315 #|this is the slisp code that was in the place of the error call above.
    1316   before floatp was put in place of shortfloatp.
    1317       ;(scale-expt-aux x (%sp-l-float 0) (%sp-l-float 1) %long-float-ten
    1318       ;                %long-float-one-tenth long-log10-of-2)))
    1319 |#
    1320 
    1321 ; this dies with floating point overflow (?) if fed least-positive-double-float
    1322 
    1323 (defun scale-expt-aux (x zero one ten one-tenth log10-of-2)
    1324   (let ((exponent (nth-value 1 (decode-float x))))
    1325     (if (= x zero)
    1326       (values zero 1)
    1327       (let* ((e (round (* exponent log10-of-2)))
    1328              (x (if (minusp e)          ;For the end ranges.
    1329                   (* x ten (expt ten (- -1 e)))
    1330                   (/ x ten (expt ten (1- e))))))
    1331         (do ((d ten (* d ten))
    1332              (y x (/ x d))
    1333              (e e (1+ e)))
    1334             ((< y one)
    1335              (do ((m ten (* m ten))
    1336                   (z y (* z m))
    1337                   (e e (1- e)))
    1338                  ((>= z one-tenth) (values x e)))))))))
    1339 |#
    1340 
    1341 (defun scale-exponent (n)
    1342   (let ((exp (nth-value 1 (decode-float n))))
    1343     (values (round (* exp long-log10-of-2)))))
    1344 
     1303(defconstant single-float-min-e
     1304  (nth-value 1 (decode-float least-positive-single-float)))
     1305(defconstant double-float-min-e
     1306  (nth-value 1 (decode-float least-positive-double-float)))
     1307
     1308;;; Adapted from CMUCL.
     1309
     1310;; This is a modified version of the scale computation from Burger and
     1311;; Dybvig's paper "Printing floating-point quickly and accurately."
     1312;; We only want the exponent, so most things not needed for the
     1313;; computation of the exponent have been removed.  We also implemented
     1314;; the floating-point log approximation given in Burger and Dybvig.
     1315;; This is very noticeably faster for large and small numbers.  It is
     1316;; slower for intermediate sized numbers.
     1317(defun accurate-scale-exponent (v)
     1318  (declare (type float v))
     1319  (if (zerop v)
     1320      1
     1321      (let ((float-radix 2)             ; b
     1322            (float-digits (float-digits v)) ; p
     1323            (min-e
     1324             (etypecase v
     1325               (single-float single-float-min-e)
     1326               (double-float double-float-min-e))))
     1327        (multiple-value-bind (f e)
     1328            (integer-decode-float v)
     1329          (let ( ;; FIXME: these even tests assume normal IEEE rounding
     1330                ;; mode.  I wonder if we should cater for non-normal?
     1331                (high-ok (evenp f)))
     1332            ;; We only want the exponent here.
     1333            (labels ((flog (x)
     1334                       (declare (type (float (0.0)) x))
     1335                       (let ((xd (etypecase x
     1336                                   (single-float
     1337                                    (float x 1d0))
     1338                                   (double-float
     1339                                    x))))
     1340                         (ceiling (- (the (double-float -400d0 400d0)
     1341                                          (log xd 10d0))
     1342                                     1d-10))))
     1343                     (fixup (r s m+ k)
     1344                       (if (if high-ok
     1345                               (>= (+ r m+) s)
     1346                               (> (+ r m+) s))
     1347                           (+ k 1)
     1348                           k))
     1349                     (scale (r s m+)
     1350                       (let* ((est (flog v))
     1351                              (scale (the integer (10-to-e (abs est)))))
     1352                         (if (>= est 0)
     1353                             (fixup r (* s scale) m+ est)
     1354                             (fixup (* r scale) s (* m+ scale) est)))))
     1355              (let (r s m+)
     1356                (if (>= e 0)
     1357                    (let* ((be (expt float-radix e))
     1358                           (be1 (* be float-radix)))
     1359                      (if (/= f (expt float-radix (1- float-digits)))
     1360                          (setf r (* f be 2)
     1361                                s 2
     1362                                m+ be)
     1363                          (setf r (* f be1 2)
     1364                                s (* float-radix 2)
     1365                                m+ be1)))
     1366                    (if (or (= e min-e)
     1367                            (/= f (expt float-radix (1- float-digits))))
     1368                        (setf r (* f 2)
     1369                              s (* (expt float-radix (- e)) 2)
     1370                              m+ 1)
     1371                        (setf r (* f float-radix 2)
     1372                              s (* (expt float-radix (- 1 e)) 2)
     1373                              m+ float-radix)))
     1374                (scale r s m+))))))))
    13451375
    13461376;;; Page  ~|
     
    19211951                  (and (minusp k)(< k (- d))))
    19221952          (format-error "incompatible values for k and d")))
    1923       (when (not exp) (setq exp (scale-exponent  number)))
     1953      (when (not exp) (setq exp (accurate-scale-exponent (abs number))))
    19241954      AGAIN
    19251955      (let* ((expt (- exp k))
Note: See TracChangeset for help on using the changeset viewer.