adding floating point utils for bytevectors

This commit is contained in:
Alex Shinn 2018-12-04 00:43:08 +08:00
parent 11ccfcb5de
commit d513bdc977
4 changed files with 285 additions and 2 deletions

View file

@ -0,0 +1,67 @@
(define-library (chibi bytevector-test)
(export run-tests)
(import (scheme base) (chibi bytevector) (chibi test))
(begin
(define floats
`(0.0 -1.0 #i1/3 1.192092896E-07 ,(+ 1 1.192092896E-07)
1e-23 -1e-23
3.40282346638528860e+38 -3.40282346638528860e+38
1.40129846432481707e-45 -1.40129846432481707e-45
3.14159265358979323846))
(define f32-le
'#u8(#x00 #x00 #x00 #x00 #x00 #x00 #x80 #xbf
#xab #xaa #xaa #x3e #x00 #x00 #x00 #x34
#x01 #x00 #x80 #x3f #x9a #x6d #x41 #x19
#x9a #x6d #x41 #x99 #xff #xff #x7f #x7f
#xff #xff #x7f #xff #x01 #x00 #x00 #x00
#x01 #x00 #x00 #x80 #xdb #x0f #x49 #x40))
(define f64-le
'#u8(#x00 #x00 #x00 #x00 #x00 #x00 #x00 #x00
#x00 #x00 #x00 #x00 #x00 #x00 #xf0 #xbf
#x55 #x55 #x55 #x55 #x55 #x55 #xd5 #x3f
#x68 #x5f #x1c #x00 #x00 #x00 #x80 #x3e
#x00 #x00 #x00 #x20 #x00 #x00 #xf0 #x3f
#x51 #xb2 #x12 #x40 #xb3 #x2d #x28 #x3b
#x51 #xb2 #x12 #x40 #xb3 #x2d #x28 #xbb
#x00 #x00 #x00 #xe0 #xff #xff #xef #x47
#x00 #x00 #x00 #xe0 #xff #xff #xef #xc7
#x00 #x00 #x00 #x00 #x00 #x00 #xa0 #x36
#x00 #x00 #x00 #x00 #x00 #x00 #xa0 #xb6
#x18 #x2d #x44 #x54 #xfb #x21 #x09 #x40))
(define (run-tests)
(test-begin "bytevector")
(test-group "reading"
(do ((ls floats (cdr ls))
(i 0 (+ i 4)))
((null? ls))
(test (car ls) (bytevector-ieee-single-native-ref f32-le i)))
(do ((ls floats (cdr ls))
(i 0 (+ i 8)))
((null? ls))
(test (car ls) (bytevector-ieee-double-native-ref f64-le i))))
(test-group "writing"
(do ((ls floats (cdr ls))
(i 0 (+ i 4)))
((null? ls))
(let ((bv (make-bytevector 4 0)))
(bytevector-ieee-single-native-set! bv 0 (car ls))
(test (bytevector-copy f32-le i (+ i 4)) (values bv))))
(do ((ls floats (cdr ls))
(i 0 (+ i 8)))
((null? ls))
(let ((bv (make-bytevector 8 0)))
(bytevector-ieee-double-native-set! bv 0 (car ls))
(test (bytevector-copy f64-le i (+ i 8)) (values bv)))))
(test-end))))

View file

@ -8,10 +8,29 @@
bytevector-pad-left
integer->bytevector bytevector->integer
integer->hex-string hex-string->integer
bytevector->hex-string hex-string->bytevector)
bytevector->hex-string hex-string->bytevector
bytevector-ieee-single-native-ref
bytevector-ieee-single-native-set!
bytevector-ieee-double-native-ref
bytevector-ieee-double-native-set!
)
(import (scheme base))
(cond-expand
(big-endian
(begin
(define-syntax native-endianness
(syntax-rules () ((_) 'big)))))
(else
(begin
(define-syntax native-endianness
(syntax-rules () ((_) 'little))))))
(cond-expand
((library (srfi 151)) (import (srfi 151)))
((library (srfi 33)) (import (srfi 33)))
(else (import (srfi 60))))
(include "bytevector.scm"))
(include "bytevector.scm")
(cond-expand
;;(chibi
;; (include-shared "ieee-754-native"))
(else
(include "ieee-754.scm"))))

195
lib/chibi/ieee-754.scm Normal file
View file

@ -0,0 +1,195 @@
;;; Copyright (c) 2004-2018 by Alex Shinn.
;; Adapted from SRFI-56.
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; syntax
(define-syntax combine
(syntax-rules ()
((combine) 0)
((combine b1) b1)
((combine b1 b2 b3 ...)
(combine (+ (arithmetic-shift b1 8) b2) b3 ...))))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; reading floating point numbers
;; Inspired by Oleg's implementation from
;; http://okmij.org/ftp/Scheme/reading-IEEE-floats.txt
;; but removes mutations and magic numbers and allows for manually
;; specifying the endianness.
;;
;; See also
;; http://www.cs.auckland.ac.nz/~jham1/07.211/floats.html
;; and
;; http://babbage.cs.qc.edu/courses/cs341/IEEE-754references.html
;; as references to IEEE 754.
(define (bytevector-ieee-single-native-ref bytevector k)
(define (mantissa expn b2 b3 b4)
(case expn ; recognize special literal exponents
((255)
;;(if (zero? (combine b2 b3 b4)) +/0. 0/0.) ; XXXX for SRFI-70
#f)
((0) ; denormalized
(inexact (* (expt 2.0 (- 1 (+ 127 23))) (combine b2 b3 b4))))
(else
(inexact
(* (expt 2.0 (- expn (+ 127 23)))
(combine (+ b2 128) b3 b4)))))) ; hidden bit
(define (exponent b1 b2 b3 b4)
(if (> b2 127) ; 1st bit of b2 is low bit of expn
(mantissa (+ (* 2 b1) 1) (- b2 128) b3 b4)
(mantissa (* 2 b1) b2 b3 b4)))
(define (sign b1 b2 b3 b4)
(if (> b1 127) ; 1st bit of b1 is sign
(cond ((exponent (- b1 128) b2 b3 b4) => -) (else #f))
(exponent b1 b2 b3 b4)))
(let* ((b1 (bytevector-u8-ref bytevector (+ k 0)))
(b2 (bytevector-u8-ref bytevector (+ k 1)))
(b3 (bytevector-u8-ref bytevector (+ k 2)))
(b4 (bytevector-u8-ref bytevector (+ k 3))))
(if (eq? (native-endianness) 'big)
(sign b1 b2 b3 b4)
(sign b4 b3 b2 b1))))
(define (bytevector-ieee-double-native-ref bytevector k)
(define (mantissa expn b2 b3 b4 b5 b6 b7 b8)
(case expn ; recognize special literal exponents
((255) #f) ; won't handle NaN and +/- Inf
((0) ; denormalized
(inexact (* (expt 2.0 (- 1 (+ 1023 52)))
(combine b2 b3 b4 b5 b6 b7 b8))))
(else
(inexact
(* (expt 2.0 (- expn (+ 1023 52)))
(combine (+ b2 16) b3 b4 b5 b6 b7 b8)))))) ; hidden bit
(define (exponent b1 b2 b3 b4 b5 b6 b7 b8)
(mantissa (bitwise-ior (arithmetic-shift b1 4) ; 7 bits
(arithmetic-shift b2 -4)) ; + 4 bits
(bitwise-and b2 #b1111)
b3 b4 b5 b6 b7 b8))
(define (sign b1 b2 b3 b4 b5 b6 b7 b8)
(if (> b1 127) ; 1st bit of b1 is sign
(cond ((exponent (- b1 128) b2 b3 b4 b5 b6 b7 b8) => -)
(else #f))
(exponent b1 b2 b3 b4 b5 b6 b7 b8)))
(let* ((b1 (bytevector-u8-ref bytevector (+ k 0)))
(b2 (bytevector-u8-ref bytevector (+ k 1)))
(b3 (bytevector-u8-ref bytevector (+ k 2)))
(b4 (bytevector-u8-ref bytevector (+ k 3)))
(b5 (bytevector-u8-ref bytevector (+ k 4)))
(b6 (bytevector-u8-ref bytevector (+ k 5)))
(b7 (bytevector-u8-ref bytevector (+ k 6)))
(b8 (bytevector-u8-ref bytevector (+ k 7))))
(if (eq? (native-endianness) 'big)
(sign b1 b2 b3 b4 b5 b6 b7 b8)
(sign b8 b7 b6 b5 b4 b3 b2 b1))))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; writing floating point numbers
;; Underflow rounds down to zero as in IEEE-754, and overflow gets
;; written as +/- Infinity.
;; Break a real number down to a normalized mantissa and exponent.
;; Default base=2, mant-size=23 (52), exp-size=8 (11) for IEEE singles
;; (doubles).
;;
;; Note: This should never be used in practice, since it can be
;; implemented much faster in C. See decode-float in ChezScheme or
;; Gauche.
(define (call-with-mantissa&exponent num base mant-size exp-size proc)
(cond
((negative? num)
(call-with-mantissa&exponent (- num) base mant-size exp-size proc))
((zero? num) (proc 0 0))
(else
(let* ((bot (expt base mant-size))
(top (* base bot)))
(let loop ((n (inexact num)) (e 0))
(cond
((>= n top)
(loop (/ n base) (+ e 1)))
((< n bot)
(loop (* n base) (- e 1)))
(else
(proc (exact (round n)) e))))))))
(define (bytevector-ieee-single-native-set! bytevector k num)
(define (bytes)
(call-with-mantissa&exponent num 2 23 8
(lambda (f e)
(let ((e0 (+ e 127 23)))
(cond
((negative? e0)
(let* ((f1 (exact (round (* f (expt 2 (- e0 1))))))
(b2 (bit-field f1 16 24)) ; mant:16-23
(b3 (bit-field f1 8 16)) ; mant:8-15
(b4 (bit-field f1 0 8))) ; mant:0-7
(list (if (negative? num) 128 0) b2 b3 b4)))
((> e0 255) ; XXXX here we just write infinity
(list (if (negative? num) 255 127) 128 0 0))
(else
(let* ((b0 (arithmetic-shift e0 -1))
(b1 (if (negative? num) (+ b0 128) b0)) ; sign + exp:1-7
(b2 (bitwise-ior
(if (odd? e0) 128 0) ; exp:0
(bit-field f 16 23))) ; + mant:16-23
(b3 (bit-field f 8 16)) ; mant:8-15
(b4 (bit-field f 0 8))) ; mant:0-7
(list b1 b2 b3 b4))))))))
(let ((result (cond
((zero? num) '(0 0 0 0))
((eq? (native-endianness) 'big) (bytes))
(else (reverse (bytes))))))
(bytevector-u8-set! bytevector (+ k 0) (list-ref result 0))
(bytevector-u8-set! bytevector (+ k 1) (list-ref result 1))
(bytevector-u8-set! bytevector (+ k 2) (list-ref result 2))
(bytevector-u8-set! bytevector (+ k 3) (list-ref result 3))))
(define (bytevector-ieee-double-native-set! bytevector k num)
(define (bytes)
(call-with-mantissa&exponent num 2 52 11
(lambda (f e)
(let ((e0 (+ e 1023 52)))
(cond
((negative? e0)
(let* ((f1 (exact (round (* f (expt 2 (- e0 1))))))
(b2 (bit-field f1 48 52))
(b3 (bit-field f1 40 48))
(b4 (bit-field f1 32 40))
(b5 (bit-field f1 24 32))
(b6 (bit-field f1 16 24))
(b7 (bit-field f1 8 16))
(b8 (bit-field f1 0 8)))
(list (if (negative? num) 128 0) b2 b3 b4 b5 b6 b7 b8)))
((> e0 4095) ; infinity
(list (if (negative? num) 255 127) 224 0 0 0 0 0 0))
(else
(let* ((b0 (bit-field e0 4 11))
(b1 (if (negative? num) (+ b0 128) b0))
(b2 (bitwise-ior (arithmetic-shift
(bit-field e0 0 4)
4)
(bit-field f 48 52)))
(b3 (bit-field f 40 48))
(b4 (bit-field f 32 40))
(b5 (bit-field f 24 32))
(b6 (bit-field f 16 24))
(b7 (bit-field f 8 16))
(b8 (bit-field f 0 8)))
(list b1 b2 b3 b4 b5 b6 b7 b8))))))))
(let ((result (cond
((zero? num) '(0 0 0 0 0 0 0 0))
((eq? (native-endianness) 'big) (bytes))
(else (reverse (bytes))))))
(bytevector-u8-set! bytevector (+ k 0) (list-ref result 0))
(bytevector-u8-set! bytevector (+ k 1) (list-ref result 1))
(bytevector-u8-set! bytevector (+ k 2) (list-ref result 2))
(bytevector-u8-set! bytevector (+ k 3) (list-ref result 3))
(bytevector-u8-set! bytevector (+ k 4) (list-ref result 4))
(bytevector-u8-set! bytevector (+ k 5) (list-ref result 5))
(bytevector-u8-set! bytevector (+ k 6) (list-ref result 6))
(bytevector-u8-set! bytevector (+ k 7) (list-ref result 7))))

View file

@ -29,6 +29,7 @@
(rename (srfi 139 test) (run-tests run-srfi-139-tests))
(rename (srfi 151 test) (run-tests run-srfi-151-tests))
(rename (chibi base64-test) (run-tests run-base64-tests))
(rename (chibi bytevector-test) (run-tests run-bytevector-tests))
(rename (chibi crypto md5-test) (run-tests run-md5-tests))
(rename (chibi crypto rsa-test) (run-tests run-rsa-tests))
(rename (chibi crypto sha2-test) (run-tests run-sha2-tests))
@ -89,6 +90,7 @@
(run-srfi-139-tests)
(run-srfi-151-tests)
(run-base64-tests)
(run-bytevector-tests)
(run-doc-tests)
(run-generic-tests)
(run-io-tests)