diff --git a/eval.c b/eval.c index 024f9ecb..734afe97 100644 --- a/eval.c +++ b/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) \ - else if (t && sexp_complexp(z)) return sexp_complex_normalize(f(ctx, z)); +#define maybe_convert_complex(z, f) \ + 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) \ - sexp name (sexp ctx, sexp self, sexp_sint_t n, sexp z) { \ +#define define_math_op(name, cname, f) \ + 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_log, log, 1, sexp_complex_log) -define_math_op(sexp_sin, sin, 1, sexp_complex_sin) -define_math_op(sexp_cos, cos, 1, sexp_complex_cos) -define_math_op(sexp_tan, tan, 1, sexp_complex_tan) -define_math_op(sexp_asin, asin, 1, sexp_complex_asin) -define_math_op(sexp_acos, acos, 1, sexp_complex_acos) -define_math_op(sexp_atan, atan, 1, sexp_complex_atan) -define_math_op(sexp_round, round, 0, sexp_complex_dummy) -define_math_op(sexp_trunc, trunc, 0, sexp_complex_dummy) -define_math_op(sexp_floor, floor, 0, sexp_complex_dummy) -define_math_op(sexp_ceiling, ceil, 0, sexp_complex_dummy) +define_math_op(sexp_exp, exp, sexp_complex_exp) +define_math_op(sexp_log, log, sexp_complex_log) +define_math_op(sexp_sin, sin, sexp_complex_sin) +define_math_op(sexp_cos, cos, sexp_complex_cos) +define_math_op(sexp_tan, tan, sexp_complex_tan) +define_math_op(sexp_asin, asin, sexp_complex_asin) +define_math_op(sexp_acos, acos, sexp_complex_acos) +define_math_op(sexp_atan, atan, sexp_complex_atan) + +#if SEXP_USE_RATIOS +#define maybe_round_ratio(ctx, q, f) \ + 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); diff --git a/include/chibi/bignum.h b/include/chibi/bignum.h index 1e3d3703..17fe7704 100644 --- a/include/chibi/bignum.h +++ b/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); diff --git a/include/chibi/sexp.h b/include/chibi/sexp.h index 862a572f..bb558cb7 100755 --- a/include/chibi/sexp.h +++ b/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)) \ diff --git a/lib/scheme/division.scm b/lib/scheme/division.scm index 08042817..24a0e065 100644 --- a/lib/scheme/division.scm +++ b/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) - (if (> n 0) (ceiling/ n m) (floor/ 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))) + +;; 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)) diff --git a/opt/bignum.c b/opt/bignum.c index 06c3b31b..23dba7e4 100644 --- a/opt/bignum.c +++ b/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);