cl-pack/ieee-floats/ieee-floats.lisp

139 lines
6.2 KiB
Common Lisp

;;; Functions for converting floating point numbers represented in
;;; IEEE 754 style to lisp numbers.
;;;
;;; See http://common-lisp.net/project/ieee-floats/
(in-package :common-lisp)
(defpackage :ieee-floats
(:use :common-lisp)
(:export :make-float-converters
:encode-float32
:decode-float32
:encode-float64
:decode-float64))
(in-package :ieee-floats)
;; The following macro may look a bit overcomplicated to the casual
;; reader. The main culprit is the fact that NaN and infinity can be
;; optionally included, which adds a bunch of conditional parts.
;;
;; Assuming you already know more or less how floating point numbers
;; are typically represented, I'll try to elaborate a bit on the more
;; confusing parts, as marked by letters:
;;
;; (A) Exponents in IEEE floats are offset by half their range, for
;; example with 8 exponent bits a number with exponent 2 has 129
;; stored in its exponent field.
;;
;; (B) The maximum possible exponent is reserved for special cases
;; (NaN, infinity).
;;
;; (C) If the exponent fits in the exponent-bits, we have to adjust
;; the significand for the hidden bit. Because decode-float will
;; return a significand between 0 and 1, and we want one between 1
;; and 2 to be able to hide the hidden bit, we double it and then
;; subtract one (the hidden bit) before converting it to integer
;; representation (to adjust for this, 1 is subtracted from the
;; exponent earlier). When the exponent is too small, we set it to
;; zero (meaning no hidden bit, exponent of 1), and adjust the
;; significand downward to compensate for this.
;;
;; (D) Here the hidden bit is added. When the exponent is 0, there is
;; no hidden bit, and the exponent is interpreted as 1.
;;
;; (E) Here the exponent offset is subtracted, but also an extra
;; factor to account for the fact that the bits stored in the
;; significand are supposed to come after the 'decimal dot'.
(defmacro make-float-converters (encoder-name
decoder-name
exponent-bits
significand-bits
support-nan-and-infinity-p)
"Writes an encoder and decoder function for floating point
numbers with the given amount of exponent and significand
bits (plus an extra sign bit). If support-nan-and-infinity-p is
true, the decoders will also understand these special cases. NaN
is represented as :not-a-number, and the infinities as
:positive-infinity and :negative-infinity. Note that this means
that the in- or output of these functions is not just floating
point numbers anymore, but also keywords."
(let* ((total-bits (+ 1 exponent-bits significand-bits))
(exponent-offset (1- (expt 2 (1- exponent-bits)))) ; (A)
(sign-part `(ldb (byte 1 ,(1- total-bits)) bits))
(exponent-part `(ldb (byte ,exponent-bits ,significand-bits) bits))
(significand-part `(ldb (byte ,significand-bits 0) bits))
(nan support-nan-and-infinity-p)
(max-exponent (1- (expt 2 exponent-bits)))) ; (B)
`(progn
(defun ,encoder-name (float)
,@(unless nan `((declare (type float float))))
(multiple-value-bind (sign significand exponent)
(cond ,@(when nan `(((eq float :not-a-number)
(values 0 1 ,max-exponent))
((eq float :positive-infinity)
(values 0 0 ,max-exponent))
((eq float :negative-infinity)
(values 1 0 ,max-exponent))))
((zerop float)
(values 0 0 0))
(t
(multiple-value-bind (significand exponent sign) (decode-float float)
(let ((exponent (+ (1- exponent) ,exponent-offset))
(sign (if (= sign 1.0) 0 1)))
(unless (< exponent ,(expt 2 exponent-bits))
(error "Floating point overflow when encoding ~A." float))
(if (< exponent 0) ; (C)
(values sign (ash (round (* ,(expt 2 significand-bits) significand)) exponent) 0)
(values sign (round (* ,(expt 2 significand-bits) (1- (* significand 2)))) exponent))))))
(let ((bits 0))
(declare (type (unsigned-byte ,total-bits) bits))
(setf ,sign-part sign
,exponent-part exponent
,significand-part significand)
bits)))
(defun ,decoder-name (bits)
(declare (type (unsigned-byte ,total-bits) bits))
(let* ((sign ,sign-part)
(exponent ,exponent-part)
(significand ,significand-part))
,@(when nan `((when (= exponent ,max-exponent)
(return-from ,decoder-name
(cond ((not (zerop significand)) :not-a-number)
((zerop sign) :positive-infinity)
(t :negative-infinity))))))
(if (zerop exponent) ; (D)
(setf exponent 1)
(setf (ldb (byte 1 ,significand-bits) significand) 1))
(unless (zerop sign)
(setf significand (- significand)))
(scale-float (float significand ,(if (> total-bits 32) 1.0d0 1.0))
(- exponent ,(+ exponent-offset significand-bits)))))))) ; (E)
;; And instances of the above for the common forms of floats.
(make-float-converters encode-float32 decode-float32 8 23 nil)
(make-float-converters encode-float64 decode-float64 11 52 nil)
;;; Copyright (c) 2006 Marijn Haverbeke
;;;
;;; This software is provided 'as-is', without any express or implied
;;; warranty. In no event will the authors be held liable for any
;;; damages arising from the use of this software.
;;;
;;; Permission is granted to anyone to use this software for any
;;; purpose, including commercial applications, and to alter it and
;;; redistribute it freely, subject to the following restrictions:
;;;
;;; 1. The origin of this software must not be misrepresented; you must
;;; not claim that you wrote the original software. If you use this
;;; software in a product, an acknowledgment in the product
;;; documentation would be appreciated but is not required.
;;;
;;; 2. Altered source versions must be plainly marked as such, and must
;;; not be misrepresented as being the original software.
;;;
;;; 3. This notice may not be removed or altered from any source
;;; distribution.