From: Taylor R Campbell Date: Mon, 5 Nov 2018 04:40:00 +0000 (+0000) Subject: Tidy up ieee754.scm, add some tests, and export some of it. X-Git-Tag: mit-scheme-pucked-10.1.2~16^2~121 X-Git-Url: https://birchwood-abbey.net/git?a=commitdiff_plain;h=6e89a22bebd9842ff4c3b95eac87b445b1d05677;p=mit-scheme.git Tidy up ieee754.scm, add some tests, and export some of it. --- diff --git a/src/runtime/ieee754.scm b/src/runtime/ieee754.scm index debcb48d2..0be451060 100644 --- a/src/runtime/ieee754.scm +++ b/src/runtime/ieee754.scm @@ -24,27 +24,42 @@ USA. ;;;; 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 -(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 @@ -57,10 +72,12 @@ USA. (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) @@ -71,82 +88,84 @@ USA. ;; 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))))))) -(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))) @@ -158,11 +177,23 @@ USA. (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)))) (define (ieee754-binary-hex-string x exponent-bits precision) (receive (base emin emax bias exp-subnormal exp-inf/nan) @@ -179,39 +210,47 @@ USA. (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))) diff --git a/src/runtime/runtime.pkg b/src/runtime/runtime.pkg index bb69b109b..fdd05d79f 100644 --- a/src/runtime/runtime.pkg +++ b/src/runtime/runtime.pkg @@ -5986,4 +5986,26 @@ USA. (export (runtime) eval-r7rs-scode-file eval-r7rs-source - syntax-r7rs-source)) \ No newline at end of file + syntax-r7rs-source)) + +(define-package (runtime ieee754) + (files "ieee754") + (parent (runtime)) + (export () + compose-ieee754 + compose-ieee754-binary128 + compose-ieee754-binary32 + compose-ieee754-binary64 + decompose-ieee754 + decompose-ieee754-binary128 + decompose-ieee754-binary32 + decompose-ieee754-binary64 + ieee754-binary-hex-string + ieee754-binary-parameters + ieee754-binary128-exact? + ieee754-binary128-hex-string + ieee754-binary32-exact? + ieee754-binary32-hex-string + ieee754-binary64-exact? + ieee754-binary64-hex-string + )) \ No newline at end of file diff --git a/tests/check.scm b/tests/check.scm index 11e6c1084..b26219b0a 100644 --- a/tests/check.scm +++ b/tests/check.scm @@ -49,6 +49,7 @@ USA. "microcode/test-keccak" "microcode/test-lookup" "runtime/test-arith" + "runtime/test-ieee754" "runtime/test-binary-port" "runtime/test-bundle" "runtime/test-bytevector" diff --git a/tests/runtime/test-ieee754.scm b/tests/runtime/test-ieee754.scm new file mode 100644 index 000000000..b04ba5405 --- /dev/null +++ b/tests/runtime/test-ieee754.scm @@ -0,0 +1,147 @@ +#| -*-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 + +(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))))