external/oct/qd-io.lisp

   1 ;;;; -*- Mode: lisp -*-
   2 ;;;;
   3 ;;;; Copyright (c) 2007 Raymond Toy
   4 ;;;;
   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
  12 ;;;; conditions:
  13 ;;;;
  14 ;;;; The above copyright notice and this permission notice shall be
  15 ;;;; included in all copies or substantial portions of the Software.
  16 ;;;;
  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.
  25
  26 (in-package #:qdi)
  27
  28 ;; Smallest exponent for a double-float.
  29 (eval-when (:compile-toplevel :load-toplevel :execute)
  30 (defconstant +double-float-min-e+
  31   -1073)
  32
  33 (defconstant +digits+
  34   "0123456789")
  35 ) ; eval-when
  36
  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.
  42   ;;
  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
  47         (float-radix 2)                 ; b
  48         (float-digits (cl:* 4 53)) ; p
  49         (min-e +double-float-min-e+))
  50     (multiple-value-bind (f e)
  51         (integer-decode-qd v)
  52       (let ( ;; FIXME: these even tests assume normal IEEE rounding
  53             ;; mode.  I wonder if we should cater for non-normal?
  54             (high-ok (evenp f))
  55             (low-ok (evenp f))
  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
  60                    (do ((k 0 (1+ k))
  61                         (s s (cl:* s print-base)))
  62                        ((not (let ((test (cl:+ r m+)))
  63                                (or (> test s)
  64                                    (and high-ok (= test s)))))
  65                         ;; k is too big.  Decrease until
  66                         (do ((k k (1- k))
  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)))
  71                                     (or (< test s)
  72                                         (and (not high-ok) (= test s)))))
  73                              ;; k is correct.  Generate the digits.
  74                              (values k (generate r s m+ m-)))))))
  75                  (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+)))
  82                                     (or (> test s)
  83                                         (and high-ok (= test s))))))
  84                          (cond
  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.
  89                             (generate r s m+ m-))
  90                            ;; pedantically keeping all the conditions
  91                            ;; in so that we can move them around.
  92                            ((and (not tc1) tc2)
  93                             (vector-push-extend (char +digits+ (1+ d)) result)
  94                             result)
  95                            ((and tc1 (not tc2))
  96                             (vector-push-extend (char +digits+ d) result)
  97                             result)
  98                            ((and tc1 tc2)
  99                             (vector-push-extend (char +digits+
 100                                                       (if (< (cl:* r 2) s) d (1+ d)))
 101                                                 result)
 102                             result)))))))
 103           (let (r s m+ m-)
 104             (if (>= e 0)
 105                 (let* ((be (expt float-radix e))
 106                        (be1 (cl:* be float-radix)))
 107                   (if (/= f (expt float-radix (1-
 108                                                float-digits)))
 109                       (setf r (cl:* f be 2)
 110                             s 2
 111                             m+ be
 112                             m- be)
 113                       (setf r (cl:* f be1 2)
 114                             s (cl:* float-radix 2)
 115                             m+ be1
 116                             m- be)))
 117                 (if (or (= e min-e)
 118                         (/= f (expt float-radix (1-
 119                                                  float-digits))))
 120                     (setf r (cl:* f 2)
 121                           s (cl:* (expt float-radix (cl:- e)) 2)
 122                           m+ 1
 123                           m- 1)
 124                     (setf r (cl:* f float-radix 2)
 125                           s (cl:* (expt float-radix (cl:- 1 e)) 2)
 126                           m+ float-radix
 127                           m- 1)))
 128             (when position
 129               (when relativep
 130                 ;;(aver (> position 0))
 131                 (do ((k 0 (1+ k))
 132                      ;; running out of letters here
 133                      (l 1 (cl:* l print-base)))
 134                     ((>= (cl:* s l) (cl:+ r m+))
 135                      ;; k is now \hat{k}
 136                      (if (< (cl:+ r (cl:* s (cl:/ (expt print-base (cl:- k
 137                                                              position)) 2)))
 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
 142                                                position)) 2)))
 143                     (high (max m+ (cl:/ (cl:* s (expt print-base
 144                                                 position)) 2))))
 145                 (when (<= m- low)
 146                   (setf m- low)
 147                   (setf low-ok t))
 148                 (when (<= m+ high)
 149                   (setf m+ high)
 150                   (setf high-ok t))))
 151             (scale r s m+ m-)))))))
 152
 153 (defun qd-print-exponent (x exp stream)
 154   (let ((*print-radix* nil))
 155     (format stream "q~D" exp)))
 156
 157 (defun qd-to-string (x &optional width fdigits scale fmin)
 158   (setf x (abs-qd x))
 159   (cond ((zerop-qd x)
 160          ;;zero is a special case which float-string cannot handle
 161          (if fdigits
 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)))
 166         (t
 167          (multiple-value-bind (e string)
 168              (if fdigits
 169                  (qd-to-digits x (min (- (+ fdigits (or scale 0)))
 170                                           (- (or fmin 0))))
 171                  (if (and width (> width 1))
 172                      (let ((w (multiple-value-list
 173                                (qd-to-digits x
 174                                              (max 0
 175                                                   (+ (1- width)
 176                                                      (if (and scale (minusp scale))
 177                                                          scale 0)))
 178                                              t)))
 179                            (f (multiple-value-list
 180                                (qd-to-digits x (- (+ (or fmin 0)
 181                                                      (if scale scale 0)))))))
 182                        (cond
 183                          ((>= (length (cadr w)) (length (cadr f)))
 184                           (values-list w))
 185                          (t (values-list f))))
 186                      (qd-to-digits x)))
 187            (let ((e (+ e (or scale 0)))
 188                  (stream (make-string-output-stream)))
 189              (if (plusp e)
 190                  (progn
 191                    (write-string string stream :end (min (length string)
 192                                                          e))
 193                    (dotimes (i (- e (length string)))
 194                      (write-char #\0 stream))
 195                    (write-char #\. stream)
 196                    (write-string string stream :start (min (length string)
 197                                                            e))
 198                    (when fdigits
 199                      (dotimes (i (- fdigits
 200                                     (- (length string)
 201                                        (min (length string) e))))
 202                        (write-char #\0 stream))))
 203                  (progn
 204                    (write-string "." stream)
 205                    (dotimes (i (- e))
 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))
 214                    (when fdigits
 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))))))))
 222
 223
 224 (defun qd-output-aux (x &optional (stream *standard-output*))
 225   (if (zerop-qd x)
 226       (write-string (if (minusp (float-sign (qd-0 x)))
 227                         "-0.0q0"
 228                         "0.0q0")
 229                     stream)
 230       (multiple-value-bind (e string)
 231           (qd-to-digits x)
 232         (when (minusp (float-sign (qd-0 x)))
 233           (write-char #\- stream))
 234         (cond ((< -3 e 8)
 235                ;; Free format
 236                (cond ((plusp e)
 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))
 245                      (t
 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))))
 251               (t
 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))))))
 261
 262 ;; Function that can be used with FORMAT ~/
 263 #-cmu
 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))
 269
 270 #+cmu
 271 (defun qd-output-infinity (x stream)
 272   (cond (*read-eval*
 273          (write-string "#." stream))
 274         (*print-readably*
 275          (error 'print-not-readable :object x))
 276         (t
 277          (write-string "#<" stream)))
 278   (write-string "QD:+QUAD-DOUBLE-FLOAT" stream)
 279   (write-string (if (plusp (qd-0 x))
 280                     "-POSITIVE-"
 281                     "-NEGATIVE-")
 282                 stream)
 283   (write-string "INFINITY+" stream)
 284   (unless *read-eval*
 285     (write-string ">" stream)))
 286
 287 #+cmu
 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)))
 293
 294 #+cmu
 295 (defun qd-format (stream arg colon-p at-sign-p &rest par)
 296   (declare (type %quad-double arg)
 297            (stream stream))
 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))
 304         (t
 305          (qd-output-aux arg stream))))
 306
 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)
 310            (fixnum scale exp))
 311   #+(or)
 312   (progn
 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))
 319                    frac-part))
 320         (power (cl:- exp scale)))
 321     #+(or)
 322     (format t "~A * ~A * 10^(~A)~%" sign int power)
 323     (let* ((len (integer-length int)))
 324       #+(or)
 325       (format t "len = ~A~%" len)
 326       (cond ((<= len 53)
 327              (let ((xx (make-qd-d (float int 1d0)))
 328                    (yy (npow (make-qd-d 10d0)
 329                              power)))
 330                #+(or)
 331                (progn
 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)))
 336                  (format t "~A~%" xx)
 337                  (format t "npow = ~A~%" yy))
 338                (if (minusp sign)
 339                    (neg-qd (mul-qd xx yy))
 340                    (mul-qd xx yy))))
 341             (t
 342              (let* #+nil
 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)
 346                                                         (cl:- len 106)))
 347                                 (cl:* sign (scale-float (float lo 1w0)
 348                                                         (cl:- len 106 106)))))
 349                 (yy (npow (make-qd-d 10d0)
 350                           power)))
 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)
 356                                             (cl:- len 53))
 357                                (scale-float (float q1 1d0)
 358                                             (cl:- len (* 2 53)))
 359                                (scale-float (float q2 1d0)
 360                                             (cl:- len (* 3 53)))
 361                                (scale-float (float q3 1d0)
 362                                             (cl:- len (* 4 53)))))
 363                 (yy (npow (make-qd-d 10d0)
 364                           power)))
 365                #+(or)
 366                (progn
 367                  (format t "xx = ~A~%" xx)
 368                  #+(or)
 369                  (format t "   = ~/qd::qd-format/~%" xx)
 370                  (format t "yy = ~A~%" yy)
 371                  #+(or)
 372                  (format t "   = ~/qd::qd-format/~%" yy)
 373                  (format t "q0 = ~X (~A)~%" q0
 374                          (scale-float (float q0 1d0)
 375                                       (cl:- len 53)))
 376                  (format t "q1 = ~X (~A)~%" q1
 377                          (scale-float (float q1 1d0)
 378                                       (cl:- len (* 2 53))))
 379                  #+(or)
 380                  (format t "~/qdi::qd-format/~%" (mul-qd xx yy)))
 381                (if (minusp sign)
 382                    (neg-qd (mul-qd xx yy))
 383                    (mul-qd xx yy))))))))
 384
 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.
 391              (let ((val 0)
 392                    (ch nil)
 393                    (count 0))
 394                (loop
 395                    (setf ch (peek-char nil s nil :eof))
 396                    (cond ((eq ch :eof)
 397                           (return))
 398                          ((digit-char-p ch)
 399                           (read-char s)
 400                           (incf count)
 401                           (setf val (cl:+ (digit-char-p ch)
 402                                        (cl:* 10 val))))
 403                          (t
 404                           (return))))
 405                (values ch val count)))
 406            (read-sign (s)
 407              (let ((maybe-sign (peek-char t s nil :eof)))
 408                (cond ((eql maybe-sign #\-)
 409                       (read-char s)
 410                       -1
 411                       )
 412                      ((eql maybe-sign #\+)
 413                       (read-char s)
 414                       +1)
 415                      ((and (characterp maybe-sign)
 416                            (digit-char-p maybe-sign))
 417                       +1)
 418                      ((eql maybe-sign #\.)
 419                       +1)
 420                      (t
 421                       maybe-sign))))
 422            (read-exp (s)
 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)
 429                    (read-digits s)
 430                  (declare (ignore char))
 431                  (cl:* exp-sign expo)))))
 432     (let ((sign (read-sign stream))
 433           (int-part 0)
 434           (frac-part 0)
 435           (exp 0)
 436           (scale 0))
 437       (when (eq :eof (peek-char t stream nil :eof))
 438         (error "Sign, but no value"))
 439       (multiple-value-bind (char int)
 440           (read-digits stream)
 441         (setf int-part int)
 442         (cond ((eql :eof char)
 443                )
 444               ((eql char #\.)
 445                (read-char stream)
 446                (multiple-value-bind (char val scale-val)
 447                    (read-digits stream)
 448                  (declare (ignore char))
 449                  (setf frac-part val)
 450                  (setf scale scale-val)
 451                  (let ((next (peek-char nil stream nil :eof)))
 452                    (when (equalp next #\q)
 453                      (read-char stream)
 454                      (let ((exp-sign (read-sign stream)))
 455                        (setf exp (read-exp stream))
 456                        (setf exp (cl:* exp exp-sign)))))))
 457               ((equalp char #\q)
 458                (read-char stream)
 459                (setf exp (read-exp stream))
 460                ))
 461         (make-float sign int-part frac-part scale exp)))))
 462
 463 (defun qd-reader (stream subchar arg)
 464   (read-qd stream))
 465
 466 (defun qd-from-string (string)
 467   (cl::with-input-from-string (s string)
 468     (read-qd s)))
 469
 470 #+nil
 471 (set-dispatch-macro-character #\# #\Q #'qd-reader)
 472