File: encode-decode-float.lisp

package info (click to toggle)
maxima-sage 5.45.1-2
  • links: PTS, VCS
  • area: main
  • in suites: bookworm, sid, trixie
  • size: 113,788 kB
  • sloc: lisp: 440,833; fortran: 14,665; perl: 14,369; tcl: 10,997; sh: 4,475; makefile: 2,520; ansic: 447; python: 262; xml: 59; awk: 37; sed: 17
file content (136 lines) | stat: -rw-r--r-- 4,553 bytes parent folder | download | duplicates (8)
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
129
130
131
132
133
134
135
136
;; encode-decode-float.lisp
;; Copyright 2007 by Robert Dodier

;; This program is free software; you can redistribute it and/or
;; modify it under the terms of the GNU General Public License.

;; This program has NO WARRANTY, not even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

;; These functions encode integers into 64 bit IEEE 754 floats
;; and decode 64 bit floats into 64 bit integers.
;; These functions cannot handle any other size of float.
;; 
;; Encode float-64 to integer:  SMASH-FLOAT-64-INTO-INTEGER
;; Decode integer to float-64:  CONSTRUCT-FLOAT-64-FROM-INTEGER
;;
;; Write float-64 to output stream:  WRITE-FLOAT-64
;; Read float-64 from input stream:  READ-FLOAT-64
;; Read an unsigned integer (of any size) from input stream: READ-UNSIGNED-INTEGER
;; Write an unsigned integer (of any size) to output stream: WRITE-UNSIGNED-INTEGER
;; Set assumed external byte order for input and output:  DEFINE-EXTERNAL-BYTE-ORDER

(in-package :maxima)

(defun smash-float-into-integer (x)
  (multiple-value-bind
    (significand exponent sign)
    (integer-decode-float x)
    ;; This logic cannot be guaranteed to work -- there is no necessary
    ;; correlation between IEEE 754 and CL floats. Oh well.
    (if (or (typep x 'double-float) (typep x 'long-float))
      (smash-decoded-float-64-into-integer significand exponent sign)
      (smash-decoded-float-32-into-integer significand exponent sign))))

(defun smash-decoded-float-32-into-integer (significand exponent sign)
  (if (= significand 0)
    0
    (dpb
      (if (> sign 0) 0 1)
      (byte 1 (+ 23 8))
      (dpb
        (+ exponent 127 23)
        (byte 8 23)
        (ldb
          (byte 23 0)
          significand)))))

(defun smash-decoded-float-64-into-integer (significand exponent sign)
  (if (= significand 0)
    0
    (dpb
      (if (> sign 0) 0 1)
      (byte 1 (+ 52 11))
      (dpb
        (+ exponent 1023 52)
        (byte 11 52)
        (ldb
          (byte 52 0)
          significand)))))

(defun construct-float-64-from-integer (x)
  (multiple-value-bind
    (significand exponent sign)
    (extract-smashed-float-64-from-integer x)
    (* sign (scale-float (float significand 1d0) exponent))))

(defun extract-smashed-float-64-from-integer (x)
  (if (eql x 0)
    (values 0 0 0)
    (let
      ((significand (dpb x (byte 52 0) #x10000000000000))
       (exponent (- (ldb (byte 11 52) x) 1023 52))
       (sign (if (eql (ldb (byte 1 63) x) 0) 1 -1)))
      (values significand exponent sign))))

;; Stream input and output

(defun write-float (x s)
  (write-unsigned-integer (smash-float-into-integer x) (size-in-bytes x) s))

(defun size-in-bytes (x)
  (if (or (typep x 'double-float) (typep x 'long-float)) 8 4)) ;; AUGHHHH!! THIS IS TERRIBLE!

(defun read-float-64 (s)
  (let ((x (read-unsigned-integer 8 s)))
    (if (eq x 'eof) 'eof (construct-float-64-from-integer x))))

;; READ-UNSIGNED-INTEGER, WRITE-UNSIGNED-INTEGER, and associated
;; byte order stuff adapted from read-bytes-standalone.lisp,
;; by Martin Raspaud and Robert Strandh,
;; which was released under terms of GNU GPL v2 or later.

(deftype external-byte-order ()
  "Defines the legal values for *EXTERNAL-BYTE-ORDER*."
  '(member :msb :lsb))

(defvar *external-byte-order* :msb
  "*EXTERNAL-BYTE-ORDER* must be either :msb or :lsb")

(defun define-external-byte-order (x)
  (check-type x external-byte-order)
  (setf *external-byte-order* x))

(defun read-unsigned-integer (nb-bytes s)
  "Read an unsigned integer of size NB-BYTES bytes from input stream S."
  (if (zerop nb-bytes) 0
    (let (bytes (y 0))
      (dotimes (i nb-bytes)
        (let ((x (read-byte s nil 'eof)))
          (if (eq x 'eof)
            (return-from read-unsigned-integer 'eof)
            (setq bytes (nconc bytes (list x))))))
      (case *external-byte-order*
        (:lsb
          (mapc #'(lambda (x) (setq y (+ x (ash y 8)))) (nreverse bytes)))
        (:msb
          (mapc #'(lambda (x) (setq y (+ x (ash y 8)))) bytes)))
      y)))

(defun write-unsigned-integer (quantity nb-bytes s)
  "Write an unsigned integer of size NB-BYTES bytes to output stream S."
  (case *external-byte-order*
    (:lsb
      (unless (zerop nb-bytes)
        (write-byte (logand quantity #xff) s)
        (write-unsigned-integer
          (ash quantity -8)
          (1- nb-bytes)
          s)))
    (:msb
      (unless (zerop nb-bytes)
        (write-unsigned-integer
          (ash quantity -8)
          (1- nb-bytes)
          s)
        (write-byte (logand quantity #xff) s)))))