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 smash-decoded-float-64-into-integer (significand exponent sign)
  49   (if (= significand 0)
  50     0
  51     (dpb
  52       (if (> sign 0) 0 1)
  53       (byte 1 (+ 52 11))
  54       (dpb
  55         (+ exponent 1023 52)
  56         (byte 11 52)
  57         (ldb
  58           (byte 52 0)
  59           significand)))))
  60
  61 (defun construct-float-64-from-integer (x)
  62   (multiple-value-bind
  63     (significand exponent sign)
  64     (extract-smashed-float-64-from-integer x)
  65     (* sign (scale-float (float significand 1d0) exponent))))
  66
  67 (defun extract-smashed-float-64-from-integer (x)
  68   (if (eql x 0)
  69     (values 0 0 0)
  70     (let
  71       ((significand (dpb x (byte 52 0) #x10000000000000))
  72        (exponent (- (ldb (byte 11 52) x) 1023 52))
  73        (sign (if (eql (ldb (byte 1 63) x) 0) 1 -1)))
  74       (values significand exponent sign))))
  75
  76 ;; Stream input and output
  77
  78 (defun write-float (x s)
  79   (write-unsigned-integer (smash-float-into-integer x) (size-in-bytes x) s))
  80
  81 (defun size-in-bytes (x)
  82   (if (or (typep x 'double-float) (typep x 'long-float)) 8 4)) ;; AUGHHHH!! THIS IS TERRIBLE!
  83
  84 (defun read-float-64 (s)
  85   (let ((x (read-unsigned-integer 8 s)))
  86     (if (eq x 'eof) 'eof (construct-float-64-from-integer x))))
  87
  88 ;; READ-UNSIGNED-INTEGER, WRITE-UNSIGNED-INTEGER, and associated
  89 ;; byte order stuff adapted from read-bytes-standalone.lisp,
  90 ;; by Martin Raspaud and Robert Strandh,
  91 ;; which was released under terms of GNU GPL v2 or later.
  92
  93 (deftype external-byte-order ()
  94   "Defines the legal values for *EXTERNAL-BYTE-ORDER*."
  95   '(member :msb :lsb))
  96
  97 (defvar *external-byte-order* :msb
  98   "*EXTERNAL-BYTE-ORDER* must be either :msb or :lsb")
  99
 100 (defun define-external-byte-order (x)
 101   (check-type x external-byte-order)
 102   (setf *external-byte-order* x))
 103
 104 (defun read-unsigned-integer (nb-bytes s)
 105   "Read an unsigned integer of size NB-BYTES bytes from input stream S."
 106   (if (zerop nb-bytes) 0
 107     (let (bytes (y 0))
 108       (dotimes (i nb-bytes)
 109         (let ((x (read-byte s nil 'eof)))
 110           (if (eq x 'eof)
 111             (return-from read-unsigned-integer 'eof)
 112             (setq bytes (nconc bytes (list x))))))
 113       (case *external-byte-order*
 114         (:lsb
 115           (mapc #'(lambda (x) (setq y (+ x (ash y 8)))) (nreverse bytes)))
 116         (:msb
 117           (mapc #'(lambda (x) (setq y (+ x (ash y 8)))) bytes)))
 118       y)))
 119
 120 (defun write-unsigned-integer (quantity nb-bytes s)
 121   "Write an unsigned integer of size NB-BYTES bytes to output stream S."
 122   (case *external-byte-order*
 123     (:lsb
 124       (unless (zerop nb-bytes)
 125         (write-byte (logand quantity #xff) s)
 126         (write-unsigned-integer
 127           (ash quantity -8)
 128           (1- nb-bytes)
 129           s)))
 130     (:msb
 131       (unless (zerop nb-bytes)
 132         (write-unsigned-integer
 133           (ash quantity -8)
 134           (1- nb-bytes)
 135           s)
 136         (write-byte (logand quantity #xff) s)))))