blob: c2df37c2b8702d1323b4364c1f86367efbfca335 (
plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
|
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
|
GitWeb