1 ;;;; -*- Mode: lisp -*-
3 ;;;; Copyright (c) 2007 Raymond Toy
5 ;;;; Permission is hereby granted, free of charge, to any person
6 ;;;; obtaining a copy of this software and associated documentation
7 ;;;; files (the "Software"), to deal in the Software without
8 ;;;; restriction, including without limitation the rights to use,
9 ;;;; copy, modify, merge, publish, distribute, sublicense, and/or sell
10 ;;;; copies of the Software, and to permit persons to whom the
11 ;;;; Software is furnished to do so, subject to the following
14 ;;;; The above copyright notice and this permission notice shall be
15 ;;;; included in all copies or substantial portions of the Software.
17 ;;;; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
18 ;;;; EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
19 ;;;; OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
20 ;;;; NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
21 ;;;; HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
22 ;;;; WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
23 ;;;; FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
24 ;;;; OTHER DEALINGS IN THE SOFTWARE.
28 ;; Smallest exponent for a double-float.
29 (eval-when (:compile-toplevel
:load-toplevel
:execute
)
30 (defconstant +double-float-min-e
+
37 (defun qd-to-digits (v &optional position relativep
)
38 ;; V is the number to be printed. If RELATIVEP is NIL, POSITION is
39 ;; the number of digits to the left of the decimal point where we
40 ;; want to stop printing. If RELATIVEP is non-NIL, POSITION is the
41 ;; total number of digits we want printed.
43 ;; Two values are returned: k, and the digit string, without a
44 ;; decimal point. k is the index into the string, before which the
45 ;; decimal point would go.
46 (let ((print-base 10) ; B
48 (float-digits (cl:* 4 53)) ; p
49 (min-e +double-float-min-e
+))
50 (multiple-value-bind (f e
)
52 (let ( ;; FIXME: these even tests assume normal IEEE rounding
53 ;; mode. I wonder if we should cater for non-normal?
56 (result (make-array 50 :element-type
'base-char
57 :fill-pointer
0 :adjustable t
)))
58 (labels ((scale (r s m
+ m-
)
59 ;; Keep increasing k until it's big enough
61 (s s
(cl:* s print-base
)))
62 ((not (let ((test (cl:+ r m
+)))
64 (and high-ok
(= test s
)))))
65 ;; k is too big. Decrease until
67 (r r
(cl:* r print-base
))
68 (m+ m
+ (cl:* m
+ print-base
))
69 (m- m-
(cl:* m- print-base
)))
70 ((not (let ((test (cl:* (cl:+ r m
+) print-base
)))
72 (and (not high-ok
) (= test s
)))))
73 ;; k is correct. Generate the digits.
74 (values k
(generate r s m
+ m-
)))))))
76 (multiple-value-bind (d r
)
77 (truncate (cl:* r print-base
) s
)
78 (let ((m+ (cl:* m
+ print-base
))
79 (m- (cl:* m- print-base
)))
80 (let ((tc1 (or (< r m-
) (and low-ok
(= r m-
))))
81 (tc2 (let ((test (cl:+ r m
+)))
83 (and high-ok
(= test s
))))))
85 ((and (not tc1
) (not tc2
))
86 (vector-push-extend (char +digits
+ d
) result
)
87 ;; FIXME sucky tail recursion. This whole
88 ;; kaboodle should be DO*/LOOPified.
90 ;; pedantically keeping all the conditions
91 ;; in so that we can move them around.
93 (vector-push-extend (char +digits
+ (1+ d
)) result
)
96 (vector-push-extend (char +digits
+ d
) result
)
99 (vector-push-extend (char +digits
+
100 (if (< (cl:* r
2) s
) d
(1+ d
)))
105 (let* ((be (expt float-radix e
))
106 (be1 (cl:* be float-radix
)))
107 (if (/= f
(expt float-radix
(1-
109 (setf r
(cl:* f be
2)
113 (setf r
(cl:* f be1
2)
114 s
(cl:* float-radix
2)
118 (/= f
(expt float-radix
(1-
121 s
(cl:* (expt float-radix
(cl:- e
)) 2)
124 (setf r
(cl:* f float-radix
2)
125 s
(cl:* (expt float-radix
(cl:-
1 e
)) 2)
130 ;;(aver (> position 0))
132 ;; running out of letters here
133 (l 1 (cl:* l print-base
)))
134 ((>= (cl:* s l
) (cl:+ r m
+))
136 (if (< (cl:+ r
(cl:* s
(cl:/ (expt print-base
(cl:- k
138 (cl:* s
(expt print-base k
)))
139 (setf position
(cl:- k position
))
140 (setf position
(cl:- k position
1))))))
141 (let ((low (max m-
(cl:/ (cl:* s
(expt print-base
143 (high (max m
+ (cl:/ (cl:* s
(expt print-base
151 (scale r s m
+ m-
)))))))
153 (defun qd-print-exponent (x exp stream
)
154 (let ((*print-radix
* nil
))
155 (format stream
"q~D" exp
)))
157 (defun qd-to-string (x &optional width fdigits scale fmin
)
160 ;;zero is a special case which float-string cannot handle
162 (let ((s (make-string (1+ fdigits
) :initial-element
#\
0)))
163 (setf (schar s
0) #\.
)
164 (values s
(length s
) t
(zerop fdigits
) 0))
165 (values "." 1 t t
0)))
167 (multiple-value-bind (e string
)
169 (qd-to-digits x
(min (- (+ fdigits
(or scale
0)))
171 (if (and width
(> width
1))
172 (let ((w (multiple-value-list
176 (if (and scale
(minusp scale
))
179 (f (multiple-value-list
180 (qd-to-digits x
(- (+ (or fmin
0)
181 (if scale scale
0)))))))
183 ((>= (length (cadr w
)) (length (cadr f
)))
185 (t (values-list f
))))
187 (let ((e (+ e
(or scale
0)))
188 (stream (make-string-output-stream)))
191 (write-string string stream
:end
(min (length string
)
193 (dotimes (i (- e
(length string
)))
194 (write-char #\
0 stream
))
195 (write-char #\. stream
)
196 (write-string string stream
:start
(min (length string
)
199 (dotimes (i (- fdigits
201 (min (length string
) e
))))
202 (write-char #\
0 stream
))))
204 (write-string "." stream
)
206 (write-char #\
0 stream
))
207 ;; If we're out of room (because fdigits is too
208 ;; small), don't print out our string. This fixes
209 ;; things like (format nil "~,2f" 0.001). We should
210 ;; print ".00", not ".001".
211 (when (or (null fdigits
)
212 (plusp (+ e fdigits
)))
213 (write-string string stream
))
215 (dotimes (i (+ fdigits e
(- (length string
))))
216 (write-char #\
0 stream
)))))
217 (let ((string (get-output-stream-string stream
)))
218 (values string
(length string
)
219 (char= (char string
0) #\.
)
220 (char= (char string
(1- (length string
))) #\.
)
221 (position #\. string
))))))))
224 (defun qd-output-aux (x &optional
(stream *standard-output
*))
226 (write-string (if (minusp (float-sign (qd-0 x
)))
230 (multiple-value-bind (e string
)
232 (when (minusp (float-sign (qd-0 x
)))
233 (write-char #\- stream
))
237 (write-string string stream
:end
(min (length string
) e
))
238 (dotimes (i (cl:- e
(length string
)))
239 (write-char #\
0 stream
))
240 (write-char #\. stream
)
241 (write-string string stream
:start
(min (length string
) e
))
242 (when (<= (length string
) e
)
243 (write-char #\
0 stream
))
244 (qd-print-exponent x
0 stream
))
246 (write-string "0." stream
)
247 (dotimes (i (cl:- e
))
248 (write-char #\
0 stream
))
249 (write-string string stream
)
250 (qd-print-exponent x
0 stream
))))
252 ;; Exponential format
253 (write-string string stream
:end
1)
254 (write-char #\. stream
)
255 (write-string string stream
:start
1)
256 ;; CLHS 22.1.3.1.3 says at least one digit must be printed
257 ;; after the decimal point.
258 (when (= (length string
) 1)
259 (write-char #\
0 stream
))
260 (qd-print-exponent x
(1- e
) stream
))))))
262 ;; Function that can be used with FORMAT ~/
264 (defun qd-format (stream arg colon-p at-sign-p
&rest par
)
265 ;; We should do something with colon-p and at-sign-p
266 (declare (ignore colon-p at-sign-p par
))
267 (write-string "#q" stream
)
268 (qd-output-aux arg stream
))
271 (defun qd-output-infinity (x stream
)
273 (write-string "#." stream
))
275 (error 'print-not-readable
:object x
))
277 (write-string "#<" stream
)))
278 (write-string "QD:+QUAD-DOUBLE-FLOAT" stream
)
279 (write-string (if (plusp (qd-0 x
))
283 (write-string "INFINITY+" stream
)
285 (write-string ">" stream
)))
288 (defun qd-output-nan (x stream
)
289 (print-unreadable-object (x stream
)
290 (write-string "QUAD-DOUBLE-FLOAT" stream
)
291 (write-string (if (float-trapping-nan-p (qd-0 x
)) " Trapping" " Quiet") stream
)
292 (write-string " NaN" stream
)))
295 (defun qd-format (stream arg colon-p at-sign-p
&rest par
)
296 (declare (type %quad-double arg
)
298 ;; We should do something with colon-p and at-sign-p
299 (declare (ignore colon-p at-sign-p par
))
300 (cond ((ext:float-infinity-p
(qd-0 arg
))
301 (qd-output-infinity arg stream
))
302 ((ext:float-nan-p
(qd-0 arg
))
303 (qd-output-nan arg stream
))
305 (qd-output-aux arg stream
))))
307 (defun make-float (sign int-part frac-part scale exp
)
308 (declare (type (member -
1 1) sign
)
309 (type unsigned-byte int-part frac-part
)
313 (format t
"sign = ~A~%" sign
)
314 (format t
"int-part = ~A~%" int-part
)
315 (format t
"frac-part = ~A~%" frac-part
)
316 (format t
"scale = ~A~%" scale
)
317 (format t
"exp = ~A~%" exp
))
318 (let ((int (cl:+ (cl:* int-part
(expt 10 scale
))
320 (power (cl:- exp scale
)))
322 (format t
"~A * ~A * 10^(~A)~%" sign int power
)
323 (let* ((len (integer-length int
)))
325 (format t
"len = ~A~%" len
)
327 (let ((xx (make-qd-d (float int
1d0
)))
328 (yy (npow (make-qd-d 10d0
)
332 (format t
"int = ~A~%" int
)
333 (format t
"fl = ~A~%" (float int
1w0
))
334 (format t
"s = ~A~%" sign
)
335 (format t
"sint = ~A~%" (cl:* sign
(float int
1w0
)))
337 (format t
"npow = ~A~%" yy
))
339 (neg-qd (mul-qd xx yy
))
343 ((hi (ldb (byte 106 (cl:- len
106)) int
))
344 (lo (ldb (byte 106 (cl:- len
212)) int
))
345 (xx (make-qd-dd (cl:* sign
(scale-float (float hi
1w0
)
347 (cl:* sign
(scale-float (float lo
1w0
)
348 (cl:- len
106 106)))))
349 (yy (npow (make-qd-d 10d0
)
351 ((q0 (ldb (byte 53 (cl:- len
53)) int
))
352 (q1 (ldb (byte 53 (cl:- len
(* 2 53))) int
))
353 (q2 (ldb (byte 53 (cl:- len
(* 3 53))) int
))
354 (q3 (ldb (byte 53 (cl:- len
(* 4 53))) int
))
355 (xx (make-qd-d (scale-float (float q0
1d0
)
357 (scale-float (float q1
1d0
)
359 (scale-float (float q2
1d0
)
361 (scale-float (float q3
1d0
)
362 (cl:- len
(* 4 53)))))
363 (yy (npow (make-qd-d 10d0
)
367 (format t
"xx = ~A~%" xx
)
369 (format t
" = ~/qd::qd-format/~%" xx
)
370 (format t
"yy = ~A~%" yy
)
372 (format t
" = ~/qd::qd-format/~%" yy
)
373 (format t
"q0 = ~X (~A)~%" q0
374 (scale-float (float q0
1d0
)
376 (format t
"q1 = ~X (~A)~%" q1
377 (scale-float (float q1
1d0
)
378 (cl:- len
(* 2 53))))
380 (format t
"~/qdi::qd-format/~%" (mul-qd xx yy
)))
382 (neg-qd (mul-qd xx yy
))
383 (mul-qd xx yy
))))))))
385 ;; This seems to work, but really needs to be rewritten!
386 (defun read-qd (stream)
387 (labels ((read-digits (s)
388 ;; Read a sequence of digits and return the decimal
389 ;; value, the character that terminated the sequence, and
390 ;; how many characters were read.
395 (setf ch
(peek-char nil s nil
:eof
))
401 (setf val
(cl:+ (digit-char-p ch
)
405 (values ch val count
)))
407 (let ((maybe-sign (peek-char t s nil
:eof
)))
408 (cond ((eql maybe-sign
#\-
)
412 ((eql maybe-sign
#\
+)
415 ((and (characterp maybe-sign
)
416 (digit-char-p maybe-sign
))
418 ((eql maybe-sign
#\.
)
423 (let ((exp-sign (read-sign s
)))
424 (when (eq :eof exp-sign
)
425 (return-from read-exp
0))
426 (when (eq :eof
(peek-char t s nil
:eof
))
427 (error "Sign of exponent, but no value"))
428 (multiple-value-bind (char expo
)
430 (declare (ignore char
))
431 (cl:* exp-sign expo
)))))
432 (let ((sign (read-sign stream
))
437 (when (eq :eof
(peek-char t stream nil
:eof
))
438 (error "Sign, but no value"))
439 (multiple-value-bind (char int
)
442 (cond ((eql :eof char
)
446 (multiple-value-bind (char val scale-val
)
448 (declare (ignore char
))
450 (setf scale scale-val
)
451 (let ((next (peek-char nil stream nil
:eof
)))
452 (when (equalp next
#\q
)
454 (let ((exp-sign (read-sign stream
)))
455 (setf exp
(read-exp stream
))
456 (setf exp
(cl:* exp exp-sign
)))))))
459 (setf exp
(read-exp stream
))
461 (make-float sign int-part frac-part scale exp
)))))
463 (defun qd-reader (stream subchar arg
)
466 (defun qd-from-string (string)
467 (cl::with-input-from-string
(s string
)
471 (set-dispatch-macro-character #\
# #\Q
#'qd-reader
)