share/numericalio/encode-decode-float.lisp

   1 ;; encode-decode-float.lisp
   2 ;; Copyright 2007 by Robert Dodier
   3
   4 ;; This program is free software; you can redistribute it and/or
   5 ;; modify it under the terms of the GNU General Public License.
   6
   7 ;; This program has NO WARRANTY, not even the implied warranty of
   8 ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
   9
  10 ;; These functions encode integers into 64 bit IEEE 754 floats
  11 ;; and decode 64 bit floats into 64 bit integers.
  12 ;; These functions cannot handle any other size of float.
  13 ;;
  14 ;; Encode float-64 to integer:  SMASH-FLOAT-64-INTO-INTEGER
  15 ;; Decode integer to float-64:  CONSTRUCT-FLOAT-64-FROM-INTEGER
  16 ;;
  17 ;; Write float-64 to output stream:  WRITE-FLOAT-64
  18 ;; Read float-64 from input stream:  READ-FLOAT-64
  19 ;; Read an unsigned integer (of any size) from input stream: READ-UNSIGNED-INTEGER
  20 ;; Write an unsigned integer (of any size) to output stream: WRITE-UNSIGNED-INTEGER
  21 ;; Set assumed external byte order for input and output:  DEFINE-EXTERNAL-BYTE-ORDER
  22
  23 (in-package :maxima)
  24
  25 (defun smash-float-into-integer (x)
  26   (multiple-value-bind
  27     (significand exponent sign)
  28     (integer-decode-float x)
  29     ;; This logic cannot be guaranteed to work -- there is no necessary
  30     ;; correlation between IEEE 754 and CL floats. Oh well.
  31     (if (or (typep x 'double-float) (typep x 'long-float))
  32       (smash-decoded-float-64-into-integer significand exponent sign)
  33       (smash-decoded-float-32-into-integer significand exponent sign))))
  34
  35 (defun smash-decoded-float-32-into-integer (significand exponent sign)
  36   (if (= significand 0)
  37     0
  38     (dpb
  39       (if (> sign 0) 0 1)
  40       (byte 1 (+ 23 8))
  41       (dpb
  42         (+ exponent 127 23)
  43         (byte 8 23)
  44         (ldb
  45           (byte 23 0)
  46           significand)))))
  47
  48 (defun construct-float-32-from-integer (x)
  49   (multiple-value-bind
  50     (significand exponent sign)
  51     (extract-smashed-float-32-from-integer x)
  52     (* sign (scale-float (float significand 1f0) exponent))))
  53
  54 (defun extract-smashed-float-32-from-integer (x)
  55   (if (eql x 0)
  56     (values 0 0 0)
  57     (let
  58       ((significand (dpb x (byte 23 0) #x800000))
  59        (exponent (- (ldb (byte 8 23) x) 127 23))
  60        (sign (if (eql (ldb (byte 1 31) x) 0) 1 -1)))
  61       (values significand exponent sign))))
  62
  63 (defun smash-decoded-float-64-into-integer (significand exponent sign)
  64   (if (= significand 0)
  65     0
  66     (dpb
  67       (if (> sign 0) 0 1)
  68       (byte 1 (+ 52 11))
  69       (dpb
  70         (+ exponent 1023 52)
  71         (byte 11 52)
  72         (ldb
  73           (byte 52 0)
  74           significand)))))
  75
  76 (defun construct-float-64-from-integer (x)
  77   (multiple-value-bind
  78     (significand exponent sign)
  79     (extract-smashed-float-64-from-integer x)
  80     (* sign (scale-float (float significand 1d0) exponent))))
  81
  82 (defun extract-smashed-float-64-from-integer (x)
  83   (if (eql x 0)
  84     (values 0 0 0)
  85     (let
  86       ((significand (dpb x (byte 52 0) #x10000000000000))
  87        (exponent (- (ldb (byte 11 52) x) 1023 52))
  88        (sign (if (eql (ldb (byte 1 63) x) 0) 1 -1)))
  89       (values significand exponent sign))))
  90
  91 ;; Stream input and output
  92
  93 (defun write-float (x s)
  94   (write-unsigned-integer (smash-float-into-integer x) (size-in-bytes x) s))
  95
  96 (defun size-in-bytes (x)
  97   (if (or (typep x 'double-float) (typep x 'long-float)) 8 4)) ;; AUGHHHH!! THIS IS TERRIBLE!
  98
  99 (defun read-float-64 (s)
 100   (let ((x (read-unsigned-integer 8 s)))
 101     (if (eq x 'eof) 'eof (construct-float-64-from-integer x))))
 102
 103 ;; READ-UNSIGNED-INTEGER, WRITE-UNSIGNED-INTEGER, and associated
 104 ;; byte order stuff adapted from read-bytes-standalone.lisp,
 105 ;; by Martin Raspaud and Robert Strandh,
 106 ;; which was released under terms of GNU GPL v2 or later.
 107
 108 (deftype external-byte-order ()
 109   "Defines the legal values for *EXTERNAL-BYTE-ORDER*."
 110   '(member :msb :lsb))
 111
 112 (defvar *external-byte-order* :msb
 113   "*EXTERNAL-BYTE-ORDER* must be either :msb or :lsb")
 114
 115 (defun define-external-byte-order (x)
 116   (check-type x external-byte-order)
 117   (setf *external-byte-order* x))
 118
 119 (defun read-unsigned-integer (nb-bytes s)
 120   "Read an unsigned integer of size NB-BYTES bytes from input stream S."
 121   (if (zerop nb-bytes) 0
 122     (let (bytes (y 0))
 123       (dotimes (i nb-bytes)
 124         (let ((x (read-byte s nil 'eof)))
 125           (if (eq x 'eof)
 126             (return-from read-unsigned-integer 'eof)
 127             (setq bytes (nconc bytes (list x))))))
 128       (case *external-byte-order*
 129         (:lsb
 130           (mapc #'(lambda (x) (setq y (+ x (ash y 8)))) (nreverse bytes)))
 131         (:msb
 132           (mapc #'(lambda (x) (setq y (+ x (ash y 8)))) bytes)))
 133       y)))
 134
 135 (defun write-unsigned-integer (quantity nb-bytes s)
 136   "Write an unsigned integer of size NB-BYTES bytes to output stream S."
 137   (case *external-byte-order*
 138     (:lsb
 139       (unless (zerop nb-bytes)
 140         (write-byte (logand quantity #xff) s)
 141         (write-unsigned-integer
 142           (ash quantity -8)
 143           (1- nb-bytes)
 144           s)))
 145     (:msb
 146       (unless (zerop nb-bytes)
 147         (write-unsigned-integer
 148           (ash quantity -8)
 149           (1- nb-bytes)
 150           s)
 151         (write-byte (logand quantity #xff) s)))))