suzanne.soy/racket
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Chez Scheme: fix flonum modulo and remainder

Browse Source 
Use fmod() instead of trying to work around imprecision in the naive
algorithm.

Closes #3469

Great catch, Xsmith!
This commit is contained in:
Matthew Flatt 2020-10-27 16:47:10 -06:00
parent fd2ffe3170
commit d70a11236f
4 changed files with 67 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
pkgs/racket-test-core/tests/racket/number.rktl
View File

@ -1110,6 +1110,19 @@
(err/rt-test (modulo 6 -0.0) exn:fail:contract:divide-by-zero?)
(err/rt-test (remainder 6 -0.0) exn:fail:contract:divide-by-zero?)
(let ()
(define (check-rem-mod a b rem mod)
(test rem remainder a b)
(test rem remainder (inexact->exact a) b)
(test rem remainder a (inexact->exact b))
(test mod modulo a b)
(test mod modulo (inexact->exact a) b)
(test mod modulo a (inexact->exact b)))
(check-rem-mod 5.842423430828094e+60 10.0 4.0 4.0)
(check-rem-mod 5.842423430828094e+60 -10.0 4.0 -6.0)
(check-rem-mod -5.842423430828094e+60 10.0 -4.0 6.0)
(check-rem-mod -5.842423430828094e+60 -10.0 -4.0 -4.0))
(define (test-qrm-inf v)
(define iv (exact->inexact v))

4
racket/src/ChezScheme/c/prim5.c
View File

@ -1320,6 +1320,9 @@ extern double log1p();
#endif /* LOG1P */
#endif /* defined(__STDC__) || defined(USE_ANSI_PROTOTYPES) */
static double s_mod PROTO((double x, double y));
static double s_mod(x, y) double x, y; { return fmod(x, y); }
static double s_exp PROTO((double x));
static double s_exp(x) double x; { return exp(x); }
@ -1756,6 +1759,7 @@ void S_prim5_init() {
Sforeign_symbol("(cs)dequeue_scheme_signals", (void *)S_dequeue_scheme_signals);
Sforeign_symbol("(cs)register_scheme_signal", (void *)S_register_scheme_signal);
Sforeign_symbol("(cs)mod", (void *)s_mod);
Sforeign_symbol("(cs)exp", (void *)s_exp);
Sforeign_symbol("(cs)log", (void *)s_log);
Sforeign_symbol("(cs)pow", (void *)s_pow);

40
racket/src/ChezScheme/mats/5_3.ms
View File

@ -2101,9 +2101,28 @@
(fl= (remainder 5.0 2.0) 1.0)
(fl= (remainder -5.0 3.0) -2.0)
(fl= (remainder 5.0 -3.0) 2.0)
(eqv? (remainder -4.0 2.0) 0.0)
(eqv? (remainder 4.0 -2.0) 0.0)
(eqv? (remainder 0 2.0) 0)
(fl= (remainder 5.842423430828094e+60 10) 4.0)
(fl= (remainder 5.842423430828094e+60 10.0) 4.0)
(fl= (remainder 5.842423430828094e+60 -10) 4.0)
(fl= (remainder 5.842423430828094e+60 -10.0) 4.0)
(fl= (remainder -5.842423430828094e+60 10) -4.0)
(fl= (remainder -5.842423430828094e+60 10.0) -4.0)
(fl= (remainder -5.842423430828094e+60 -10) -4.0)
(fl= (remainder -5.842423430828094e+60 -10.0) -4.0)
(fl= (remainder (exact 5.842423430828094e+60) 10.0) 4.0)
(fl= (remainder (exact 5.842423430828094e+60) -10.0) 4.0)
(fl= (remainder (exact -5.842423430828094e+60) 10.0) -4.0)
(fl= (remainder (exact -5.842423430828094e+60) -10.0) -4.0)
(eqv? (remainder (exact 5.842423430828094e+60) 10) 4)
(eqv? (remainder (exact 5.842423430828094e+60) -10) 4)
(eqv? (remainder (exact -5.842423430828094e+60) 10) -4)
(eqv? (remainder (exact -5.842423430828094e+60) -10) -4)
;; following returns incorrect result with naive algorithm,
;; i.e., remainder = (lambda (x,y) (- x (* (quotient x y) y)))
(fl= (remainder 1e194 10.0) 0.0)
(fl= (remainder 1e194 10.0) 8.0)
;; following returns incorrect result in all versions prior to 5.9b
(eq? (remainder (most-negative-fixnum) (- (most-negative-fixnum))) 0)
)
@ -2130,6 +2149,25 @@
(fl= (modulo 5.0 2) 1.0)
(fl= (modulo 5.0 2.0) 1.0)
(fl= (modulo 5.0 2.0) 1.0)
(eqv? (modulo -4.0 2.0) 0.0)
(eqv? (modulo 4.0 -2.0) 0.0)
(eqv? (modulo 0 2.0) 0)
(fl= (modulo 5.842423430828094e+60 10) 4.0)
(fl= (modulo 5.842423430828094e+60 10.0) 4.0)
(fl= (modulo -5.842423430828094e+60 10) 6.0)
(fl= (modulo -5.842423430828094e+60 10.0) 6.0)
(fl= (modulo 5.842423430828094e+60 -10) -6.0)
(fl= (modulo 5.842423430828094e+60 -10.0) -6.0)
(fl= (modulo -5.842423430828094e+60 -10) -4.0)
(fl= (modulo -5.842423430828094e+60 -10.0) -4.0)
(fl= (modulo (exact 5.842423430828094e+60) 10.0) 4.0)
(fl= (modulo (exact -5.842423430828094e+60) 10.0) 6.0)
(fl= (modulo (exact 5.842423430828094e+60) -10.0) -6.0)
(fl= (modulo (exact -5.842423430828094e+60) -10.0) -4.0)
(eqv? (modulo (exact 5.842423430828094e+60) 10) 4)
(eqv? (modulo (exact -5.842423430828094e+60) 10) 6)
(eqv? (modulo (exact 5.842423430828094e+60) -10) -6)
(eqv? (modulo (exact -5.842423430828094e+60) -10) -4)
)
(mat truncate

18
racket/src/ChezScheme/s/5_3.ss
View File

@ -1951,13 +1951,17 @@
[else (domain-error who y)])))
(set-who! remainder
(let ([f (lambda (x y)
(let ([r (- x (* (quotient x y) y))])
;;; filter out outrageous results
;;; try (remainder 1e194 10.0) without this hack...
(if (if (negative? y) (> r y) (< r y))
r
0.0)))])
(let* ([fmod (cflop2 "(cs)mod")]
[f (lambda (x y)
(cond
[(eqv? x 0) 0]
[else
(let ([r (fmod (real->flonum x) (real->flonum y))])
(if (fl= r 0.0)
;; Always return positive 0.0 --- not sure why,
;; but Racket and other Schemes seem to agree
0.0
r))]))])
(lambda (x y)
(type-case y
[(fixnum?)