From 77e27ff729d170f22b0c5f52d816f6d44fc47784 Mon Sep 17 00:00:00 2001 From: Stephen Adams Date: Tue, 8 Jul 1997 06:04:02 +0000 Subject: [PATCH] Speed up RATIO->FLONUM and INT:->INEXACT by using (INTEGER->FLONUM N #b11) when N is small enough that it has an `exact' flonum representation. A prior log message comments on this, but the code does not seem to take advantage. DISCUSSION: Question: why is (INT:->INEXACT N) not simply (INTEGER->FLONUM N #b10) ? (EXACT->INEXACT (expt 2 3000)) used to fail in INTEGER->FLONUM, but now it returns MAXDOUBLE (due to ldexp), but an Infinity or an error would seem better. R4RS says: If an exact argument has no reasonably close inexact equivalent, the a violation of an implementation restriction may be reported. I would read this as NOT returning MAXDOUBLE as, say, 2^3000 is not `reasonably close' to any FP number. A previous log entry says INTEGER->FLONUM does not round reliably. This is because bignum_to_double in "bignum.c" accumulates, which may cause error due to intermediate rounding. Perhaps bignum_to_double should be changed to extract the top 53 bits and explicitly calculate the exponent; another test would be required in place of bignum_fits_in_word_p, which does not (and should not) understand rounding carry. Currently, calling bignum_to_double on (- (expt 2 1024) 1) signals a floating point overflow or returns an infinity depending on which FPU exceptions are enabled. If bignum_to_double was fixed it could be a lot faster than all the current bignum arithmetic. --- v7/src/runtime/arith.scm | 109 ++++++++++++++++++++------------------- 1 file changed, 57 insertions(+), 52 deletions(-) diff --git a/v7/src/runtime/arith.scm b/v7/src/runtime/arith.scm index 0c3ec67de..641d23607 100644 --- a/v7/src/runtime/arith.scm +++ b/v7/src/runtime/arith.scm @@ -1,6 +1,6 @@ #| -*-Scheme-*- -$Id: arith.scm,v 1.42 1997/07/08 01:22:28 adams Exp $ +$Id: arith.scm,v 1.43 1997/07/08 06:04:02 adams Exp $ Copyright (c) 1989-97 Massachusetts Institute of Technology @@ -772,59 +772,64 @@ MIT in each case. |# (int:->inexact q))) (define (ratio->flonum n d) - (let ((n>0 - (lambda (n d) - (let ((k (int:- (integer-length-in-bits n) - (integer-length-in-bits d))) - (p flo:significand-digits-base-2)) - (letrec - ((step1 - (lambda (n d) - ;; (assert (< (expt 2 (- k 1)) (/ n d) (expt 2 (+ k 1)))) - (if (int:negative? k) - (step2 (integer-shift-left n (int:negate k)) d) - (step2 n (integer-shift-left d k))))) - (step2 - (lambda (n d) - ;; (assert (< 1/2 (/ n d) 2)) - (if (int:< n d) - (step3 n d (int:- k p)) - (step3 n (int:* 2 d) (int:- (int:1+ k) p))))) - (step3 - (lambda (n d e) - ;; (assert (and (<= 1/2 (/ n d)) (< (/ n d) 1))) - (let ((n (int:round (integer-shift-left n p) d))) - (if (int:= n int:flonum-integer-limit) - (step4 (int:quotient n 2) (int:1+ e)) - (step4 n e))))) - (step4 - (lambda (n e) - (flo:denormalize (integer->flonum n #b11) e)))) - (step1 n d)))))) - (cond ((fix:zero? n) flo:0) - ((int:positive? n) (n>0 n d)) - (else (flo:negate (n>0 (int:negate n) d)))))) + (define (n>0 n d) + (if (and (int:< n int:flonum-integer-limit) ; integer->flonum `exact'? + (int:< d int:flonum-integer-limit)) ; integer->flonum `exact'? + (flo:/ (integer->flonum n #b11) (integer->flonum d #b11)) ; flo:/ rounds + (let ((k (int:- (integer-length-in-bits n) + (integer-length-in-bits d))) + (p flo:significand-digits-base-2)) + (letrec + ((step1 + (lambda (n d) + ;; (assert (< (expt 2 (- k 1)) (/ n d) (expt 2 (+ k 1)))) + (if (int:negative? k) + (step2 (integer-shift-left n (int:negate k)) d) + (step2 n (integer-shift-left d k))))) + (step2 + (lambda (n d) + ;; (assert (< 1/2 (/ n d) 2)) + (if (int:< n d) + (step3 n d (int:- k p)) + (step3 n (int:* 2 d) (int:- (int:1+ k) p))))) + (step3 + (lambda (n d e) + ;; (assert (and (<= 1/2 (/ n d)) (< (/ n d) 1))) + (let ((n (int:round (integer-shift-left n p) d))) + (if (int:= n int:flonum-integer-limit) + (step4 (int:quotient n 2) (int:1+ e)) + (step4 n e))))) + (step4 + (lambda (n e) + (flo:denormalize (integer->flonum n #b11) e)))) + (step1 n d))))) + + (cond ((fix:zero? n) flo:0) + ((int:positive? n) (n>0 n d)) + (else (flo:negate (n>0 (int:negate n) d))))) (define (int:->inexact n) - (let ((n>0 - (lambda (n) - (let ((e (int:- (integer-length-in-bits n) - flo:significand-digits-base-2)) - (finish - (lambda (n e) - (flo:denormalize (integer->flonum n #b11) e)))) - (cond ((fix:zero? e) - (finish n e)) - ((int:positive? e) - (let ((n (int:round n (integer-shift-left 1 e)))) - (if (int:= n int:flonum-integer-limit) - (finish (int:quotient n 2) (int:1+ e)) - (finish n e)))) - (else - (finish (integer-shift-left n (int:negate e)) e))))))) - (cond ((fix:zero? n) flo:0) - ((int:positive? n) (n>0 n)) - (else (flo:negate (n>0 (int:negate n))))))) + (define (n>0 n) + (if (int:< n int:flonum-integer-limit) ; The flonum is `exact' + (integer->flonum n #b11) + (let ((e (int:- (integer-length-in-bits n) + flo:significand-digits-base-2)) + (finish + (lambda (n e) + (flo:denormalize (integer->flonum n #b11) e)))) + (cond ((fix:zero? e) + (finish n e)) + ((int:positive? e) + (let ((n (int:round n (integer-shift-left 1 e)))) + (if (int:= n int:flonum-integer-limit) + (finish (int:quotient n 2) (int:1+ e)) + (finish n e)))) + (else + (finish (integer-shift-left n (int:negate e)) e)))))) + + (cond ((fix:zero? n) flo:0) + ((int:positive? n) (n>0 n)) + (else (flo:negate (n>0 (int:negate n)))))) (define (flo:significand-digits radix) (cond ((int:= radix 10) -- 2.25.1