From 6695286ffb887a994f3b06956625ebae02861c7b Mon Sep 17 00:00:00 2001 From: Stephen Adams Date: Thu, 10 Jul 1997 09:16:23 +0000 Subject: [PATCH] Bug fixes: . "1/" returns #F instead of signalling divide by zero. . "1.2345e-306" and "1e-400" with FLONUM-PARSER-FAST? true no longer signal arithmetic errors. . STRING->NUMBER and SUBSTRING->NUMBER now check their arguments Other changes: FINISH-REAL now exploits the fact that some integers have exact floating point representations and so flonum arithmetic can be used to compute the correct result. Note that to maintain correctness, i*10^-e must be calculated as i/10^e. The results are a modest improvement on free-formatted random flonums (typically 16 digit) in the range 1e-7 to 1e35: Times for 10000 flonums generated by (random range) range 1e-100 1e-10 1e-5 1. 1e5 1e10 1e35 1e40 old 10180 8200 8210 7490 5610 3970 3990 4030 new 10280 8350 4690 4390 3930 3230 3180 4170 With FLONUM-PARSER-FAST? true, using the division trick gives better error rates over a larger range. Number of errors in 10000 random flonums range 1e-20 1e-10 1e-6 1. 1e10 1e20 1e35 1e40 old 3573 2400 3309 2798 3025 685 722 2730 new 2594 2778 1907 1010 761 685 772 2730 --- v7/src/runtime/numpar.scm | 111 ++++++++++++++++++++++++++++++++------ 1 file changed, 96 insertions(+), 15 deletions(-) diff --git a/v7/src/runtime/numpar.scm b/v7/src/runtime/numpar.scm index 3052c4782..9de3dad06 100644 --- a/v7/src/runtime/numpar.scm +++ b/v7/src/runtime/numpar.scm @@ -1,6 +1,6 @@ #| -*-Scheme-*- -$Id: numpar.scm,v 14.12 1997/05/02 05:46:25 cph Exp $ +$Id: numpar.scm,v 14.13 1997/07/10 09:16:23 adams Exp $ Copyright (c) 1989-97 Massachusetts Institute of Technology @@ -38,11 +38,23 @@ MIT in each case. |# (declare (usual-integrations)) (define (string->number string #!optional radix) + (if (not (string? string)) + (error:wrong-type-argument string "string" 'STRING->NUMBER)) (parse-number string 0 (string-length string) (if (default-object? radix) #f radix) 'STRING->NUMBER)) (define (substring->number string start end #!optional radix) + (if (not (string? string)) + (error:wrong-type-argument string "string" 'SUBSTRING->NUMBER)) + (if (not (index-fixnum? start)) + (error:wrong-type-argument start "string index" 'SUBSTRING->NUMBER)) + (if (not (index-fixnum? end)) + (error:wrong-type-argument end "string index" 'SUBSTRING->NUMBER)) + (if (not (< end (string-length string))) + (error:bad-range-argument end 'SUBSTRING->NUMBER)) + (if (not (<= start end)) + (error:bad-range-argument start 'SUBSTRING->NUMBER)) (parse-number string start end (if (default-object? radix) #f radix) 'SUBSTRING->NUMBER)) @@ -125,8 +137,9 @@ MIT in each case. |# (let ((char (string-ref string start)) (start+1 (fix:+ start 1))) (cond ((char=? #\/ char) - (parse-denominator-1 string start+1 end - integer exactness radix sign)) + (and (fix:< start+1 end) + (parse-denominator-1 string start+1 end + integer exactness radix sign))) ((char=? #\. char) (and (fix:= radix 10) (if sharp? @@ -288,20 +301,80 @@ MIT in each case. |# (define (finish-rational numerator denominator exactness sign) ;; State: result is rational, apply exactness and sign. (finish (/ numerator denominator) exactness sign)) + +;; (finish-real integer exponent exactness sign) +;; +;; magnitude is (* INTEGER (EXPT 10 EXPONENT)) +;; +;; In the general case for an inexact result, to obtain a correctly +;; rounded result, it is necessary to work with exact or high +;; precision numbers and convert to the rounded result at the last +;; moment. +;; +;; Sometimes flonum arithmetic is sufficient to obtain a correct result. +;; This is true when all the operations are known, by properties of +;; the numbers they operate on, to give exact results, except possibly +;; for the final operation which must then round correctly. +;; +;; Certain integers can be represented exactly by floating point numbers, +;; for example, IEEE 64 bit fp numbers can represent the integers 0 +;; through 9007199254740991 (lets call these floating point integers), +;; and powers of 10 from 10^0 up to 10^22 (because 5^22 = +;; 2384185791015625 < 9007199254740991). +;; +;; This means that all 15 and fewer digit numbers and 90% of 16 digit +;; numbers with relatively small exponents can be converted correctly +;; using flonum arithmetic. +;; +;; (INTEGER->FLONUM N #b01) acts as both a conversion and a predicate for +;; integers that are also floating point integers. (It might be +;; useful to have an extra flag that tests for N being a floating +;; point integer scaled by a power of two, e.g. 10^20.) +;; +;; Reciprocals of powers of 10 cannot be represented exactly as floating +;; point numbers because 1/10 is a continued fraction in binary. +;; Instead of +;; (* INTEGER (EXPT 10 EXPONENT)) +;; we compute +;; (/ INTEGER (EXPT 10 (- EXPONENT))) +;; This method also benfits accuracy when FLONUM-PARSER-FAST? is true and +;; the reciprocal is exact. + +(define exact-flonum-powers-of-10) ; a vector, i -> 10.^i (define (finish-real integer exponent exactness sign) ;; State: result is integer, apply exactness and sign. - (if (and (or (eq? 'INEXACT exactness) (eq? 'IMPLICIT-INEXACT exactness)) - flonum-parser-fast?) - (if (eq? #\- sign) - (flo:- 0. - (flo:* (int:->flonum integer) - (flo:expt 10. (int:->flonum exponent)))) - (flo:* (int:->flonum integer) - (flo:expt 10. (int:->flonum exponent)))) - (apply-exactness exactness - (* (apply-sign sign integer) - (expt 10 exponent))))) + + (define (high-precision-method) + (apply-exactness exactness + (* (apply-sign sign integer) + (expt 10 exponent)))) + + (if (or (eq? 'INEXACT exactness) (eq? 'IMPLICIT-INEXACT exactness)) + (let ((abs-exponent (if (< exponent 0) (- exponent) exponent)) + (powers-of-10 exact-flonum-powers-of-10)) + (define-integrable (finish-flonum x power-of-10) + (if (eq? #\- sign) + (if (eq? exponent abs-exponent) + (flo:- 0. (flo:* x power-of-10)) + (flo:- 0. (flo:/ x power-of-10))) + (if (eq? exponent abs-exponent) + (flo:* x power-of-10) + (flo:/ x power-of-10)))) + (cond ((and flonum-parser-fast? + (<= abs-exponent 308)) ; this aught to be defined somewhere + (if (< abs-exponent (vector-length powers-of-10)) + (finish-flonum (int:->flonum integer) + (vector-ref powers-of-10 abs-exponent)) + (finish-flonum (int:->flonum integer) + (flo:expt 10. (int:->flonum abs-exponent))))) + ((and (< abs-exponent (vector-length powers-of-10)) + ((ucode-primitive integer->flonum 2) integer #b1)) + => (lambda (exact-flonum-integer) + (finish-flonum exact-flonum-integer + (vector-ref powers-of-10 abs-exponent)))) + (else (high-precision-method)))) + (high-precision-method))) (define flonum-parser-fast? #f) @@ -330,4 +403,12 @@ MIT in each case. |# (or (char=? #\+ char) (char=? #\- char))) (define-integrable (i? char) - (or (char=? #\i char) (char=? #\I char))) \ No newline at end of file + (or (char=? #\i char) (char=? #\I char))) + +(define (initialize-package!) + (set! exact-flonum-powers-of-10 + (let loop ((i 0) (power 1) (powers '())) + (if (= (inexact->exact (exact->inexact power)) power) + (loop (+ i 1) (* power 10) (cons (exact->inexact power) powers)) + (list->vector (reverse! powers))))) + unspecific) -- 2.25.1