Changeset 16639 for trunk/source/lib/format.lisp
 Timestamp:
 Nov 17, 2015, 6:21:44 PM (5 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/source/lib/format.lisp
r16338 r16639 1295 1295 # 1296 1296 1297 1298 1299 1300 1297 ;;; Given a nonnegative floating point number, SCALEEXPONENT returns a 1301 1298 ;;; new floating point number Z in the range (0.1, 1.0] and and exponent … … 1304 1301 ;;; JUST do the EXPONENT since thats all we use 1305 1302 1306 1307 (defconstant longlog10of2 0.30103d0) 1308 1309 # 1310 (defun scaleexponent (x) 1311 (if (floatp x ) 1312 (scaleexptaux (abs x) 0.0d0 1.0d0 1.0d1 1.0d1 longlog10of2) 1313 (reportbadarg 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 ;(scaleexptaux x (%splfloat 0) (%splfloat 1) %longfloatten 1318 ; %longfloatonetenth longlog10of2))) 1319 # 1320 1321 ; this dies with floating point overflow (?) if fed leastpositivedoublefloat 1322 1323 (defun scaleexptaux (x zero one ten onetenth log10of2) 1324 (let ((exponent (nthvalue 1 (decodefloat x)))) 1325 (if (= x zero) 1326 (values zero 1) 1327 (let* ((e (round (* exponent log10of2))) 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 onetenth) (values x e))))))))) 1339 # 1340 1341 (defun scaleexponent (n) 1342 (let ((exp (nthvalue 1 (decodefloat n)))) 1343 (values (round (* exp longlog10of2))))) 1344 1303 (defconstant singlefloatmine 1304 (nthvalue 1 (decodefloat leastpositivesinglefloat))) 1305 (defconstant doublefloatmine 1306 (nthvalue 1 (decodefloat leastpositivedoublefloat))) 1307 1308 ;;; Adapted from CMUCL. 1309 1310 ;; This is a modified version of the scale computation from Burger and 1311 ;; Dybvig's paper "Printing floatingpoint 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 floatingpoint 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 accuratescaleexponent (v) 1318 (declare (type float v)) 1319 (if (zerop v) 1320 1 1321 (let ((floatradix 2) ; b 1322 (floatdigits (floatdigits v)) ; p 1323 (mine 1324 (etypecase v 1325 (singlefloat singlefloatmine) 1326 (doublefloat doublefloatmine)))) 1327 (multiplevaluebind (f e) 1328 (integerdecodefloat v) 1329 (let ( ;; FIXME: these even tests assume normal IEEE rounding 1330 ;; mode. I wonder if we should cater for nonnormal? 1331 (highok (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 (singlefloat 1337 (float x 1d0)) 1338 (doublefloat 1339 x)))) 1340 (ceiling ( (the (doublefloat 400d0 400d0) 1341 (log xd 10d0)) 1342 1d10)))) 1343 (fixup (r s m+ k) 1344 (if (if highok 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 (10toe (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 floatradix e)) 1358 (be1 (* be floatradix))) 1359 (if (/= f (expt floatradix (1 floatdigits))) 1360 (setf r (* f be 2) 1361 s 2 1362 m+ be) 1363 (setf r (* f be1 2) 1364 s (* floatradix 2) 1365 m+ be1))) 1366 (if (or (= e mine) 1367 (/= f (expt floatradix (1 floatdigits)))) 1368 (setf r (* f 2) 1369 s (* (expt floatradix ( e)) 2) 1370 m+ 1) 1371 (setf r (* f floatradix 2) 1372 s (* (expt floatradix ( 1 e)) 2) 1373 m+ floatradix))) 1374 (scale r s m+)))))))) 1345 1375 1346 1376 ;;; Page ~ … … 1921 1951 (and (minusp k)(< k ( d)))) 1922 1952 (formaterror "incompatible values for k and d"))) 1923 (when (not exp) (setq exp ( scaleexponent number)))1953 (when (not exp) (setq exp (accuratescaleexponent (abs number)))) 1924 1954 AGAIN 1925 1955 (let* ((expt ( exp k))
Note: See TracChangeset
for help on using the changeset viewer.