Rounding functions now support rational arguments exactly.

Fixing the division operators.

eval.c

@@ -1194,15 +1194,15 @@ sexp sexp_register_optimization (sexp ctx, sexp self, sexp_sint_t n, sexp f, sex
 #endif
 
 #if SEXP_USE_COMPLEX
-#define maybe_convert_complex(z, t, f)                                  \
+#define maybe_convert_complex(z, f)                                     \
-  else if (t && sexp_complexp(z)) return sexp_complex_normalize(f(ctx, z));
+  else if (sexp_complexp(z)) return sexp_complex_normalize(f(ctx, z));
 #define sexp_complex_dummy(ctx, z) 0
 #else
-#define maybe_convert_complex(z, t, f)
+#define maybe_convert_complex(z, f)
 #endif
 
-#define define_math_op(name, cname, t, f)                               \
+#define define_math_op(name, cname, f)                                  \
-  sexp name (sexp ctx, sexp self, sexp_sint_t n, sexp z) {               \
+  sexp name (sexp ctx, sexp self, sexp_sint_t n, sexp z) {              \
     double d;                                                           \
     if (sexp_flonump(z))                                                \
       d = sexp_flonum_value(z);                                         \
@@ -1210,24 +1210,49 @@ sexp sexp_register_optimization (sexp ctx, sexp self, sexp_sint_t n, sexp f, sex
       d = (double)sexp_unbox_fixnum(z);                                 \
     maybe_convert_ratio(z)                                              \
     maybe_convert_bignum(z)                                             \
-    maybe_convert_complex(z, t, f)                                      \
+    maybe_convert_complex(z, f)                                         \
     else                                                                \
       return sexp_type_exception(ctx, self, SEXP_NUMBER, z);            \
     return sexp_make_flonum(ctx, cname(d));                             \
   }
 
-define_math_op(sexp_exp, exp, 1, sexp_complex_exp)
+define_math_op(sexp_exp, exp, sexp_complex_exp)
-define_math_op(sexp_log, log, 1, sexp_complex_log)
+define_math_op(sexp_log, log, sexp_complex_log)
-define_math_op(sexp_sin, sin, 1, sexp_complex_sin)
+define_math_op(sexp_sin, sin, sexp_complex_sin)
-define_math_op(sexp_cos, cos, 1, sexp_complex_cos)
+define_math_op(sexp_cos, cos, sexp_complex_cos)
-define_math_op(sexp_tan, tan, 1, sexp_complex_tan)
+define_math_op(sexp_tan, tan, sexp_complex_tan)
-define_math_op(sexp_asin, asin, 1, sexp_complex_asin)
+define_math_op(sexp_asin, asin, sexp_complex_asin)
-define_math_op(sexp_acos, acos, 1, sexp_complex_acos)
+define_math_op(sexp_acos, acos, sexp_complex_acos)
-define_math_op(sexp_atan, atan, 1, sexp_complex_atan)
+define_math_op(sexp_atan, atan, sexp_complex_atan)
-define_math_op(sexp_round, round, 0, sexp_complex_dummy)
+
-define_math_op(sexp_trunc, trunc, 0, sexp_complex_dummy)
+#if SEXP_USE_RATIOS
-define_math_op(sexp_floor, floor, 0, sexp_complex_dummy)
+#define maybe_round_ratio(ctx, q, f)            \
-define_math_op(sexp_ceiling, ceil, 0, sexp_complex_dummy)
+  if (sexp_ratiop(q)) return f(ctx, q);
+#else
+#define maybe_round_ratio(ctx, q, f)
+#endif
+
+#define define_math_rounder(name, cname, f)                             \
+  sexp name (sexp ctx, sexp self, sexp_sint_t n, sexp z) {              \
+    maybe_round_ratio(ctx, z, f)                                        \
+    if (sexp_flonump(z))                                                \
+      return sexp_make_flonum(ctx, cname(sexp_flonum_value(z)));        \
+    else if (sexp_fixnump(z) || sexp_bignump(z))                        \
+      return z;                                                         \
+    return sexp_type_exception(ctx, self, SEXP_NUMBER, z);              \
+  }
+
+static double even_round (double d) {
+  double res = round(d);
+  if (fabs(d - res) == 0.5 && ((long)res & 1))
+    res += (res < 0) ? 1 : -1;
+  return res;
+}
+
+define_math_rounder(sexp_round, even_round, sexp_ratio_round)
+define_math_rounder(sexp_trunc, trunc, sexp_ratio_trunc)
+define_math_rounder(sexp_floor, floor, sexp_ratio_floor)
+define_math_rounder(sexp_ceiling, ceil, sexp_ratio_ceiling)
 
 sexp sexp_sqrt (sexp ctx, sexp self, sexp_sint_t n, sexp z) {
 #if SEXP_USE_COMPLEX
@@ -1241,7 +1266,7 @@ sexp sexp_sqrt (sexp ctx, sexp self, sexp_sint_t n, sexp z) {
     d = (double)sexp_unbox_fixnum(z);
   maybe_convert_bignum(z)       /* XXXX add bignum sqrt */
   maybe_convert_ratio(z)        /* XXXX add ratio sqrt */
-  maybe_convert_complex(z, 1, sexp_complex_sqrt)
+  maybe_convert_complex(z, sexp_complex_sqrt)
   else
     return sexp_type_exception(ctx, self, SEXP_NUMBER, z);
   sexp_gc_preserve1(ctx, res);

include/chibi/bignum.h

@@ -42,6 +42,10 @@ sexp sexp_remainder (sexp ctx, sexp a, sexp b);
 double sexp_ratio_to_double (sexp rat);
 sexp sexp_make_ratio (sexp ctx, sexp num, sexp den);
 sexp sexp_ratio_normalize (sexp ctx, sexp rat, sexp in);
+sexp sexp_ratio_round (sexp ctx, sexp a);
+sexp sexp_ratio_trunc (sexp ctx, sexp a);
+sexp sexp_ratio_floor (sexp ctx, sexp a);
+sexp sexp_ratio_ceiling (sexp ctx, sexp a);
 #endif
 #if SEXP_USE_COMPLEX
 sexp sexp_make_complex (sexp ctx, sexp real, sexp image);

include/chibi/sexp.h

@@ -698,6 +698,15 @@ SEXP_API sexp sexp_make_unsigned_integer(sexp ctx, sexp_luint_t x);
                                  && (sexp_bignum_sign(x) < 0))
 #define sexp_negativep(x) (sexp_exact_negativep(x) ||                   \
                            (sexp_flonump(x) && sexp_flonum_value(x) < 0))
+#define sexp_positivep(x) (!(sexp_negativep(x)))
+
+#if SEXP_USE_BIGNUMS
+#define sexp_oddp(x) (sexp_fixnump(x) ? sexp_unbox_fixnum(x) & 1 : \
+                      sexp_bignump(x) && (sexp_bignum_data(x)[0] & 1))
+#else
+#define sexp_oddp(x) (sexp_fixnump(x) && (sexp_unbox_fixnum(x) & 1))
+#endif
+#define sexp_evenp(x) (!(sexp_oddp(x)))
 
 #define sexp_negate_exact(x)                            \
   if (sexp_bignump(x))                                  \

lib/scheme/division.scm

@@ -1,22 +1,23 @@
 
+;; The builtin quotient and remainder implement truncation - the
+;; fractional part is always discarded.
+
 (define truncate-quotient quotient)
 (define truncate-remainder remainder)
 (define (truncate/ n m)
   (values (truncate-quotient n m) (truncate-remainder n m)))
 
-(define floor-remainder modulo)
+;; Floor, ceiling and round just compose their corresponding function
+;; with division to determine the quotient, and compute the remainder
+;; from that.
+
 (define (floor-quotient n m)
-  (quotient (- n (floor-remainder n m)) m))
+  (inexact->exact (floor (/ n m))))
+(define (floor-remainder n m)
+  (- n (* m (floor-quotient n m))))
 (define (floor/ n m)
   (values (floor-quotient n m) (floor-remainder n m)))
 
-(define (round-quotient n m)
-  (inexact->exact (round (/ n m))))
-(define (round-remainder n m)
-  (- n (* m (round-quotient n m))))
-(define (round/ n m)
-  (values (round-quotient n m) (round-remainder n m)))
-
 (define (ceiling-quotient n m)
   (inexact->exact (ceiling (/ n m))))
 (define (ceiling-remainder n m)
@@ -24,12 +25,25 @@
 (define (ceiling/ n m)
   (values (ceiling-quotient n m) (ceiling-remainder n m)))
 
-(define (euclidean/ n m)
+(define (round-quotient n m)
-  (if (> n 0) (ceiling/ n m) (floor/ n m)))
+  (inexact->exact (round (/ n m))))
+(define (round-remainder n m)
+  (- n (* m (round-quotient n m))))
+(define (round/ n m)
+  (values (round-quotient n m) (round-remainder n m)))
+
+;; Euclidean is defined as floor if the divisor is negative, and
+;; ceiling otherwise.
+
 (define (euclidean-quotient n m)
-  (if (> n 0) (ceiling-quotient n m) (floor-quotient n m)))
+  (if (> m 0) (floor-quotient n m) (ceiling-quotient n m)))
 (define (euclidean-remainder n m)
-  (if (> n 0) (ceiling-remainder n m) (floor-remainder n m)))
+  (- n (* m (euclidean-quotient n m))))
+(define (euclidean/ n m)
+  (values (euclidean-quotient n m) (euclidean-remainder n m)))
+
+;; Centered places the remainder in the half-open interval
+;; [-m/2, m/2).
 
 (define (centered-remainder n m)
   (let ((r (euclidean-remainder n m))

opt/bignum.c

@@ -538,6 +538,47 @@ sexp sexp_ratio_compare (sexp ctx, sexp a, sexp b) {
   return a2;
 }
 
+sexp sexp_ratio_round (sexp ctx, sexp a) {
+  sexp_gc_var2(q, r);
+  sexp_gc_preserve2(ctx, q, r);
+  q = sexp_quotient(ctx, sexp_ratio_numerator(a), sexp_ratio_denominator(a));
+  if ((sexp_ratio_denominator(a) == SEXP_TWO) && sexp_oddp(q)) {
+    q = sexp_add(ctx, q, (sexp_positivep(q) ? SEXP_ONE : SEXP_NEG_ONE));
+  } else {
+    r = sexp_remainder(ctx, sexp_ratio_numerator(a), sexp_ratio_denominator(a));
+    r = sexp_mul(ctx, r, SEXP_TWO);
+    if (sexp_negativep(r)) {sexp_negate(r);}
+    if (sexp_unbox_fixnum(sexp_compare(ctx, r, sexp_ratio_denominator(a))) > 0)
+      q = sexp_add(ctx, q, (sexp_positivep(q) ? SEXP_ONE : SEXP_NEG_ONE));
+  }
+  sexp_gc_release2(ctx);
+  return q;
+}
+
+sexp sexp_ratio_trunc (sexp ctx, sexp a) {
+  return sexp_quotient(ctx, sexp_ratio_numerator(a), sexp_ratio_denominator(a));
+}
+
+sexp sexp_ratio_floor (sexp ctx, sexp a) {
+  sexp_gc_var1(q);
+  sexp_gc_preserve1(ctx, q);
+  q = sexp_quotient(ctx, sexp_ratio_numerator(a), sexp_ratio_denominator(a));
+  if (sexp_negativep(sexp_ratio_numerator(a)))
+    q = sexp_add(ctx, q, SEXP_NEG_ONE);
+  sexp_gc_release1(ctx);
+  return q;
+}
+
+sexp sexp_ratio_ceiling (sexp ctx, sexp a) {
+  sexp_gc_var1(q);
+  sexp_gc_preserve1(ctx, q);
+  q = sexp_quotient(ctx, sexp_ratio_numerator(a), sexp_ratio_denominator(a));
+  if (sexp_positivep(sexp_ratio_numerator(a)))
+    q = sexp_add(ctx, q, SEXP_ONE);
+  sexp_gc_release1(ctx);
+  return q;
+}
+
 #endif
 
 /************************ complex numbers ****************************/
@@ -1401,7 +1442,7 @@ sexp sexp_compare (sexp ctx, sexp a, sexp b) {
       r = sexp_make_fixnum(f > 0.0 ? 1 : f == 0.0 ? 0 : -1);
       break;
     case SEXP_NUM_FIX_BIG:
-      r = sexp_make_fixnum(-1);
+      r = sexp_make_fixnum(sexp_bignum_sign(b) < 0 ? 1 : -1);
       break;
     case SEXP_NUM_FLO_FLO:
       f = sexp_flonum_value(a) - sexp_flonum_value(b);