Fixes and improvements to number-theoretic division operators
authorMark H Weaver <mhw@netris.org>
Thu, 10 Feb 2011 20:40:57 +0000 (15:40 -0500)
committerAndy Wingo <wingo@pobox.com>
Sat, 12 Feb 2011 12:00:43 +0000 (13:00 +0100)
* libguile/numbers.c (scm_euclidean_quotient, scm_euclidean_divide,
  scm_centered_quotient, scm_centered_divide): Fix bug in inum/inum
  case, where (quotient most-negative-fixnum -1) would not be converted
  to a bignum.

  (scm_euclidean_quotient): Be more anal-retentive about calling
  scm_remember_upto_here_1 after mpz_sgn, (even though mpz_sgn is
  documented as being implemented as a macro and certainly won't
  do any allocation).  It's better to be safe than sorry here.

  (scm_euclidean_quotient, scm_centered_quotient): In the bignum/inum
  case, check if the divisor is 1, since this will allow us to avoid
  allocating a new bignum.

  (scm_euclidean_divide, scm_centered_quotient, scm_centered_divide):
  When computing the intermediate truncated quotient (xx / yy) and
  remainder, use (xx % yy) instead of (xx - qq * yy), on the theory that
  the compiler is more likely to handle this case intelligently and
  maybe combine the operations.

  (scm_euclidean_divide): In the bignum/inum case, we know that the
  remainder will fit in an fixnum, so don't bother allocating a bignum
  for it.

  (scm_euclidean_quotient, scm_euclidean_remainder,
  scm_euclidean_divide, scm_centered_quotient, scm_centered_remainder,
  scm_centered_divide): Minor stylistic changes.

* test-suite/tests/numbers.test: Rework testing framework for
  number-theoretic division operators to be more efficient and
  comprehensive in its testing of code paths and problem cases.

libguile/numbers.c
test-suite/tests/numbers.test

index 1aed0c2..05840ef 100644 (file)
@@ -1089,6 +1089,7 @@ SCM_PRIMITIVE_GENERIC (scm_euclidean_quotient, "euclidean-quotient", 2, 0, 0,
 {
   if (SCM_LIKELY (SCM_I_INUMP (x)))
     {
+      scm_t_inum xx = SCM_I_INUM (x);
       if (SCM_LIKELY (SCM_I_INUMP (y)))
        {
          scm_t_inum yy = SCM_I_INUM (y);
@@ -1096,7 +1097,6 @@ SCM_PRIMITIVE_GENERIC (scm_euclidean_quotient, "euclidean-quotient", 2, 0, 0,
            scm_num_overflow (s_scm_euclidean_quotient);
          else
            {
-             scm_t_inum xx = SCM_I_INUM (x);
              scm_t_inum qq = xx / yy;
              if (xx < qq * yy)
                {
@@ -1105,19 +1105,25 @@ SCM_PRIMITIVE_GENERIC (scm_euclidean_quotient, "euclidean-quotient", 2, 0, 0,
                  else
                    qq++;
                }
-             return SCM_I_MAKINUM (qq);
+             if (SCM_LIKELY (SCM_FIXABLE (qq)))
+               return SCM_I_MAKINUM (qq);
+             else
+               return scm_i_inum2big (qq);
            }
        }
       else if (SCM_BIGP (y))
        {
-         if (SCM_I_INUM (x) >= 0)
+         if (xx >= 0)
            return SCM_INUM0;
          else
-           return SCM_I_MAKINUM (- mpz_sgn (SCM_I_BIG_MPZ (y)));
+           {
+             scm_t_inum qq = - mpz_sgn (SCM_I_BIG_MPZ (y));
+             scm_remember_upto_here_1 (y);
+             return SCM_I_MAKINUM (qq);
+           }
        }
       else if (SCM_REALP (y))
-       return scm_i_inexact_euclidean_quotient
-         (SCM_I_INUM (x), SCM_REAL_VALUE (y));
+       return scm_i_inexact_euclidean_quotient (xx, SCM_REAL_VALUE (y));
       else if (SCM_FRACTIONP (y))
        return scm_i_slow_exact_euclidean_quotient (x, y);
       else
@@ -1131,6 +1137,8 @@ SCM_PRIMITIVE_GENERIC (scm_euclidean_quotient, "euclidean-quotient", 2, 0, 0,
          scm_t_inum yy = SCM_I_INUM (y);
          if (SCM_UNLIKELY (yy == 0))
            scm_num_overflow (s_scm_euclidean_quotient);
+         else if (SCM_UNLIKELY (yy == 1))
+           return x;
          else
            {
              SCM q = scm_i_mkbig ();
@@ -1246,6 +1254,7 @@ SCM_PRIMITIVE_GENERIC (scm_euclidean_remainder, "euclidean-remainder", 2, 0, 0,
 {
   if (SCM_LIKELY (SCM_I_INUMP (x)))
     {
+      scm_t_inum xx = SCM_I_INUM (x);
       if (SCM_LIKELY (SCM_I_INUMP (y)))
        {
          scm_t_inum yy = SCM_I_INUM (y);
@@ -1253,7 +1262,7 @@ SCM_PRIMITIVE_GENERIC (scm_euclidean_remainder, "euclidean-remainder", 2, 0, 0,
            scm_num_overflow (s_scm_euclidean_remainder);
          else
            {
-             scm_t_inum rr = SCM_I_INUM (x) % yy;
+             scm_t_inum rr = xx % yy;
              if (rr >= 0)
                return SCM_I_MAKINUM (rr);
              else if (yy > 0)
@@ -1264,7 +1273,6 @@ SCM_PRIMITIVE_GENERIC (scm_euclidean_remainder, "euclidean-remainder", 2, 0, 0,
        }
       else if (SCM_BIGP (y))
        {
-         scm_t_inum xx = SCM_I_INUM (x);
          if (xx >= 0)
            return x;
          else if (mpz_sgn (SCM_I_BIG_MPZ (y)) > 0)
@@ -1284,8 +1292,7 @@ SCM_PRIMITIVE_GENERIC (scm_euclidean_remainder, "euclidean-remainder", 2, 0, 0,
            }
        }
       else if (SCM_REALP (y))
-       return scm_i_inexact_euclidean_remainder
-         (SCM_I_INUM (x), SCM_REAL_VALUE (y));
+       return scm_i_inexact_euclidean_remainder (xx, SCM_REAL_VALUE (y));
       else if (SCM_FRACTIONP (y))
        return scm_i_slow_exact_euclidean_remainder (x, y);
       else
@@ -1420,6 +1427,7 @@ SCM_PRIMITIVE_GENERIC (scm_euclidean_divide, "euclidean/", 2, 0, 0,
 {
   if (SCM_LIKELY (SCM_I_INUMP (x)))
     {
+      scm_t_inum xx = SCM_I_INUM (x);
       if (SCM_LIKELY (SCM_I_INUMP (y)))
        {
          scm_t_inum yy = SCM_I_INUM (y);
@@ -1427,9 +1435,10 @@ SCM_PRIMITIVE_GENERIC (scm_euclidean_divide, "euclidean/", 2, 0, 0,
            scm_num_overflow (s_scm_euclidean_divide);
          else
            {
-             scm_t_inum xx = SCM_I_INUM (x);
              scm_t_inum qq = xx / yy;
-             scm_t_inum rr = xx - qq * yy;
+             scm_t_inum rr = xx % yy;
+             SCM q;
+
              if (rr < 0)
                {
                  if (yy > 0)
@@ -1437,13 +1446,15 @@ SCM_PRIMITIVE_GENERIC (scm_euclidean_divide, "euclidean/", 2, 0, 0,
                  else
                    { rr -= yy; qq++; }
                }
-             return scm_values (scm_list_2 (SCM_I_MAKINUM (qq),
-                                            SCM_I_MAKINUM (rr)));
+             if (SCM_LIKELY (SCM_FIXABLE (qq)))
+               q = SCM_I_MAKINUM (qq);
+             else
+               q = scm_i_inum2big (qq);
+             return scm_values (scm_list_2 (q, SCM_I_MAKINUM (rr)));
            }
        }
       else if (SCM_BIGP (y))
        {
-         scm_t_inum xx = SCM_I_INUM (x);
          if (xx >= 0)
            return scm_values (scm_list_2 (SCM_INUM0, x));
          else if (mpz_sgn (SCM_I_BIG_MPZ (y)) > 0)
@@ -1464,8 +1475,7 @@ SCM_PRIMITIVE_GENERIC (scm_euclidean_divide, "euclidean/", 2, 0, 0,
            }
        }
       else if (SCM_REALP (y))
-       return scm_i_inexact_euclidean_divide
-         (SCM_I_INUM (x), SCM_REAL_VALUE (y));
+       return scm_i_inexact_euclidean_divide (xx, SCM_REAL_VALUE (y));
       else if (SCM_FRACTIONP (y))
        return scm_i_slow_exact_euclidean_divide (x, y);
       else
@@ -1482,19 +1492,19 @@ SCM_PRIMITIVE_GENERIC (scm_euclidean_divide, "euclidean/", 2, 0, 0,
          else
            {
              SCM q = scm_i_mkbig ();
-             SCM r = scm_i_mkbig ();
+             scm_t_inum rr;
              if (yy > 0)
-               mpz_fdiv_qr_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r),
-                               SCM_I_BIG_MPZ (x), yy);
+               rr = mpz_fdiv_q_ui (SCM_I_BIG_MPZ (q),
+                                   SCM_I_BIG_MPZ (x), yy);
              else
                {
-                 mpz_fdiv_qr_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r),
-                                 SCM_I_BIG_MPZ (x), -yy);
+                 rr = mpz_fdiv_q_ui (SCM_I_BIG_MPZ (q),
+                                     SCM_I_BIG_MPZ (x), -yy);
                  mpz_neg (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (q));
                }
              scm_remember_upto_here_1 (x);
              return scm_values (scm_list_2 (scm_i_normbig (q),
-                                            scm_i_normbig (r)));
+                                            SCM_I_MAKINUM (rr)));
            }
        }
       else if (SCM_BIGP (y))
@@ -1607,6 +1617,7 @@ SCM_PRIMITIVE_GENERIC (scm_centered_quotient, "centered-quotient", 2, 0, 0,
 {
   if (SCM_LIKELY (SCM_I_INUMP (x)))
     {
+      scm_t_inum xx = SCM_I_INUM (x);
       if (SCM_LIKELY (SCM_I_INUMP (y)))
        {
          scm_t_inum yy = SCM_I_INUM (y);
@@ -1614,9 +1625,8 @@ SCM_PRIMITIVE_GENERIC (scm_centered_quotient, "centered-quotient", 2, 0, 0,
            scm_num_overflow (s_scm_centered_quotient);
          else
            {
-             scm_t_inum xx = SCM_I_INUM (x);
              scm_t_inum qq = xx / yy;
-             scm_t_inum rr = xx - qq * yy;
+             scm_t_inum rr = xx % yy;
              if (SCM_LIKELY (xx > 0))
                {
                  if (SCM_LIKELY (yy > 0))
@@ -1643,19 +1653,20 @@ SCM_PRIMITIVE_GENERIC (scm_centered_quotient, "centered-quotient", 2, 0, 0,
                        qq++;
                    }
                }
-             return SCM_I_MAKINUM (qq);
+             if (SCM_LIKELY (SCM_FIXABLE (qq)))
+               return SCM_I_MAKINUM (qq);
+             else
+               return scm_i_inum2big (qq);
            }
        }
       else if (SCM_BIGP (y))
        {
          /* Pass a denormalized bignum version of x (even though it
             can fit in a fixnum) to scm_i_bigint_centered_quotient */
-         return scm_i_bigint_centered_quotient
-           (scm_i_long2big (SCM_I_INUM (x)), y);
+         return scm_i_bigint_centered_quotient (scm_i_long2big (xx), y);
        }
       else if (SCM_REALP (y))
-       return scm_i_inexact_centered_quotient
-         (SCM_I_INUM (x), SCM_REAL_VALUE (y));
+       return scm_i_inexact_centered_quotient (xx, SCM_REAL_VALUE (y));
       else if (SCM_FRACTIONP (y))
        return scm_i_slow_exact_centered_quotient (x, y);
       else
@@ -1669,6 +1680,8 @@ SCM_PRIMITIVE_GENERIC (scm_centered_quotient, "centered-quotient", 2, 0, 0,
          scm_t_inum yy = SCM_I_INUM (y);
          if (SCM_UNLIKELY (yy == 0))
            scm_num_overflow (s_scm_centered_quotient);
+         else if (SCM_UNLIKELY (yy == 1))
+           return x;
          else
            {
              SCM q = scm_i_mkbig ();
@@ -1833,6 +1846,7 @@ SCM_PRIMITIVE_GENERIC (scm_centered_remainder, "centered-remainder", 2, 0, 0,
 {
   if (SCM_LIKELY (SCM_I_INUMP (x)))
     {
+      scm_t_inum xx = SCM_I_INUM (x);
       if (SCM_LIKELY (SCM_I_INUMP (y)))
        {
          scm_t_inum yy = SCM_I_INUM (y);
@@ -1840,7 +1854,6 @@ SCM_PRIMITIVE_GENERIC (scm_centered_remainder, "centered-remainder", 2, 0, 0,
            scm_num_overflow (s_scm_centered_remainder);
          else
            {
-             scm_t_inum xx = SCM_I_INUM (x);
              scm_t_inum rr = xx % yy;
              if (SCM_LIKELY (xx > 0))
                {
@@ -1875,12 +1888,10 @@ SCM_PRIMITIVE_GENERIC (scm_centered_remainder, "centered-remainder", 2, 0, 0,
        {
          /* Pass a denormalized bignum version of x (even though it
             can fit in a fixnum) to scm_i_bigint_centered_remainder */
-         return scm_i_bigint_centered_remainder
-           (scm_i_long2big (SCM_I_INUM (x)), y);
+         return scm_i_bigint_centered_remainder (scm_i_long2big (xx), y);
        }
       else if (SCM_REALP (y))
-       return scm_i_inexact_centered_remainder
-         (SCM_I_INUM (x), SCM_REAL_VALUE (y));
+       return scm_i_inexact_centered_remainder (xx, SCM_REAL_VALUE (y));
       else if (SCM_FRACTIONP (y))
        return scm_i_slow_exact_centered_remainder (x, y);
       else
@@ -2062,6 +2073,7 @@ SCM_PRIMITIVE_GENERIC (scm_centered_divide, "centered/", 2, 0, 0,
 {
   if (SCM_LIKELY (SCM_I_INUMP (x)))
     {
+      scm_t_inum xx = SCM_I_INUM (x);
       if (SCM_LIKELY (SCM_I_INUMP (y)))
        {
          scm_t_inum yy = SCM_I_INUM (y);
@@ -2069,9 +2081,10 @@ SCM_PRIMITIVE_GENERIC (scm_centered_divide, "centered/", 2, 0, 0,
            scm_num_overflow (s_scm_centered_divide);
          else
            {
-             scm_t_inum xx = SCM_I_INUM (x);
              scm_t_inum qq = xx / yy;
-             scm_t_inum rr = xx - qq * yy;
+             scm_t_inum rr = xx % yy;
+             SCM q;
+
              if (SCM_LIKELY (xx > 0))
                {
                  if (SCM_LIKELY (yy > 0))
@@ -2098,20 +2111,21 @@ SCM_PRIMITIVE_GENERIC (scm_centered_divide, "centered/", 2, 0, 0,
                        { qq++; rr -= yy; }
                    }
                }
-             return scm_values (scm_list_2 (SCM_I_MAKINUM (qq),
-                                            SCM_I_MAKINUM (rr)));
+             if (SCM_LIKELY (SCM_FIXABLE (qq)))
+               q = SCM_I_MAKINUM (qq);
+             else
+               q = scm_i_inum2big (qq);
+             return scm_values (scm_list_2 (q, SCM_I_MAKINUM (rr)));
            }
        }
       else if (SCM_BIGP (y))
        {
          /* Pass a denormalized bignum version of x (even though it
             can fit in a fixnum) to scm_i_bigint_centered_divide */
-         return scm_i_bigint_centered_divide
-           (scm_i_long2big (SCM_I_INUM (x)), y);
+         return scm_i_bigint_centered_divide (scm_i_long2big (xx), y);
        }
       else if (SCM_REALP (y))
-       return scm_i_inexact_centered_divide
-         (SCM_I_INUM (x), SCM_REAL_VALUE (y));
+       return scm_i_inexact_centered_divide (xx, SCM_REAL_VALUE (y));
       else if (SCM_FRACTIONP (y))
        return scm_i_slow_exact_centered_divide (x, y);
       else
index 1c4630e..f738189 100644 (file)
   (pass-if "-100i swings back to 45deg down"
     (eqv-loosely? +7.071-7.071i (sqrt -100.0i))))
 
+
 ;;;
 ;;; euclidean/
 ;;; euclidean-quotient
 
 (with-test-prefix "Number-theoretic division"
 
-  ;; Tests that (lo <= x < hi),
+  ;; Tests that (lo <1 x <2 hi),
   ;; but allowing for imprecision
   ;; if x is inexact.
-  (define (test-within-range? lo hi x)
+  (define (test-within-range? lo <1 x <2 hi)
     (if (exact? x)
-        (and (<= lo x) (< x hi))
+        (and (<1 lo x) (<2 x hi))
         (let ((lo (- lo test-epsilon))
               (hi (+ hi test-epsilon)))
           (<= lo x hi))))
 
-  ;; (cartesian-product-map list '(a b) '(1 2))
-  ;; ==> ((a 1) (a 2) (b 1) (b 2))
-  (define (cartesian-product-map f . lsts)
-    (define (cartmap rev-head lsts)
-      (if (null? lsts)
-          (list (apply f (reverse rev-head)))
-          (append-map (lambda (x) (cartmap (cons x rev-head) (cdr lsts)))
-                      (car lsts))))
-    (cartmap '() lsts))
-
-  (define (cartesian-product-for-each f . lsts)
-    (define (cartfor rev-head lsts)
-      (if (null? lsts)
-          (apply f (reverse rev-head))
-          (for-each (lambda (x) (cartfor (cons x rev-head) (cdr lsts)))
-                    (car lsts))))
-    (cartfor '() lsts))
-
-  (define (safe-euclidean-quotient x y)
-    (cond ((not (and (real? x) (real? y))) (throw 'wrong-type-arg))
-          ((zero? y) (throw 'divide-by-zero))
-          ((nan?  y) (nan))
-          ((positive? y) (floor   (/ x y)))
-          ((negative? y) (ceiling (/ x y)))
-          (else (throw 'unknown-problem))))
-
-  (define (safe-euclidean-remainder x y)
-    (let ((q (safe-euclidean-quotient x y)))
-      (- x (* y q))))
-
   (define (valid-euclidean-answer? x y q r)
-    (if (and (finite? x) (finite? y))
-        (and (eq? (exact? q)
-                  (exact? r)
-                  (and (exact? x) (exact? y)))
-             (integer? q)
-             (test-eqv? r (- x (* q y)))
-             (test-within-range? 0 (abs y) r))
-        (and (test-eqv? q (safe-euclidean-quotient  x y))
-             (test-eqv? r (safe-euclidean-remainder x y)))))
-
-  (define (safe-centered-quotient x y)
-    (cond ((not (and (real? x) (real? y))) (throw 'wrong-type-arg))
-          ((zero? y) (throw 'divide-by-zero))
-          ((nan?  y) (nan))
-          ((positive? y) (floor   (+  1/2 (/ x y))))
-          ((negative? y) (ceiling (+ -1/2 (/ x y))))
-          (else (throw 'unknown-problem))))
-
-  (define (safe-centered-remainder x y)
-    (let ((q (safe-centered-quotient x y)))
-      (- x (* y q))))
+    (and (eq? (exact? q)
+              (exact? r)
+              (and (exact? x) (exact? y)))
+         (test-eqv? r (- x (* q y)))
+         (if (and (finite? x) (finite? y))
+             (and (integer? q)
+                  (test-within-range? 0 <= r < (abs y)))
+             (test-eqv? q (/ x y)))))
 
   (define (valid-centered-answer? x y q r)
-    (if (and (finite? x) (finite? y))
-        (and (eq? (exact? q)
-                  (exact? r)
-                  (and (exact? x) (exact? y)))
-             (integer? q)
-             (test-eqv? r (- x (* q y)))
-             (test-within-range? (* -1/2 (abs y))
-                                 (* +1/2 (abs y))
-                                 r))
-        (and (test-eqv? q (safe-centered-quotient  x y))
-             (test-eqv? r (safe-centered-remainder x y)))))
-
-  (define test-numerators
-    (append (cartesian-product-map * '(1 -1)
-                                   '(123 125 127 130 3 5 10
-                                         123.2 125.0 127.2 130.0
-                                         123/7 125/7 127/7 130/7))
-            (cartesian-product-map * '(1 -1)
-                                   '(123 125 127 130 3 5 10)
-                                   (list 1
-                                         (+ 1 most-positive-fixnum)
-                                         (+ 2 most-positive-fixnum)))
-            (list 0 +0.0 -0.0 +inf.0 -inf.0 +nan.0
-                  most-negative-fixnum
-                  (1+ most-positive-fixnum)
-                  (1- most-negative-fixnum))))
-
-  (define test-denominators
-    (list   10  5  10/7  127/2  10.0  63.5
-            -10 -5 -10/7 -127/2 -10.0 -63.5
-            +inf.0 -inf.0 +nan.0 most-negative-fixnum
-            (+ 1 most-positive-fixnum) (+ -1 most-negative-fixnum)
-            (+ 2 most-positive-fixnum) (+ -2 most-negative-fixnum)))
+    (and (eq? (exact? q)
+              (exact? r)
+              (and (exact? x) (exact? y)))
+         (test-eqv? r (- x (* q y)))
+         (if (and (finite? x) (finite? y))
+             (and (integer? q)
+                  (test-within-range?
+                   (* -1/2 (abs y)) <= r < (* +1/2 (abs y))))
+             (test-eqv? q (/ x y)))))
+
+  (define (for lsts f) (apply for-each f lsts))
+
+  (define big (expt 10 (1+ (inexact->exact (ceiling (log10 fixnum-max))))))
+
+  (define (run-division-tests quo+rem quo rem valid-answer?)
+    (define (test n d)
+      (run-test (list n d) #t
+                (lambda ()
+                  (let-values (((q r) (quo+rem n d)))
+                    (and (test-eqv? q (quo n d))
+                         (test-eqv? r (rem n d))
+                         (valid-answer? n d q r))))))
+    (define (test+/- n d)
+      (test n    d )
+      (test n (- d))
+      (cond ((not (zero? n))
+             (test (- n)    d )
+             (test (- n) (- d)))))
+
+    (define (test-for-exception n d exception)
+      (let ((name (list n d)))
+        (pass-if-exception name exception (quo+rem n d))
+        (pass-if-exception name exception (quo n d))
+        (pass-if-exception name exception (rem n d))))
+
+    (run-test "documented?" #t
+              (lambda ()
+                (and (documented? quo+rem)
+                     (documented? quo)
+                     (documented? rem))))
+
+    (with-test-prefix "inum / inum"
+      (with-test-prefix "fixnum-min / -1"
+        (test fixnum-min -1))
+      (for '((1 2 5 10))  ;; denominators
+           (lambda (d)
+             (for '((0 1 2 5 10))  ;; multiples
+                  (lambda (m)
+                    (for '((-2 -1 0 1 2 3 4 5 7 10
+                               12 15 16 19 20))  ;; offsets
+                         (lambda (b)
+                           (test+/- (+ b (* m d))
+                                    d))))))))
+
+    (with-test-prefix "inum / big"
+      (with-test-prefix "fixnum-min / -fixnum-min"
+        (test fixnum-min (- fixnum-min)))
+      (with-test-prefix "fixnum-max / (2*fixnum-max)"
+        (test+/- fixnum-max (* 2 fixnum-max)))
+      (for `((0 1 2 10 ,(1- fixnum-max) ,fixnum-max))
+           (lambda (n)
+             (test    n  (1+ fixnum-max))
+             (test (- n) (1+ fixnum-max))
+             (test    n  (1- fixnum-min))
+             (test (- n) (1- fixnum-min)))))
+
+    (with-test-prefix "big / inum"
+      (with-test-prefix "-fixnum-min / fixnum-min"
+        (test (- fixnum-min) fixnum-min))
+      (for '((1 4 5 10))  ;; denominators
+           (lambda (d)
+             (for `((1 2 5 ,@(if (even? d)
+                                 '(1/2 3/2 5/2)
+                                 '())))  ;; multiples
+                  (lambda (m)
+                    (for '((-2 -1 0 1 2))  ;; offsets
+                         (lambda (b)
+                           (test+/- (+ b (* m d big))
+                                    d))))))))
+
+    (with-test-prefix "big / big"
+      (for `((,big ,(1+ big)))  ;; denominators
+           (lambda (d)
+             (for `((1 2 5 ,@(if (even? d)
+                                 '(1/2 3/2 5/2)
+                                 '())))  ;; multiples
+                  (lambda (m)
+                    (for '((-2 -1 0 1 2))  ;; offsets
+                         (lambda (b)
+                           (test+/- (+ b (* m d))
+                                    d))))))))
+
+    (with-test-prefix "inexact"
+      (for '((0.5 1.5 2.25 5.75))  ;; denominators
+           (lambda (d)
+             (for '((0 1 2 5 1/2 3/2 5/2))  ;; multiples
+                  (lambda (m)
+                    (for '((-2 -1 0 1 2))  ;; offsets
+                         (lambda (b)
+                           (test+/- (+ b (* m d))
+                                    d))))))))
+
+    (with-test-prefix "fractions"
+      (for '((1/10 16/3 10/7))  ;; denominators
+           (lambda (d)
+             (for '((0 1 2 5 1/2 3/2 5/2))  ;; multiples
+                  (lambda (m)
+                    (for '((-2/9 -1/11 0 1/3 2/3))  ;; offsets
+                         (lambda (b)
+                           (test+/- (+ b (* m d))
+                                    d))))))))
+
+    (with-test-prefix "mixed types"
+      (for `((10 ,big 12.0 10/7 +inf.0 -inf.0 +nan.0))  ;; denominators
+           (lambda (d)
+             (for `((25 ,(* 3/2 big) 130.0 15/7
+                        0 0.0 -0.0 +inf.0 -inf.0 +nan.0))  ;; numerators
+                  (lambda (n)
+                    (test+/- n d))))))
+
+    (with-test-prefix "divide by zero"
+      (for `((0 0.0 +0.0))  ;; denominators
+           (lambda (d)
+             (for `((15 ,(* 3/2 big) 18.0 33/7
+                        0 0.0 -0.0 +inf.0 -inf.0 +nan.0))  ;; numerators
+                  (lambda (n)
+                    (test-for-exception
+                     n d exception:numerical-overflow)))))))
 
   (with-test-prefix "euclidean/"
-    (pass-if (documented? euclidean/))
-    (pass-if (documented? euclidean-quotient))
-    (pass-if (documented? euclidean-remainder))
-
-    (cartesian-product-for-each
-     (lambda (n d)
-       (run-test (list 'euclidean/ n d) #t
-                 (lambda ()
-                   (let-values (((q r) (euclidean/ n d)))
-                     (and (test-eqv? q (euclidean-quotient n d))
-                          (test-eqv? r (euclidean-remainder n d))
-                          (valid-euclidean-answer? n d q r))))))
-     test-numerators test-denominators))
+    (run-division-tests euclidean/
+                        euclidean-quotient
+                        euclidean-remainder
+                        valid-euclidean-answer?))
 
   (with-test-prefix "centered/"
-    (pass-if (documented? centered/))
-    (pass-if (documented? centered-quotient))
-    (pass-if (documented? centered-remainder))
-
-    (cartesian-product-for-each
-     (lambda (n d)
-       (run-test (list 'centered/ n d) #t
-                 (lambda ()
-                   (let-values (((q r) (centered/ n d)))
-                     (and (test-eqv? q (centered-quotient n d))
-                          (test-eqv? r (centered-remainder n d))
-                          (valid-centered-answer? n d q r))))))
-     test-numerators test-denominators)))
+    (run-division-tests centered/
+                        centered-quotient
+                        centered-remainder
+                        valid-centered-answer?)))