From a3edb54b6207a832521160322386c5e1784cda99 Mon Sep 17 00:00:00 2001 From: Brett Saiki Date: Fri, 6 Dec 2024 17:50:00 -0800 Subject: [PATCH 1/3] try fixing again --- private/gfl.rkt | 4 ++-- private/mpfr.rkt | 51 ++++++++++++++++++++++++++++-------------------- tests/test.rkt | 3 ++- 3 files changed, 34 insertions(+), 24 deletions(-) diff --git a/private/gfl.rkt b/private/gfl.rkt index 298ec48..938d3de 100644 --- a/private/gfl.rkt +++ b/private/gfl.rkt @@ -75,9 +75,9 @@ (gfl-exponent) (gfl-bits))) (define (gfl->ordinal x) - (define sig (- (gfl-bits) (gfl-exponent))) + (define sig (- (gflonum-nb x) (gflonum-ex x))) (define-values (emin emax) (ex->ebounds (gflonum-ex x) sig)) - ((mpfr-eval emin emax sig) (curryr mpfr->ordinal (gfl-exponent) sig) (gflonum-val x))) + ((mpfr-eval emin emax sig) (curryr mpfr->ordinal (gflonum-ex x) sig) (gflonum-val x))) (define (string->gfl x) (define sig (- (gfl-bits) (gfl-exponent))) diff --git a/private/mpfr.rkt b/private/mpfr.rkt index 1280cdc..de922e7 100644 --- a/private/mpfr.rkt +++ b/private/mpfr.rkt @@ -104,46 +104,55 @@ (define (infinite-ordinal es sig) (arithmetic-shift (sub1 (expt 2 es)) (sub1 sig))) -(define (ordinal->mpfr x es sig) - (define bound (invalid-ordinal es sig)) +(define (ordinal->mpfr x es p) + (define bound (invalid-ordinal es p)) (unless (< (- bound) x bound) (error 'ordinal->mpfr "ordinal out of bounds ~a" x)) (let loop ([x x]) (cond [(zero? x) 0.bf] [(negative? x) (bf- (loop (- x)))] - [(> x (infinite-ordinal es sig)) +nan.bf] - [(= x (infinite-ordinal es sig)) +inf.bf] + [(> x (infinite-ordinal es p)) +nan.bf] + [(= x (infinite-ordinal es p)) +inf.bf] [else ; non-zero, real number - (define msize (sub1 sig)) (define expmin (sub1 (mpfr-get-emin))) - (define ebits (arithmetic-shift x (- msize))) - (define mbits (bitwise-and x (sub1 (expt 2 msize)))) + (define mask (sub1 (expt 2 (sub1 p)))) + (define mbits (bitwise-and x mask)) + (define ebits (arithmetic-shift x (- (sub1 p)))) (cond [(zero? ebits) ; subnormal (bf mbits expmin)] [else ; normal number - (define c (+ (expt 2 msize) mbits)) + (define c (+ (expt 2 (sub1 p)) mbits)) (define exp (+ (sub1 ebits) expmin)) (bf c exp)])]))) -(define (mpfr->ordinal x es sig) +(define (mpfr->ordinal x es p) (let loop ([x x]) (cond [(bfzero? x) 0] [(bfnegative? x) (- (loop (bf- x)))] - [(bfnan? x) (add1 (infinite-ordinal es sig))] - [(bfinfinite? x) (infinite-ordinal es sig)] + [(bfnan? x) (add1 (infinite-ordinal es p))] + [(bfinfinite? x) (infinite-ordinal es p)] [else - (define-values (c exp) (bigfloat->sig+exp x)) - (define e (+ exp (bigfloat-precision x) -1)) + ; format constants + ; with subnormalization MPFR's emin is not what you think it is (define expmin (sub1 (mpfr-get-emin))) - (define emin (+ expmin sig -1)) + (define emin (+ expmin p -1)) + ; extract fields + ; per MPFR documentation, `mpfr_get_z_2exp(x)` extracts the integer + ; using the (full) initialization precision of `x`. + (define-values (c exp) (bigfloat->sig+exp x)) + (define e (+ exp (sub1 p))) (cond - [(< e emin) ; subnormal - (define shift (- exp expmin)) - (arithmetic-shift c shift)] - [else ; normal - (define ebits (add1 (- exp expmin))) - (define mbits (bitwise-and c (sub1 (expt 2 (sub1 sig))))) - (+ (arithmetic-shift ebits (sub1 sig)) mbits)])]))) + [(< e emin) + ; subnormal number + ; since `x` is fully normalized `exp < expmin` + (define shift (- expmin exp)) + (arithmetic-shift c (- shift))] + [else + ; normal number + (define mask (sub1 (expt 2 (sub1 p)))) + (define mbits (bitwise-and mask c)) + (define ebits (add1 (- e emin))) + (+ (arithmetic-shift ebits (sub1 p)) mbits)])]))) diff --git a/tests/test.rkt b/tests/test.rkt index f8755eb..fb18e57 100644 --- a/tests/test.rkt +++ b/tests/test.rkt @@ -101,7 +101,8 @@ ; (check-equal? (gfl- 0.gfl) -0.gfl) ; (check-equal? (gfl- -0.gfl) 0.gfl) -(for ([bits '(13 15 29 27 43 75 139 14 16 18 22 24 140)]) +; (for ([bits '(13 15 29 27 43 75 139 14 16 18 22 24 140)]) +(for ([bits '(13)]) (parameterize ([gfl-bits bits]) ;; +max.gfl/-max.gfl ; (check-equal? (gfl- +max.gfl) -max.gfl) From a25242c9fb1ab02323e9b19fc628fd3dfa5a1358 Mon Sep 17 00:00:00 2001 From: Brett Saiki Date: Fri, 6 Dec 2024 18:46:47 -0800 Subject: [PATCH 2/3] fix precision issues --- private/gfl.rkt | 85 ++++++++++++++++++++++++------------------------- 1 file changed, 41 insertions(+), 44 deletions(-) diff --git a/private/gfl.rkt b/private/gfl.rkt index 938d3de..d20b3b1 100644 --- a/private/gfl.rkt +++ b/private/gfl.rkt @@ -279,66 +279,63 @@ ;;;;;;;;;;;;;;;;;;; Variadic operators ;;;;;;;;;;;;;;;; -(define (mpfrv+ emin emax sig . xs) - (cond - [(null? xs) 0.bf] - [else - (define xs1 (cdr xs)) - (cond - [(null? xs1) (car xs)] - [else - (define xs2 (cdr xs1)) - (cond - [(null? xs2) ((mpfr-eval emin emax sig) mpfr-add (car xs) (car xs1))] - [else ((mpfr-eval emin emax sig) mpfr-sum xs)])])])) - -(define (mpfrv- emin emax sig x . xs) - (cond - [(null? xs) (mpfr-neg x)] - [(null? (cdr xs)) ((mpfr-eval emin emax sig) mpfr-sub x (car xs))] - [else (mpfr-neg (apply mpfrv+ emin emax sig (mpfr-neg x) xs))])) - -(define (mpfrv* emin emax sig . xs) - (cond - [(null? xs) 1.bf] - [else - (let loop ([x (car xs)] [xs (cdr xs)]) - (cond - [(null? xs) x] - [else (loop ((mpfr-eval emin emax sig) mpfr-mul x (car xs)) (cdr xs))]))])) - -(define (mpfrv/ emin emax sig x . xs) - (cond - [(null? xs) ((mpfr-eval emin emax sig) mpfr-div 1.bf x)] - [else - (let loop ([x x] [xs xs]) - (cond - [(null? xs) x] - [else (loop ((mpfr-eval emin emax sig) mpfr-div x (car xs)) (cdr xs))]))])) +(define (mpfrv+ . xs) + (match xs + [(list) 0.bf] + [(list x) (bfcopy x)] + [(list x y rest ...) + (for/fold ([sum (mpfr-add x y)]) ([z (in-list rest)]) + (mpfr-add sum z))])) + +(define (mpfrv- x . xs) + (match xs + [(list) (mpfr-neg x)] + [(list y rest ...) + (for/fold ([diff (mpfr-sub x y)]) ([z (in-list rest)]) + (mpfr-sub diff z))])) + +(define (mpfrv* . xs) + (match xs + [(list) 1.bf] + [(list x) (bfcopy x)] + [(list x y rest ...) + (for/fold ([sum (mpfr-mul x y)]) ([z (in-list rest)]) + (mpfr-mul sum z))])) + +(define (mpfrv/ x . xs) + (match xs + [(list) (mpfr-div 1.bf x)] + [(list y rest ...) + (for/fold ([quo (mpfr-div x y)]) ([z (in-list rest)]) + (mpfr-div quo z))])) (define (gfl+ . args) (define sig (- (gfl-bits) (gfl-exponent))) (define-values (emin emax) (ex->ebounds (gfl-exponent) sig)) - (gflonum (apply mpfrv+ emin emax sig (map gflonum-val args)) - (gfl-exponent) (gfl-bits))) + (gflonum (apply (mpfr-eval emin emax sig) mpfrv+ (map gflonum-val args)) + (gfl-exponent) + (gfl-bits))) (define (gfl- head . rest) (define sig (- (gfl-bits) (gfl-exponent))) (define-values (emin emax) (ex->ebounds (gfl-exponent) sig)) - (gflonum (apply mpfrv- emin emax sig (gflonum-val head) (map gflonum-val rest)) - (gfl-exponent) (gfl-bits))) + (gflonum (apply (mpfr-eval emin emax sig) mpfrv- (gflonum-val head) (map gflonum-val rest)) + (gfl-exponent) + (gfl-bits))) (define (gfl* . args) (define sig (- (gfl-bits) (gfl-exponent))) (define-values (emin emax) (ex->ebounds (gfl-exponent) sig)) - (gflonum (apply mpfrv* emin emax sig (map gflonum-val args)) - (gfl-exponent) (gfl-bits))) + (gflonum (apply (mpfr-eval emin emax sig) mpfrv* (map gflonum-val args)) + (gfl-exponent) + (gfl-bits))) (define (gfl/ head . rest) (define sig (- (gfl-bits) (gfl-exponent))) (define-values (emin emax) (ex->ebounds (gfl-exponent) sig)) - (gflonum (apply mpfrv/ emin emax sig (gflonum-val head) (map gflonum-val rest)) - (gfl-exponent) (gfl-bits))) + (gflonum (apply (mpfr-eval emin emax sig) mpfrv/ (gflonum-val head) (map gflonum-val rest)) + (gfl-exponent) + (gfl-bits))) (define (gflmax2 x y) (cond From 6143d0bb5ef1ec3751c71a64cedb096a94ecc82e Mon Sep 17 00:00:00 2001 From: Brett Saiki Date: Mon, 9 Dec 2024 09:29:51 -0800 Subject: [PATCH 3/3] cleanup --- private/ffi.rkt | 11 ----------- private/gfl.rkt | 4 ++-- private/mpfr.rkt | 1 - 3 files changed, 2 insertions(+), 14 deletions(-) diff --git a/private/ffi.rkt b/private/ffi.rkt index ee137d0..5c1d065 100644 --- a/private/ffi.rkt +++ b/private/ffi.rkt @@ -146,17 +146,6 @@ (mpfr-subnormalize r t (mpfr-rounding-mode))) r) -(define mpfr-sum-fun - (get-mpfr-fun 'mpfr_sum (_fun _mpfr-pointer (_list i _mpfr-pointer) _ulong _rnd_t -> _int))) - -(define (mpfr-sum xs) - (define r (bf 0)) - (define t (mpfr-sum-fun r xs (length xs) (mpfr-rounding-mode))) - (mpfr-check-range r 0 (mpfr-rounding-mode)) - (when (mpfr-subnormalize?) - (mpfr-subnormalize r t (mpfr-rounding-mode))) - r) - (define (mpfr-set x) (define v (if (bigfloat? x) (bfcopy x) (bf x))) (mpfr-check-range v 0 (mpfr-rounding-mode)) diff --git a/private/gfl.rkt b/private/gfl.rkt index d20b3b1..296a209 100644 --- a/private/gfl.rkt +++ b/private/gfl.rkt @@ -299,8 +299,8 @@ [(list) 1.bf] [(list x) (bfcopy x)] [(list x y rest ...) - (for/fold ([sum (mpfr-mul x y)]) ([z (in-list rest)]) - (mpfr-mul sum z))])) + (for/fold ([prod (mpfr-mul x y)]) ([z (in-list rest)]) + (mpfr-mul prod z))])) (define (mpfrv/ x . xs) (match xs diff --git a/private/mpfr.rkt b/private/mpfr.rkt index de922e7..198965a 100644 --- a/private/mpfr.rkt +++ b/private/mpfr.rkt @@ -11,7 +11,6 @@ mpfr-lgamma mpfr-jn mpfr-yn - mpfr-sum mpfr-set mpfr-set-ebounds!)