Fixes bug 180991: round-off error in bignum->flonum

This commit is contained in:
Abdulaziz Ghuloum 2008-01-12 20:52:23 -05:00
parent 2dc4542148
commit d9cdcb8959
5 changed files with 162 additions and 198 deletions

View File

@ -429,201 +429,8 @@
sin cos tan asin acos atan sqrt truncate fltruncate sin cos tan asin acos atan sqrt truncate fltruncate
flmax random)) flmax random))
(define (bignum->flonum x)
(module (bignum->flonum) (foreign-call "ikrt_bignum_to_flonum" x 0 ($make-flonum)))
; sbe f6 f5 f4 f3 f2 f1 f0
;SEEEEEEE|EEEEmmmm|mmmmmmmm|mmmmmmmm|mmmmmmmm|mmmmmmmm|mmmmmmmm|mmmmmmmm
; | | | | | | |
; v0 v1 v2 v3 v4 v5 v6 v7
(define ($flonum pos? e f6 f5 f4 f3 f2 f1 f0)
(let ([be (fx+ e 1075)])
(let ([v ($make-flonum)])
(cond
[(fx< be 2047)
(let ([sbe (if pos? be (fxlogor be (fxsll 1 11)))])
($flonum-set! v 0 (fxsra sbe 4))
($flonum-set! v 1 (fxlogor (fxsll sbe 4) (fxlogand f6 #b1111)))
($flonum-set! v 2 f5)
($flonum-set! v 3 f4)
($flonum-set! v 4 f3)
($flonum-set! v 5 f2)
($flonum-set! v 6 f1)
($flonum-set! v 7 f0))]
[else ;;; inf
(let ([sbe (if pos? 2047 (fxlogor 2047 (fxsll 1 11)))])
($flonum-set! v 0 (fxsra sbe 4))
($flonum-set! v 1 (fxsll sbe 4))
($flonum-set! v 2 0)
($flonum-set! v 3 0)
($flonum-set! v 4 0)
($flonum-set! v 5 0)
($flonum-set! v 6 0)
($flonum-set! v 7 0))])
v)))
(define ($flonum/c0 pos? e f6 f5 f4 f3 f2 f1 f0 c)
(define ($fxeven? x)
(fxzero? (fxlogand x 1)))
(define-syntax cond*
(syntax-rules (else)
[(_ [test conseq] [else val])
(if test conseq val)]
[(_ [test conseq] [var val] rest ...)
(if test conseq (let ([var val]) (cond* rest ...)))]))
(cond*
[($fxeven? c) ($flonum pos? e f6 f5 f4 f3 f2 f1 f0)]
[f0 (fx+ (fxlogand f0 255) 1)]
[(fx< f0 256) ($flonum pos? e f6 f5 f4 f3 f2 f1 f0)]
[f1 (fx+ (fxlogand f1 255) 1)]
[(fx< f1 256) ($flonum pos? e f6 f5 f4 f3 f2 f1 0)]
[f2 (fx+ (fxlogand f2 255) 1)]
[(fx< f2 256) ($flonum pos? e f6 f5 f4 f3 f2 0 0)]
[f3 (fx+ (fxlogand f3 255) 1)]
[(fx< f3 256) ($flonum pos? e f6 f5 f4 f3 0 0 0)]
[f4 (fx+ (fxlogand f4 255) 1)]
[(fx< f4 256) ($flonum pos? e f6 f5 f4 0 0 0 0)]
[f5 (fx+ (fxlogand f5 255) 1)]
[(fx< f5 256) ($flonum pos? e f6 f5 0 0 0 0 0)]
[f6 (fx+ (fxlogand f6 #b1111) 1)]
[(fx< f6 16) ($flonum pos? e f6 0 0 0 0 0 0)]
[else ($flonum pos? (+ e 1) 0 0 0 0 0 0 0)]))
(define ($flonum/aux pos? e b7 b6 b5 b4 b3 b2 b1 b0)
(cond
[(fx>= b7 #x80)
($flonum/c0 pos? (fx+ e 3)
(fxsra b7 3)
(fxlogor (fxsll b7 5) (fxsra b6 3))
(fxlogor (fxsll b6 5) (fxsra b5 3))
(fxlogor (fxsll b5 5) (fxsra b4 3))
(fxlogor (fxsll b4 5) (fxsra b3 3))
(fxlogor (fxsll b3 5) (fxsra b2 3))
(fxlogor (fxsll b2 5) (fxsra b1 3))
(fxsra b1 2))]
[(fx>= b7 #x40)
($flonum/c0 pos? (fx+ e 2)
(fxsra b7 2)
(fxlogor (fxsll b7 6) (fxsra b6 2))
(fxlogor (fxsll b6 6) (fxsra b5 2))
(fxlogor (fxsll b5 6) (fxsra b4 2))
(fxlogor (fxsll b4 6) (fxsra b3 2))
(fxlogor (fxsll b3 6) (fxsra b2 2))
(fxlogor (fxsll b2 6) (fxsra b1 2))
(fxsra b1 1))]
[(fx>= b7 #x20)
($flonum/c0 pos? (fx+ e 1)
(fxsra b7 1)
(fxlogor (fxsll b7 7) (fxsra b6 1))
(fxlogor (fxsll b6 7) (fxsra b5 1))
(fxlogor (fxsll b5 7) (fxsra b4 1))
(fxlogor (fxsll b4 7) (fxsra b3 1))
(fxlogor (fxsll b3 7) (fxsra b2 1))
(fxlogor (fxsll b2 7) (fxsra b1 1))
b1)]
[(fx>= b7 #x10)
($flonum/c0 pos? e b7 b6 b5 b4 b3 b2 b1
(fxsra b0 7))]
[(fx>= b7 #x08)
($flonum/c0 pos? (fx- e 1)
(fxlogor (fxsll b7 1) (fxsra b6 7))
(fxlogor (fxsll b6 1) (fxsra b5 7))
(fxlogor (fxsll b5 1) (fxsra b4 7))
(fxlogor (fxsll b4 1) (fxsra b3 7))
(fxlogor (fxsll b3 1) (fxsra b2 7))
(fxlogor (fxsll b2 1) (fxsra b1 7))
(fxlogor (fxsll b1 1) (fxsra b0 7))
(fxsra b0 6))]
[(fx>= b7 #x04)
($flonum/c0 pos? (fx- e 2)
(fxlogor (fxsll b7 2) (fxsra b6 6))
(fxlogor (fxsll b6 2) (fxsra b5 6))
(fxlogor (fxsll b5 2) (fxsra b4 6))
(fxlogor (fxsll b4 2) (fxsra b3 6))
(fxlogor (fxsll b3 2) (fxsra b2 6))
(fxlogor (fxsll b2 2) (fxsra b1 6))
(fxlogor (fxsll b1 2) (fxsra b0 6))
(fxsra b0 5))]
[(fx>= b7 #x02)
($flonum/c0 pos? (fx- e 3)
(fxlogor (fxsll b7 3) (fxsra b6 5))
(fxlogor (fxsll b6 3) (fxsra b5 5))
(fxlogor (fxsll b5 3) (fxsra b4 5))
(fxlogor (fxsll b4 3) (fxsra b3 5))
(fxlogor (fxsll b3 3) (fxsra b2 5))
(fxlogor (fxsll b2 3) (fxsra b1 5))
(fxlogor (fxsll b1 3) (fxsra b0 5))
(fxsra b0 4))]
[(fx>= b7 #x01)
($flonum/c0 pos? (fx- e 4)
(fxlogor (fxsll b7 4) (fxsra b6 4))
(fxlogor (fxsll b6 4) (fxsra b5 4))
(fxlogor (fxsll b5 4) (fxsra b4 4))
(fxlogor (fxsll b4 4) (fxsra b3 4))
(fxlogor (fxsll b3 4) (fxsra b2 4))
(fxlogor (fxsll b2 4) (fxsra b1 4))
(fxlogor (fxsll b1 4) (fxsra b0 4))
(fxsra b0 3))]
[else (die '$float/aux "BUG: invalid b7" b7)]))
(define (bignum->flonum x)
(define (bignum/4->flonum x)
($flonum/aux ($bignum-positive? x) -24
($bignum-byte-ref x 3)
($bignum-byte-ref x 2)
($bignum-byte-ref x 1)
($bignum-byte-ref x 0)
0 0 0 0))
(define (bignum/8->flonum x)
;;; bignum: [b0 b1 b2 b3 b4 b5 b6 b7]
(let ([b0 ($bignum-byte-ref x 0)]
[b1 ($bignum-byte-ref x 1)]
[b2 ($bignum-byte-ref x 2)]
[b3 ($bignum-byte-ref x 3)]
[b4 ($bignum-byte-ref x 4)]
[b5 ($bignum-byte-ref x 5)]
[b6 ($bignum-byte-ref x 6)]
[b7 ($bignum-byte-ref x 7)]
[pos? ($bignum-positive? x)])
(if (fx= b7 0)
(if (fx= b6 0)
(if (fx= b5 0)
(if (fx= b4 0)
(die 'bignum8->flonum "malformed bignum")
($flonum/aux pos? -16 b4 b3 b2 b1 b0 0 0 0))
($flonum/aux pos? -8 b5 b4 b3 b2 b1 b0 0 0))
($flonum/aux pos? 0 b6 b5 b4 b3 b2 b1 b0 0))
($flonum/aux pos? 8 b7 b6 b5 b4 b3 b2 b1 b0))))
(define (bignum/n->flonum x bytes)
(define (aux x b7 bytes)
($flonum/aux ($bignum-positive? x) (+ (* bytes 8) -48)
b7
($bignum-byte-ref x (fx- bytes 1))
($bignum-byte-ref x (fx- bytes 2))
($bignum-byte-ref x (fx- bytes 3))
($bignum-byte-ref x (fx- bytes 4))
($bignum-byte-ref x (fx- bytes 5))
($bignum-byte-ref x (fx- bytes 6))
($bignum-byte-ref x (fx- bytes 7))))
;;; bignum: [b0 b1 b2 b3 ... b_{bytes-1}]
(let* ([bytes (fxsub1 bytes)] [bn ($bignum-byte-ref x bytes)])
(if (fx= bn 0)
(let* ([bytes (fxsub1 bytes)] [bn ($bignum-byte-ref x bytes)])
(if (fx= bn 0)
(let* ([bytes (fxsub1 bytes)] [bn ($bignum-byte-ref x bytes)])
(if (fx= bn 0)
(let* ([bytes (fxsub1 bytes)] [bn ($bignum-byte-ref x bytes)])
(if (fx= bn 0)
(die 'bignum/n->flonum "malformed bignum")
(aux x bn bytes)))
(aux x bn bytes)))
(aux x bn bytes)))
(aux x bn bytes))))
(unless (bignum? x)
(die 'bignum->flonum "not a bignum" x))
(let ([bytes ($bignum-size x)])
(case bytes
[(4) (bignum/4->flonum x)]
[(8) (bignum/8->flonum x)]
[else (bignum/n->flonum x bytes)]))))
;;; (define (ratnum->flonum x) ;;; (define (ratnum->flonum x)

View File

@ -1 +1 @@
1338 1339

View File

@ -57,6 +57,7 @@
(test-strings) (test-strings)
(test-exact-integer-sqrt) (test-exact-integer-sqrt)
(test-bignum-to-flonum) (test-bignum-to-flonum)
(test-bignum->flonum)
(test-string-to-number) (test-string-to-number)
(test-div-and-mod) (test-div-and-mod)
(test-bignums) (test-bignums)

View File

@ -1,6 +1,6 @@
(library (tests bignum-to-flonum) (library (tests bignum-to-flonum)
(export test-bignum-to-flonum) (export test-bignum-to-flonum test-bignum->flonum)
(import (ikarus) (tests framework)) (import (ikarus) (tests framework))
(define (t x s) (define (t x s)
(let ([fl (format "~a" (exact->inexact x))]) (let ([fl (format "~a" (exact->inexact x))])
@ -13,6 +13,41 @@
(define-tests name (define-tests name
[(lambda (x) (string=? str (number->string x))) (exact->inexact num)] [(lambda (x) (string=? str (number->string x))) (exact->inexact num)]
...)])) ...)]))
(define (testnum x)
(define precision 53)
(assert (bignum? x))
(let ([fl (inexact x)])
(let ([n (if (> x 0) x (- x))])
(let ([bits (bitwise-length n)])
(printf "bits = ~s\n" bits)
(cond
[(<= bits precision)
(unless (= x (exact fl))
(error #f "should be exactly equal" x fl (exact fl)))]
[else
(let ([hi53 (sra n (- bits precision))]
[lo (bitwise-and n (- (sll 1 (- bits precision)) 1))]
[breakpoint (sll 1 (- bits precision 1))])
(assert (= n (+ lo (sll hi53 (- bits precision)))))
(let ([fl2
(cond
[(or (< lo breakpoint)
(and (= lo breakpoint) (even? hi53)))
(* (inexact hi53) (sll 1 (- bits precision)))]
[else
(* (inexact (+ hi53 1)) (sll 1 (- bits precision)))])])
(let ([fl2 (if (> x 0) fl2 (* fl2 -1))])
(printf "x=~s fl=~s\n" x fl)
(unless (fl=? fl fl2)
(error #f "should be equal" x fl fl2)))))])))))
(define (test-pos-neg x)
(testnum x)
(testnum (- x)))
(test* test-bignum-to-flonum (test* test-bignum-to-flonum
(1000000000 "1e9") (1000000000 "1e9")
(2000000000 "2e9") (2000000000 "2e9")
@ -50,5 +85,40 @@
(100000000000000000000000000000000000000000000000 "1e47") (100000000000000000000000000000000000000000000000 "1e47")
(-1000000000000000000000 "-1e21") (-1000000000000000000000 "-1e21")
(-100000000000000000000000000000 "-1e29") (-100000000000000000000000000000 "-1e29")
(-100000000000000000000000000000000000000000000000 "-1e47"))) (-100000000000000000000000000000000000000000000000 "-1e47"))
(define (test-bignum->flonum)
(test-pos-neg 34872389478)
(test-pos-neg 34872389479)
(test-pos-neg 3487238948347878)
(test-pos-neg 3487238948347879)
(test-pos-neg 5487238948347878)
(test-pos-neg 5487238948347879)
(test-pos-neg 543877238948347878)
(test-pos-neg 543877238948347879)
(test-pos-neg 5438748878948347878)
(test-pos-neg 5438748878948347879)
(test-pos-neg 13874887238948347878)
(test-pos-neg 13874887238948347879)
(test-pos-neg 543874887238948347878)
(test-pos-neg 543874887238948347879)
(test-pos-neg 5433847834874887238948347878)
(test-pos-neg 5433847834874887238948347879)
(test-pos-neg 329847892374892374895433847834874887238948347878)
(test-pos-neg 329847892374892374895433847834874887238948347879)
(test-pos-neg
13407807929942598588139732355608757972494524375225679733981068131349151486565474898751136354405850399729303719974268319295398132445078977825297784408899585)
(test-pos-neg
13407807929942598588139732355608757972494524375225679733981068131349151486565474898751136354405850399729303719974268319295398132445078977825297784408899584)
(test-pos-neg
13407807929942598588139732355608757972494524375225679733981068131349151486565474898751136354405850399729303719974268319295398132445078977825297784408899586)
(test-pos-neg
1340780792994259858813973235560875797249452437522567973398106813134915148656547489875113635440585039972930371997426831929539813244507897782529778440889958413407807929942598588139732355608757972494524375225679733981068131349151486565474898751136354405850399729303719974268319295398132445078977825297784408899584))
)

View File

@ -1746,3 +1746,89 @@ ikrt_fxrandom(ikptr x){
} }
} }
static int
limb_size(mp_limb_t x){
int i = 0;
while(x){
i++;
x = x>>1;
}
return i;
}
static int
all_zeros(mp_limb_t* start, mp_limb_t* end){
while(start <= end){
if(*end) return 0;
end--;
}
return 1;
}
#define PRECISION 53
ikptr
ikrt_bignum_to_flonum(ikptr bn, ikptr more_bits, ikptr fl){
if(mp_bits_per_limb != 32){
fprintf(stderr, "ikarus BUG: bignum_to_flonum only works in 32bit now\n");
exit(-1);
}
ikptr fst = ref(bn, -vector_tag);
long int limb_count = bnfst_limb_count(fst);
mp_limb_t* sp = (mp_limb_t*)(long)(bn+off_bignum_data);
double pos_result;
if(limb_count == 1){
pos_result = sp[0];
} else if (limb_count == 2){
mp_limb_t lo = sp[0];
mp_limb_t hi = sp[1];
pos_result = hi;
pos_result = pos_result * 4294967296.0;
pos_result = pos_result + lo;
} else {
mp_limb_t hi = sp[limb_count-1];
mp_limb_t mi = sp[limb_count-2];
int bc = limb_size(hi);
if(bc < 32){
mp_limb_t lo = sp[limb_count-3];
hi = (hi << (32-bc)) | (mi >> bc);
mi = (mi << (32-bc)) | (lo >> bc);
}
/* now hi has 32 full bits, and mi has 32 full bits */
mp_limb_t mask = ((1<<(64-PRECISION)) - 1);
if((mi & mask) == ((mask+1)>>1)){
/* exactly at break point */
if(((sp[limb_count-3] << (32-bc)) == 0) &&
all_zeros(sp, sp+limb_count-4) &&
(more_bits == 0)){
if(mi & (1<<(64-PRECISION))){
/* odd number, round to even */
mi = mi | mask;
}
} else {
/* round up */
mi = mi | mask;
}
} else if ((mi & mask) > ((mask+1)>>1)){
/* also round up */
mi = mi | mask;
} else {
/* keep it to round down */
}
pos_result = hi;
pos_result = pos_result * 4294967296.0;
pos_result = pos_result + mi;
int bignum_bits = bc + (mp_bits_per_limb * (limb_count-1));
int exponent = bignum_bits - (2 * mp_bits_per_limb);
while(exponent){
pos_result *= 2.0;
exponent -= 1;
}
}
if(bnfst_negative(fst)){
flonum_data(fl) = - pos_result;
} else {
flonum_data(fl) = pos_result;
}
return fl;
}