;;;; IEEE 754 Format
(declare (usual-integrations))
+
+;;; Nomenclature:
+;;;
+;;; base, b
+;;; exponent width, w: number of bits in exponent field
+;;; precision, p: number of bits in significand before and after radix point
+;;; significand: integer in [b^{p - 1}, b^b)
+;;; trailing significand: integer in [0, b^{p - 1}) = [0, b^t)
+;;; trailing significand width: number of digits in trailing significand
\f
-(define (decompose-ieee754-double x)
+(define (decompose-ieee754-binary32 x)
+ (decompose-ieee754-binary x 8 24))
+
+(define (decompose-ieee754-binary64 x)
(decompose-ieee754-binary x 11 53))
+(define (decompose-ieee754-binary128 x)
+ (decompose-ieee754-binary x 15 113))
+
(define (decompose-ieee754-binary x exponent-bits precision)
(receive (base emin emax bias exp-subnormal exp-inf/nan)
(ieee754-binary-parameters exponent-bits precision)
(decompose-ieee754 x base emax precision
(lambda (sign) ;if-zero
(values sign 0 0))
- (lambda (sign scaled-significand) ;if-subnormal
- (assert (= 0 (shift-right scaled-significand precision)))
- (values sign exp-subnormal scaled-significand))
- (lambda (sign exponent scaled-significand) ;if-normal
+ (lambda (sign significand) ;if-subnormal
+ (assert (= 0 (shift-right significand (- precision 1))))
+ (values sign (+ exp-subnormal bias) significand))
+ (lambda (sign exponent significand) ;if-normal
(assert (<= emin exponent emax))
;; The integer part is always 1. Strip it for the binary
;; interchange format.
- (assert (= 1 (shift-right scaled-significand (- precision 1))))
+ (assert (= 1 (shift-right significand (- precision 1))))
(values sign
(+ exponent bias)
- (extract-bit-field (- precision 1) 0 scaled-significand)))
+ (extract-bit-field (- precision 1) 0 significand)))
(lambda (sign) ;if-infinite
(values sign exp-inf/nan 0))
(lambda (sign quiet payload) ;if-nan
(define (ieee754-sign x)
(cond ((< 0 x) 0)
((< x 0) 1)
- ;; Zero -- can't use < directly to detect sign. Elicit a
- ;; computational difference.
- ((negative? (atan x -1)) 1)
- (else 0)))
+ (else
+ ;; Zero -- can't use < directly to detect sign. Elicit a
+ ;; computational difference.
+ (if (negative? (atan x -1))
+ 1
+ 0))))
(define (decompose-ieee754 x base emax precision
if-zero if-subnormal if-normal if-infinite if-nan)
;; bits, so we'll just assume every NaN is a trivial positive
;; signalling NaN and hope the caller has a good story...
(if-nan 0 0 1))
- ((and (< 1. (abs x)) (= x (/ x 2)))
+ ((and (< 1 (abs x)) (= x (/ x 2)))
(if-infinite (if (< 0. x) 0 1)))
(else
- (let ((sign (ieee754-sign x)) (x (abs x)) (emin (- 1 emax)))
+ (let ((sign (ieee754-sign x))
+ (x (abs x))
+ (emin (- 1 emax)))
(define (significand x)
- (truncate->exact (* x (expt base (- precision 1)))))
+ (round->exact (* x (expt base (- precision 1)))))
(cond ((<= 1 x) ;Nonnegative exponent (normal)
(let loop ((exponent 0) (x x))
(cond ((< emax exponent) (if-infinite sign))
((<= base x) (loop (+ exponent 1) (/ x base)))
(else (if-normal sign exponent (significand x))))))
- ((< (expt base emin) x) ;Negative exponent, normal
+ ((<= (expt base emin) x) ;Negative exponent, normal
(let loop ((exponent 0) (x x))
(assert (<= emin exponent))
(if (<= 1 x)
(if-normal sign exponent (significand x))
(loop (- exponent 1) (* x base)))))
((< 0 x) ;Negative exponent, subnormal
- (if (<= x (- (expt base emin) (expt base (- 0 precision))))
+ (if (<= x (/ (expt base (- emin (+ precision 1))) 2))
(if-zero sign)
- (if-subnormal
- sign
- (significand (/ x (expt base emin))))))
+ (if-subnormal sign
+ (significand (/ x (expt base emin))))))
(else
(if-zero sign)))))))
\f
-(define (compose-ieee754-double sign biased-exponent trailing-significand)
+(define (compose-ieee754-binary32 sign biased-exponent trailing-significand)
+ (compose-ieee754-binary sign biased-exponent trailing-significand 8 24))
+
+(define (compose-ieee754-binary64 sign biased-exponent trailing-significand)
(compose-ieee754-binary sign biased-exponent trailing-significand 11 53))
+(define (compose-ieee754-binary128 sign biased-exponent trailing-significand)
+ (compose-ieee754-binary sign biased-exponent trailing-significand 15 113))
+
(define (compose-ieee754-binary sign biased-exponent trailing-significand
exponent-bits precision)
(receive (base emin emax bias exp-subnormal exp-inf/nan)
(ieee754-binary-parameters exponent-bits precision)
- (let ((exponent (- biased-exponent bias)))
+ (let ((exponent (- biased-exponent bias))
+ (t (- precision 1)))
(cond ((= exponent exp-subnormal)
(if (zero? trailing-significand)
- (compose-ieee754-zero sign base emax precision)
+ (compose-ieee754-zero sign)
(compose-ieee754-subnormal sign trailing-significand
- base emax precision)))
+ base emin precision)))
((= exponent exp-inf/nan)
(if (zero? trailing-significand)
- (compose-ieee754-infinity sign base emax precision)
- (let ((p-1 (- precision 1))
- (t trailing-significand))
- (let ((quiet (extract-bit-field 1 p-1 t))
- (payload (extract-bit-field p-1 0 t)))
- (compose-ieee754-nan sign quiet payload
- base emax precision)))))
+ (compose-ieee754-infinity sign)
+ (let ((quiet (extract-bit-field 1 t trailing-significand))
+ (payload (extract-bit-field t 0 trailing-significand)))
+ (compose-ieee754-nan sign quiet payload))))
(else
(assert (<= emin exponent emax))
- (let ((scaled-significand
+ (let ((significand
;; Add the implied integer part of 1.
- (replace-bit-field 1 (- precision 1) trailing-significand
- 1)))
- (compose-ieee754-normal sign exponent scaled-significand
- base emax precision)))))))
+ (replace-bit-field 1 t trailing-significand 1)))
+ (compose-ieee754-normal sign exponent significand
+ base precision)))))))
-(define (compose-ieee754-zero sign base emax precision)
- base emax precision ;ignore
- (* (expt -1 sign) 0))
+(define (compose-ieee754-zero sign)
+ (* (expt -1 sign) 0.))
-(define (compose-ieee754-subnormal sign significand base emax precision)
+(define (compose-ieee754-subnormal sign significand base emin precision)
(* (expt -1 sign)
- (* significand (expt base (- precision emax)))))
+ (* significand (expt base (- emin (- precision 1))))))
-(define (compose-ieee754-normal sign exponent significand base emax precision)
- (assert (<= (- 1 emax) exponent emax))
+(define (compose-ieee754-normal sign exponent significand base precision)
(* (expt -1 sign)
- (expt base exponent)
- (/ significand (expt base (- precision 1)))))
+ (* significand (expt base (- exponent (- precision 1))))))
(define (compose-ieee754-infinity sign)
- (error "Can't compose an IEEE754 infinity!" sign))
+ (* (expt -1 sign)
+ (flo:+inf.0)))
(define (compose-ieee754-nan sign quiet payload)
- (error "Can't compose an IEEE754 NaN!" sign quiet payload))
+ (flo:nan.0))
(define (ieee754-binary-parameters exponent-bits precision)
(assert (zero? (modulo (+ exponent-bits precision) 32)))
(exp-inf/nan (+ emax 1)))
(values base emin emax bias exp-subnormal exp-inf/nan)))))
-(define (ieee754-double-recomposable? x)
+(define (ieee754-binary32-exact? x)
+ (= x
+ (receive (sign biased-exponent trailing-significand)
+ (decompose-ieee754-binary32 x)
+ (compose-ieee754-binary32 sign biased-exponent trailing-significand))))
+
+(define (ieee754-binary64-exact? x)
+ (= x
+ (receive (sign biased-exponent trailing-significand)
+ (decompose-ieee754-binary64 x)
+ (compose-ieee754-binary64 sign biased-exponent trailing-significand))))
+
+(define (ieee754-binary128-exact? x)
(= x
(receive (sign biased-exponent trailing-significand)
- (decompose-ieee754-double x)
- (compose-ieee754-double sign biased-exponent trailing-significand))))
+ (decompose-ieee754-binary64 x)
+ (compose-ieee754-binary128 sign biased-exponent trailing-significand))))
\f
(define (ieee754-binary-hex-string x exponent-bits precision)
(receive (base emin emax bias exp-subnormal exp-inf/nan)
(assert (<= 0 fractional))
(let ((sign (if (zero? sign) "" "-"))
(integer (if (zero? integer) "0" "1"))
+ (dot (if (zero? fractional) "" "."))
(fractional
(if (zero? fractional) "" (number->string fractional #x10)))
- (exponent (number->string exponent #d10)))
- (string-append sign "0x" integer "." fractional "p" exponent)))
+ (expsign (if (< exponent 0) "-" "+"))
+ (exponent (number->string (abs exponent) #d10)))
+ (string-append sign "0x" integer dot fractional "p" expsign exponent)))
(decompose-ieee754 x base emax precision
(lambda (sign) ;if-zero
(numeric sign 0 0 0))
- (lambda (sign scaled-significand) ;if-subnormal
- (assert (< 0 scaled-significand))
- (let ((start (first-set-bit scaled-significand))
- (end (integer-length scaled-significand)))
- (let ((bits (- (- end 1) start)))
+ (lambda (sign significand) ;if-subnormal
+ (assert (< 0 significand))
+ (assert (= 0 (shift-right significand (- precision 1))))
+ (let ((start (first-set-bit significand))
+ (end (integer-length significand)))
+ (let ((fracbits (- (- end 1) start)))
(let ((exponent (- emin (- precision end)))
;; Strip the integer part (1) and the trailing zeros.
(fractional
- (extract-bit-field bits start scaled-significand)))
- ;; Align to a multiple of four bits (hex digits).
- (let* ((alignment (modulo (- 4 (modulo bits 4)) 4))
- (aligned-fractional
- (replace-bit-field bits alignment 0 fractional)))
- (numeric sign 1 aligned-fractional exponent))))))
- (lambda (sign exponent scaled-significand)
- (assert (< 0 scaled-significand))
- (assert (= 1 (shift-right scaled-significand (- precision 1))))
- (let ((useless-zeros (round-down (first-set-bit scaled-significand) 4))
+ (extract-bit-field fracbits start significand)))
+ (numeric sign 1 fractional exponent)))))
+ (lambda (sign exponent significand)
+ (assert (< 0 significand))
+ (assert (= 1 (shift-right significand (- precision 1))))
+ (let ((useless-zeros (round-down (first-set-bit significand) 4))
(fractional
- (extract-bit-field (- precision 1) 0 scaled-significand)))
+ (extract-bit-field (- precision 1) 0 significand)))
(numeric sign 1 (shift-right fractional useless-zeros) exponent)))
(lambda (sign)
(symbolic sign "inf" 0))
(lambda (sign quiet payload)
(symbolic sign (if (zero? quiet) "sNaN" "qNaN") payload)))))
+(define (ieee754-binary32-hex-string x)
+ (ieee754-binary-hex-string x 8 24))
+
+(define (ieee754-binary64-hex-string x)
+ (ieee754-binary-hex-string x 11 53))
+
+(define (ieee754-binary128-hex-string x)
+ (ieee754-binary-hex-string x 15 113))
+
(define (round-up x n)
(* n (quotient (+ x (- n 1)) n)))
--- /dev/null
+#| -*-Scheme-*-
+
+Copyright (C) 1986, 1987, 1988, 1989, 1990, 1991, 1992, 1993, 1994,
+ 1995, 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005,
+ 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016,
+ 2017, 2018 Massachusetts Institute of Technology
+
+This file is part of MIT/GNU Scheme.
+
+MIT/GNU Scheme is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or (at
+your option) any later version.
+
+MIT/GNU Scheme is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with MIT/GNU Scheme; if not, write to the Free Software
+Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301,
+USA.
+
+|#
+
+;;;; Test of IEEE 754 utilities
+\f
+(define (define-enumerated-test prefix elements procedure)
+ (do ((i 0 (+ i 1))
+ (elements elements (cdr elements)))
+ ((not (pair? elements)))
+ (let ((element (car elements)))
+ (define-test (symbol prefix '/ element)
+ (lambda ()
+ (procedure element))))))
+
+(define ((test-ieee754-roundtrip w t bexp-inf/nan compose exact? decompose)
+ bits)
+ (let ((sign (extract-bit-field 1 (+ w t) bits))
+ (biased-exponent (extract-bit-field w t bits))
+ (trailing-significand (extract-bit-field t 0 bits)))
+ (if (not (= biased-exponent bexp-inf/nan))
+ (let ((x (compose sign biased-exponent trailing-significand)))
+ (assert (exact? x))
+ (receive (sign* biased-exponent* trailing-significand*)
+ (decompose x)
+ (assert-= sign* sign)
+ (assert-= biased-exponent* biased-exponent)
+ (assert-= trailing-significand* trailing-significand))))))
+
+(define-test 'binary32-roundtrip-exhaustive
+ (lambda ()
+ (define test
+ (test-ieee754-roundtrip 8 23 255
+ compose-ieee754-binary32
+ ieee754-binary32-exact?
+ decompose-ieee754-binary32))
+ (do ((i 0 (+ i (if keep-it-fast!? 347911 1))))
+ ((>= i (expt 2 32)))
+ (test i))))
+
+(define-enumerated-test 'binary64-roundtrip-selective
+ '(#x0000000000000000
+ #xffffffffffffffff
+ #x0000000000000001
+ #x1000000000000000
+ #x1000000000000001
+ #x7ff0000000000000
+ #xfff0000000000000
+ #x0123456789abcdef
+ #xfedcba9876543210)
+ (test-ieee754-roundtrip 11 52 2047
+ compose-ieee754-binary64
+ ieee754-binary64-exact?
+ decompose-ieee754-binary64))
+
+(define-enumerated-test 'decompose-ieee754-binary64
+ `((0. (zero +))
+ (-0. (zero -))
+ (,(expt 2 -1074) (subnormal + 1))
+ (,(- (expt 2 -1074)) (subnormal - 1))
+ (,(expt 2 -1050) (subnormal + #x1000000))
+ (,(- (expt 2 -1050)) (subnormal - #x1000000))
+ (,(- (expt 2 -1022) (expt 2 -1074)) (subnormal + #xfffffffffffff))
+ (,(- (expt 2 -1074) (expt 2 -1022)) (subnormal - #xfffffffffffff))
+ (,(expt 2 -1022) (normal + -1022 #x10000000000000))
+ (,(- (expt 2 -1022)) (normal - -1022 #x10000000000000))
+ (,(+ (expt 2 -1022) (expt 2 -1074)) (normal + -1022 #x10000000000001))
+ (,(- (+ (expt 2 -1022) (expt 2 -1074))) (normal - -1022 #x10000000000001))
+ (1 (normal + 0 #x10000000000000))
+ (-1 (normal - 0 #x10000000000000))
+ (3/2 (normal + 0 #x18000000000000))
+ (-3/2 (normal - 0 #x18000000000000))
+ (2 (normal + 1 #x10000000000000))
+ (-2 (normal - 1 #x10000000000000))
+ (,(flo:+inf.0) (infinity +))
+ (,(flo:-inf.0) (infinity -))
+ (,(flo:nan.0) (nan + s 1)))
+ (lambda (c)
+ (let ((x (list-ref c 0))
+ (y (list-ref c 1)))
+ (define (signify sign)
+ (case sign
+ ((0) '+)
+ ((1) '-)
+ (else (error "Invalid sign:" sign))))
+ (flo:with-trapped-exceptions 0
+ (lambda ()
+ ((lambda (z)
+ (assert-equal z y)
+ (flo:clear-exceptions! (flo:supported-exceptions)))
+ (let ((exponent-bits 11)
+ (precision 53))
+ (receive (base emin emax bias exp-subnormal exp-inf/nan)
+ (ieee754-binary-parameters exponent-bits precision)
+ emin bias exp-subnormal exp-inf/nan ;ignore
+ (decompose-ieee754 x base emax precision
+ (lambda (sign) `(zero ,(signify sign)))
+ (lambda (sign significand)
+ `(subnormal ,(signify sign) ,significand))
+ (lambda (sign exponent significand)
+ `(normal ,(signify sign) ,exponent ,significand))
+ (lambda (sign)
+ `(infinity ,(signify sign)))
+ (lambda (sign quiet payload)
+ `(nan ,(signify sign)
+ ,(case quiet
+ ((0) 's)
+ ((1) 'q)
+ (else (error "Quiet bit:" quiet)))
+ ,payload)))))))))))
+
+(define-enumerated-test 'ieee754-binary64-hex
+ '((0 "0x0p+0")
+ (-0. "-0x0p+0")
+ (1 "0x1p+0")
+ (-1 "-0x1p+0")
+ (1/2 "0x1p-1")
+ (-1/2 "-0x1p-1")
+ (12345 "0x1.81c8p+13")
+ (123456 "0x1.e24p+16")
+ (1.2061684984132626e-11 "0x1.a862p-37"))
+ (lambda (c)
+ (let ((x (list-ref c 0))
+ (s (list-ref c 1)))
+ (assert-string= (ieee754-binary64-hex-string x) s))))
projects / mit-scheme.git / commitdiff