X-Git-Url: http://git.hcoop.net/bpt/guile.git/blobdiff_plain/9d427b2cc3a7062403d75074760a2b0feb01cb21..78f0ef20a77dd745a8fa4c1ac1fae4d42b481adf:/libguile/numbers.c diff --git a/libguile/numbers.c b/libguile/numbers.c index 85ca0fdb6..24ae2bc31 100644 --- a/libguile/numbers.c +++ b/libguile/numbers.c @@ -45,6 +45,8 @@ # include #endif +#include + #include #include #include @@ -72,6 +74,9 @@ #ifndef M_LOG10E #define M_LOG10E 0.43429448190325182765 #endif +#ifndef M_LN2 +#define M_LN2 0.69314718055994530942 +#endif #ifndef M_PI #define M_PI 3.14159265358979323846 #endif @@ -111,6 +116,7 @@ typedef scm_t_signed_bits scm_t_inum; static SCM flo0; static SCM exactly_one_half; +static SCM flo_log10e; #define SCM_SWAP(x, y) do { SCM __t = x; x = y; y = __t; } while (0) @@ -130,9 +136,9 @@ static double acosh (double x) { return log (x + sqrt (x * x - 1)); } static double atanh (double x) { return 0.5 * log ((1 + x) / (1 - x)); } #endif -/* mpz_cmp_d in gmp 4.1.3 doesn't recognise infinities, so xmpz_cmp_d uses - an explicit check. In some future gmp (don't know what version number), - mpz_cmp_d is supposed to do this itself. */ +/* mpz_cmp_d in GMP before 4.2 didn't recognise infinities, so + xmpz_cmp_d uses an explicit check. Starting with GMP 4.2 (released + in March 2006), mpz_cmp_d now handles infinities properly. */ #if 1 #define xmpz_cmp_d(z, d) \ (isinf (d) ? (d < 0.0 ? 1 : -1) : mpz_cmp_d (z, d)) @@ -142,7 +148,7 @@ static double atanh (double x) { return 0.5 * log ((1 + x) / (1 - x)); } #if defined (GUILE_I) -#if HAVE_COMPLEX_DOUBLE +#if defined HAVE_COMPLEX_DOUBLE /* For an SCM object Z which is a complex number (ie. satisfies SCM_COMPLEXP), return its value as a C level "complex double". */ @@ -312,16 +318,15 @@ scm_i_dbl2num (double u) we need to use mpz_getlimbn. mpz_tstbit is not right, it treats negatives as twos complement. - In current gmp 4.1.3, mpz_get_d rounding is unspecified. It ends up - following the hardware rounding mode, but applied to the absolute value - of the mpz_t operand. This is not what we want so we put the high - DBL_MANT_DIG bits into a temporary. In some future gmp, don't know when, - mpz_get_d is supposed to always truncate towards zero. + In GMP before 4.2, mpz_get_d rounding was unspecified. It ended up + following the hardware rounding mode, but applied to the absolute + value of the mpz_t operand. This is not what we want so we put the + high DBL_MANT_DIG bits into a temporary. Starting with GMP 4.2 + (released in March 2006) mpz_get_d now always truncates towards zero. - ENHANCE-ME: The temporary init+clear to force the rounding in gmp 4.1.3 - is a slowdown. It'd be faster to pick out the relevant high bits with - mpz_getlimbn if we could be bothered coding that, and if the new - truncating gmp doesn't come out. */ + ENHANCE-ME: The temporary init+clear to force the rounding in GMP + before 4.2 is a slowdown. It'd be faster to pick out the relevant + high bits with mpz_getlimbn. */ double scm_i_big2dbl (SCM b) @@ -333,7 +338,12 @@ scm_i_big2dbl (SCM b) #if 1 { - /* Current GMP, eg. 4.1.3, force truncation towards zero */ + /* For GMP earlier than 4.2, force truncation towards zero */ + + /* FIXME: DBL_MANT_DIG is the number of base-`FLT_RADIX' digits, + _not_ the number of bits, so this code will break badly on a + system with non-binary doubles. */ + mpz_t tmp; if (bits > DBL_MANT_DIG) { @@ -349,7 +359,7 @@ scm_i_big2dbl (SCM b) } } #else - /* Future GMP */ + /* GMP 4.2 or later */ result = mpz_get_d (SCM_I_BIG_MPZ (b)); #endif @@ -522,7 +532,7 @@ SCM_PRIMITIVE_GENERIC (scm_exact_p, "exact?", 1, 0, 0, else if (SCM_NUMBERP (x)) return SCM_BOOL_T; else - SCM_WTA_DISPATCH_1 (g_scm_exact_p, x, 1, s_scm_exact_p); + return scm_wta_dispatch_1 (g_scm_exact_p, x, 1, s_scm_exact_p); } #undef FUNC_NAME @@ -538,7 +548,7 @@ SCM_PRIMITIVE_GENERIC (scm_inexact_p, "inexact?", 1, 0, 0, else if (SCM_NUMBERP (x)) return SCM_BOOL_F; else - SCM_WTA_DISPATCH_1 (g_scm_inexact_p, x, 1, s_scm_inexact_p); + return scm_wta_dispatch_1 (g_scm_inexact_p, x, 1, s_scm_inexact_p); } #undef FUNC_NAME @@ -572,7 +582,7 @@ SCM_PRIMITIVE_GENERIC (scm_odd_p, "odd?", 1, 0, 0, return SCM_BOOL_F; } } - SCM_WTA_DISPATCH_1 (g_scm_odd_p, n, 1, s_scm_odd_p); + return scm_wta_dispatch_1 (g_scm_odd_p, n, 1, s_scm_odd_p); } #undef FUNC_NAME @@ -606,7 +616,7 @@ SCM_PRIMITIVE_GENERIC (scm_even_p, "even?", 1, 0, 0, return SCM_BOOL_T; } } - SCM_WTA_DISPATCH_1 (g_scm_even_p, n, 1, s_scm_even_p); + return scm_wta_dispatch_1 (g_scm_even_p, n, 1, s_scm_even_p); } #undef FUNC_NAME @@ -621,7 +631,7 @@ SCM_PRIMITIVE_GENERIC (scm_finite_p, "finite?", 1, 0, 0, else if (scm_is_real (x)) return SCM_BOOL_T; else - SCM_WTA_DISPATCH_1 (g_scm_finite_p, x, 1, s_scm_finite_p); + return scm_wta_dispatch_1 (g_scm_finite_p, x, 1, s_scm_finite_p); } #undef FUNC_NAME @@ -636,7 +646,7 @@ SCM_PRIMITIVE_GENERIC (scm_inf_p, "inf?", 1, 0, 0, else if (scm_is_real (x)) return SCM_BOOL_F; else - SCM_WTA_DISPATCH_1 (g_scm_inf_p, x, 1, s_scm_inf_p); + return scm_wta_dispatch_1 (g_scm_inf_p, x, 1, s_scm_inf_p); } #undef FUNC_NAME @@ -651,7 +661,7 @@ SCM_PRIMITIVE_GENERIC (scm_nan_p, "nan?", 1, 0, 0, else if (scm_is_real (x)) return SCM_BOOL_F; else - SCM_WTA_DISPATCH_1 (g_scm_nan_p, x, 1, s_scm_nan_p); + return scm_wta_dispatch_1 (g_scm_nan_p, x, 1, s_scm_nan_p); } #undef FUNC_NAME @@ -778,7 +788,7 @@ SCM_PRIMITIVE_GENERIC (scm_abs, "abs", 1, 0, 0, SCM_FRACTION_DENOMINATOR (x)); } else - SCM_WTA_DISPATCH_1 (g_scm_abs, x, 1, s_scm_abs); + return scm_wta_dispatch_1 (g_scm_abs, x, 1, s_scm_abs); } #undef FUNC_NAME @@ -788,78 +798,15 @@ SCM_PRIMITIVE_GENERIC (scm_quotient, "quotient", 2, 0, 0, "Return the quotient of the numbers @var{x} and @var{y}.") #define FUNC_NAME s_scm_quotient { - 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); - if (SCM_UNLIKELY (yy == 0)) - scm_num_overflow (s_scm_quotient); - else - { - scm_t_inum z = xx / yy; - if (SCM_LIKELY (SCM_FIXABLE (z))) - return SCM_I_MAKINUM (z); - else - return scm_i_inum2big (z); - } - } - else if (SCM_BIGP (y)) - { - if ((SCM_I_INUM (x) == SCM_MOST_NEGATIVE_FIXNUM) - && (mpz_cmp_ui (SCM_I_BIG_MPZ (y), - - SCM_MOST_NEGATIVE_FIXNUM) == 0)) - { - /* Special case: x == fixnum-min && y == abs (fixnum-min) */ - scm_remember_upto_here_1 (y); - return SCM_I_MAKINUM (-1); - } - else - return SCM_INUM0; - } - else - SCM_WTA_DISPATCH_2 (g_scm_quotient, x, y, SCM_ARG2, s_scm_quotient); - } - else if (SCM_BIGP (x)) + if (SCM_LIKELY (scm_is_integer (x))) { - if (SCM_LIKELY (SCM_I_INUMP (y))) - { - scm_t_inum yy = SCM_I_INUM (y); - if (SCM_UNLIKELY (yy == 0)) - scm_num_overflow (s_scm_quotient); - else if (SCM_UNLIKELY (yy == 1)) - return x; - else - { - SCM result = scm_i_mkbig (); - if (yy < 0) - { - mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result), - SCM_I_BIG_MPZ (x), - - yy); - mpz_neg (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (result)); - } - else - mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (x), yy); - scm_remember_upto_here_1 (x); - return scm_i_normbig (result); - } - } - else if (SCM_BIGP (y)) - { - SCM result = scm_i_mkbig (); - mpz_tdiv_q (SCM_I_BIG_MPZ (result), - SCM_I_BIG_MPZ (x), - SCM_I_BIG_MPZ (y)); - scm_remember_upto_here_2 (x, y); - return scm_i_normbig (result); - } + if (SCM_LIKELY (scm_is_integer (y))) + return scm_truncate_quotient (x, y); else - SCM_WTA_DISPATCH_2 (g_scm_quotient, x, y, SCM_ARG2, s_scm_quotient); + return scm_wta_dispatch_2 (g_scm_quotient, x, y, SCM_ARG2, s_scm_quotient); } else - SCM_WTA_DISPATCH_2 (g_scm_quotient, x, y, SCM_ARG1, s_scm_quotient); + return scm_wta_dispatch_2 (g_scm_quotient, x, y, SCM_ARG1, s_scm_quotient); } #undef FUNC_NAME @@ -872,69 +819,15 @@ SCM_PRIMITIVE_GENERIC (scm_remainder, "remainder", 2, 0, 0, "@end lisp") #define FUNC_NAME s_scm_remainder { - if (SCM_LIKELY (SCM_I_INUMP (x))) - { - if (SCM_LIKELY (SCM_I_INUMP (y))) - { - scm_t_inum yy = SCM_I_INUM (y); - if (SCM_UNLIKELY (yy == 0)) - scm_num_overflow (s_scm_remainder); - else - { - /* C99 specifies that "%" is the remainder corresponding to a - quotient rounded towards zero, and that's also traditional - for machine division, so z here should be well defined. */ - scm_t_inum z = SCM_I_INUM (x) % yy; - return SCM_I_MAKINUM (z); - } - } - else if (SCM_BIGP (y)) - { - if ((SCM_I_INUM (x) == SCM_MOST_NEGATIVE_FIXNUM) - && (mpz_cmp_ui (SCM_I_BIG_MPZ (y), - - SCM_MOST_NEGATIVE_FIXNUM) == 0)) - { - /* Special case: x == fixnum-min && y == abs (fixnum-min) */ - scm_remember_upto_here_1 (y); - return SCM_INUM0; - } - else - return x; - } - else - SCM_WTA_DISPATCH_2 (g_scm_remainder, x, y, SCM_ARG2, s_scm_remainder); - } - else if (SCM_BIGP (x)) + if (SCM_LIKELY (scm_is_integer (x))) { - if (SCM_LIKELY (SCM_I_INUMP (y))) - { - scm_t_inum yy = SCM_I_INUM (y); - if (SCM_UNLIKELY (yy == 0)) - scm_num_overflow (s_scm_remainder); - else - { - SCM result = scm_i_mkbig (); - if (yy < 0) - yy = - yy; - mpz_tdiv_r_ui (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ(x), yy); - scm_remember_upto_here_1 (x); - return scm_i_normbig (result); - } - } - else if (SCM_BIGP (y)) - { - SCM result = scm_i_mkbig (); - mpz_tdiv_r (SCM_I_BIG_MPZ (result), - SCM_I_BIG_MPZ (x), - SCM_I_BIG_MPZ (y)); - scm_remember_upto_here_2 (x, y); - return scm_i_normbig (result); - } + if (SCM_LIKELY (scm_is_integer (y))) + return scm_truncate_remainder (x, y); else - SCM_WTA_DISPATCH_2 (g_scm_remainder, x, y, SCM_ARG2, s_scm_remainder); + return scm_wta_dispatch_2 (g_scm_remainder, x, y, SCM_ARG2, s_scm_remainder); } else - SCM_WTA_DISPATCH_2 (g_scm_remainder, x, y, SCM_ARG1, s_scm_remainder); + return scm_wta_dispatch_2 (g_scm_remainder, x, y, SCM_ARG1, s_scm_remainder); } #undef FUNC_NAME @@ -948,181 +841,171 @@ SCM_PRIMITIVE_GENERIC (scm_modulo, "modulo", 2, 0, 0, "@end lisp") #define FUNC_NAME s_scm_modulo { - 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); - if (SCM_UNLIKELY (yy == 0)) - scm_num_overflow (s_scm_modulo); - else - { - /* C99 specifies that "%" is the remainder corresponding to a - quotient rounded towards zero, and that's also traditional - for machine division, so z here should be well defined. */ - scm_t_inum z = xx % yy; - scm_t_inum result; - - if (yy < 0) - { - if (z > 0) - result = z + yy; - else - result = z; - } - else - { - if (z < 0) - result = z + yy; - else - result = z; - } - return SCM_I_MAKINUM (result); - } - } - else if (SCM_BIGP (y)) - { - int sgn_y = mpz_sgn (SCM_I_BIG_MPZ (y)); - { - mpz_t z_x; - SCM result; - - if (sgn_y < 0) - { - SCM pos_y = scm_i_clonebig (y, 0); - /* do this after the last scm_op */ - mpz_init_set_si (z_x, xx); - result = pos_y; /* re-use this bignum */ - mpz_mod (SCM_I_BIG_MPZ (result), - z_x, - SCM_I_BIG_MPZ (pos_y)); - scm_remember_upto_here_1 (pos_y); - } - else - { - result = scm_i_mkbig (); - /* do this after the last scm_op */ - mpz_init_set_si (z_x, xx); - mpz_mod (SCM_I_BIG_MPZ (result), - z_x, - SCM_I_BIG_MPZ (y)); - scm_remember_upto_here_1 (y); - } - - if ((sgn_y < 0) && mpz_sgn (SCM_I_BIG_MPZ (result)) != 0) - mpz_add (SCM_I_BIG_MPZ (result), - SCM_I_BIG_MPZ (y), - SCM_I_BIG_MPZ (result)); - scm_remember_upto_here_1 (y); - /* and do this before the next one */ - mpz_clear (z_x); - return scm_i_normbig (result); - } - } - else - SCM_WTA_DISPATCH_2 (g_scm_modulo, x, y, SCM_ARG2, s_scm_modulo); - } - else if (SCM_BIGP (x)) + if (SCM_LIKELY (scm_is_integer (x))) { - if (SCM_LIKELY (SCM_I_INUMP (y))) - { - scm_t_inum yy = SCM_I_INUM (y); - if (SCM_UNLIKELY (yy == 0)) - scm_num_overflow (s_scm_modulo); - else - { - SCM result = scm_i_mkbig (); - mpz_mod_ui (SCM_I_BIG_MPZ (result), - SCM_I_BIG_MPZ (x), - (yy < 0) ? - yy : yy); - scm_remember_upto_here_1 (x); - if ((yy < 0) && (mpz_sgn (SCM_I_BIG_MPZ (result)) != 0)) - mpz_sub_ui (SCM_I_BIG_MPZ (result), - SCM_I_BIG_MPZ (result), - - yy); - return scm_i_normbig (result); - } - } - else if (SCM_BIGP (y)) - { - SCM result = scm_i_mkbig (); - int y_sgn = mpz_sgn (SCM_I_BIG_MPZ (y)); - SCM pos_y = scm_i_clonebig (y, y_sgn >= 0); - mpz_mod (SCM_I_BIG_MPZ (result), - SCM_I_BIG_MPZ (x), - SCM_I_BIG_MPZ (pos_y)); - - scm_remember_upto_here_1 (x); - if ((y_sgn < 0) && (mpz_sgn (SCM_I_BIG_MPZ (result)) != 0)) - mpz_add (SCM_I_BIG_MPZ (result), - SCM_I_BIG_MPZ (y), - SCM_I_BIG_MPZ (result)); - scm_remember_upto_here_2 (y, pos_y); - return scm_i_normbig (result); - } + if (SCM_LIKELY (scm_is_integer (y))) + return scm_floor_remainder (x, y); else - SCM_WTA_DISPATCH_2 (g_scm_modulo, x, y, SCM_ARG2, s_scm_modulo); + return scm_wta_dispatch_2 (g_scm_modulo, x, y, SCM_ARG2, s_scm_modulo); } else - SCM_WTA_DISPATCH_2 (g_scm_modulo, x, y, SCM_ARG1, s_scm_modulo); + return scm_wta_dispatch_2 (g_scm_modulo, x, y, SCM_ARG1, s_scm_modulo); +} +#undef FUNC_NAME + +/* two_valued_wta_dispatch_2 is a version of SCM_WTA_DISPATCH_2 for + two-valued functions. It is called from primitive generics that take + two arguments and return two values, when the core procedure is + unable to handle the given argument types. If there are GOOPS + methods for this primitive generic, it dispatches to GOOPS and, if + successful, expects two values to be returned, which are placed in + *rp1 and *rp2. If there are no GOOPS methods, it throws a + wrong-type-arg exception. + + FIXME: This obviously belongs somewhere else, but until we decide on + the right API, it is here as a static function, because it is needed + by the *_divide functions below. +*/ +static void +two_valued_wta_dispatch_2 (SCM gf, SCM a1, SCM a2, int pos, + const char *subr, SCM *rp1, SCM *rp2) +{ + SCM vals = scm_wta_dispatch_2 (gf, a1, a2, pos, subr); + + scm_i_extract_values_2 (vals, rp1, rp2); +} + +SCM_DEFINE (scm_euclidean_quotient, "euclidean-quotient", 2, 0, 0, + (SCM x, SCM y), + "Return the integer @var{q} such that\n" + "@math{@var{x} = @var{q}*@var{y} + @var{r}}\n" + "where @math{0 <= @var{r} < abs(@var{y})}.\n" + "@lisp\n" + "(euclidean-quotient 123 10) @result{} 12\n" + "(euclidean-quotient 123 -10) @result{} -12\n" + "(euclidean-quotient -123 10) @result{} -13\n" + "(euclidean-quotient -123 -10) @result{} 13\n" + "(euclidean-quotient -123.2 -63.5) @result{} 2.0\n" + "(euclidean-quotient 16/3 -10/7) @result{} -3\n" + "@end lisp") +#define FUNC_NAME s_scm_euclidean_quotient +{ + if (scm_is_false (scm_negative_p (y))) + return scm_floor_quotient (x, y); + else + return scm_ceiling_quotient (x, y); +} +#undef FUNC_NAME + +SCM_DEFINE (scm_euclidean_remainder, "euclidean-remainder", 2, 0, 0, + (SCM x, SCM y), + "Return the real number @var{r} such that\n" + "@math{0 <= @var{r} < abs(@var{y})} and\n" + "@math{@var{x} = @var{q}*@var{y} + @var{r}}\n" + "for some integer @var{q}.\n" + "@lisp\n" + "(euclidean-remainder 123 10) @result{} 3\n" + "(euclidean-remainder 123 -10) @result{} 3\n" + "(euclidean-remainder -123 10) @result{} 7\n" + "(euclidean-remainder -123 -10) @result{} 7\n" + "(euclidean-remainder -123.2 -63.5) @result{} 3.8\n" + "(euclidean-remainder 16/3 -10/7) @result{} 22/21\n" + "@end lisp") +#define FUNC_NAME s_scm_euclidean_remainder +{ + if (scm_is_false (scm_negative_p (y))) + return scm_floor_remainder (x, y); + else + return scm_ceiling_remainder (x, y); +} +#undef FUNC_NAME + +SCM_DEFINE (scm_i_euclidean_divide, "euclidean/", 2, 0, 0, + (SCM x, SCM y), + "Return the integer @var{q} and the real number @var{r}\n" + "such that @math{@var{x} = @var{q}*@var{y} + @var{r}}\n" + "and @math{0 <= @var{r} < abs(@var{y})}.\n" + "@lisp\n" + "(euclidean/ 123 10) @result{} 12 and 3\n" + "(euclidean/ 123 -10) @result{} -12 and 3\n" + "(euclidean/ -123 10) @result{} -13 and 7\n" + "(euclidean/ -123 -10) @result{} 13 and 7\n" + "(euclidean/ -123.2 -63.5) @result{} 2.0 and 3.8\n" + "(euclidean/ 16/3 -10/7) @result{} -3 and 22/21\n" + "@end lisp") +#define FUNC_NAME s_scm_i_euclidean_divide +{ + if (scm_is_false (scm_negative_p (y))) + return scm_i_floor_divide (x, y); + else + return scm_i_ceiling_divide (x, y); } #undef FUNC_NAME -static SCM scm_i_inexact_euclidean_quotient (double x, double y); -static SCM scm_i_slow_exact_euclidean_quotient (SCM x, SCM y); +void +scm_euclidean_divide (SCM x, SCM y, SCM *qp, SCM *rp) +{ + if (scm_is_false (scm_negative_p (y))) + return scm_floor_divide (x, y, qp, rp); + else + return scm_ceiling_divide (x, y, qp, rp); +} + +static SCM scm_i_inexact_floor_quotient (double x, double y); +static SCM scm_i_exact_rational_floor_quotient (SCM x, SCM y); -SCM_PRIMITIVE_GENERIC (scm_euclidean_quotient, "euclidean-quotient", 2, 0, 0, +SCM_PRIMITIVE_GENERIC (scm_floor_quotient, "floor-quotient", 2, 0, 0, (SCM x, SCM y), - "Return the integer @var{q} such that\n" - "@math{@var{x} = @var{q}*@var{y} + @var{r}}\n" - "where @math{0 <= @var{r} < abs(@var{y})}.\n" + "Return the floor of @math{@var{x} / @var{y}}.\n" "@lisp\n" - "(euclidean-quotient 123 10) @result{} 12\n" - "(euclidean-quotient 123 -10) @result{} -12\n" - "(euclidean-quotient -123 10) @result{} -13\n" - "(euclidean-quotient -123 -10) @result{} 13\n" - "(euclidean-quotient -123.2 -63.5) @result{} 2.0\n" - "(euclidean-quotient 16/3 -10/7) @result{} -3\n" + "(floor-quotient 123 10) @result{} 12\n" + "(floor-quotient 123 -10) @result{} -13\n" + "(floor-quotient -123 10) @result{} -13\n" + "(floor-quotient -123 -10) @result{} 12\n" + "(floor-quotient -123.2 -63.5) @result{} 1.0\n" + "(floor-quotient 16/3 -10/7) @result{} -4\n" "@end lisp") -#define FUNC_NAME s_scm_euclidean_quotient +#define FUNC_NAME s_scm_floor_quotient { 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); - if (SCM_UNLIKELY (yy == 0)) - scm_num_overflow (s_scm_euclidean_quotient); - else + scm_t_inum xx1 = xx; + scm_t_inum qq; + if (SCM_LIKELY (yy > 0)) { - scm_t_inum xx = SCM_I_INUM (x); - scm_t_inum qq = xx / yy; - if (xx < qq * yy) - { - if (yy > 0) - qq--; - else - qq++; - } - return SCM_I_MAKINUM (qq); + if (SCM_UNLIKELY (xx < 0)) + xx1 = xx - yy + 1; } + else if (SCM_UNLIKELY (yy == 0)) + scm_num_overflow (s_scm_floor_quotient); + else if (xx > 0) + xx1 = xx - yy - 1; + qq = xx1 / yy; + 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) - return SCM_INUM0; + int sign = mpz_sgn (SCM_I_BIG_MPZ (y)); + scm_remember_upto_here_1 (y); + if (sign > 0) + return SCM_I_MAKINUM ((xx < 0) ? -1 : 0); else - return SCM_I_MAKINUM (- mpz_sgn (SCM_I_BIG_MPZ (y))); + return SCM_I_MAKINUM ((xx > 0) ? -1 : 0); } else if (SCM_REALP (y)) - return scm_i_inexact_euclidean_quotient - (SCM_I_INUM (x), SCM_REAL_VALUE (y)); + return scm_i_inexact_floor_quotient (xx, SCM_REAL_VALUE (y)); else if (SCM_FRACTIONP (y)) - return scm_i_slow_exact_euclidean_quotient (x, y); + return scm_i_exact_rational_floor_quotient (x, y); else - SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient, x, y, SCM_ARG2, - s_scm_euclidean_quotient); + return scm_wta_dispatch_2 (g_scm_floor_quotient, x, y, SCM_ARG2, + s_scm_floor_quotient); } else if (SCM_BIGP (x)) { @@ -1130,7 +1013,9 @@ 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); + scm_num_overflow (s_scm_floor_quotient); + else if (SCM_UNLIKELY (yy == 1)) + return x; else { SCM q = scm_i_mkbig (); @@ -1138,7 +1023,7 @@ SCM_PRIMITIVE_GENERIC (scm_euclidean_quotient, "euclidean-quotient", 2, 0, 0, mpz_fdiv_q_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (x), yy); else { - mpz_fdiv_q_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (x), -yy); + mpz_cdiv_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); @@ -1148,149 +1033,139 @@ SCM_PRIMITIVE_GENERIC (scm_euclidean_quotient, "euclidean-quotient", 2, 0, 0, else if (SCM_BIGP (y)) { SCM q = scm_i_mkbig (); - if (mpz_sgn (SCM_I_BIG_MPZ (y)) > 0) - mpz_fdiv_q (SCM_I_BIG_MPZ (q), - SCM_I_BIG_MPZ (x), - SCM_I_BIG_MPZ (y)); - else - mpz_cdiv_q (SCM_I_BIG_MPZ (q), - SCM_I_BIG_MPZ (x), - SCM_I_BIG_MPZ (y)); + mpz_fdiv_q (SCM_I_BIG_MPZ (q), + SCM_I_BIG_MPZ (x), + SCM_I_BIG_MPZ (y)); scm_remember_upto_here_2 (x, y); return scm_i_normbig (q); } else if (SCM_REALP (y)) - return scm_i_inexact_euclidean_quotient + return scm_i_inexact_floor_quotient (scm_i_big2dbl (x), SCM_REAL_VALUE (y)); else if (SCM_FRACTIONP (y)) - return scm_i_slow_exact_euclidean_quotient (x, y); + return scm_i_exact_rational_floor_quotient (x, y); else - SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient, x, y, SCM_ARG2, - s_scm_euclidean_quotient); + return scm_wta_dispatch_2 (g_scm_floor_quotient, x, y, SCM_ARG2, + s_scm_floor_quotient); } else if (SCM_REALP (x)) { if (SCM_REALP (y) || SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) - return scm_i_inexact_euclidean_quotient + return scm_i_inexact_floor_quotient (SCM_REAL_VALUE (x), scm_to_double (y)); else - SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient, x, y, SCM_ARG2, - s_scm_euclidean_quotient); + return scm_wta_dispatch_2 (g_scm_floor_quotient, x, y, SCM_ARG2, + s_scm_floor_quotient); } else if (SCM_FRACTIONP (x)) { if (SCM_REALP (y)) - return scm_i_inexact_euclidean_quotient + return scm_i_inexact_floor_quotient (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_exact_rational_floor_quotient (x, y); else - return scm_i_slow_exact_euclidean_quotient (x, y); + return scm_wta_dispatch_2 (g_scm_floor_quotient, x, y, SCM_ARG2, + s_scm_floor_quotient); } else - SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient, x, y, SCM_ARG1, - s_scm_euclidean_quotient); + return scm_wta_dispatch_2 (g_scm_floor_quotient, x, y, SCM_ARG1, + s_scm_floor_quotient); } #undef FUNC_NAME static SCM -scm_i_inexact_euclidean_quotient (double x, double y) +scm_i_inexact_floor_quotient (double x, double y) { - if (SCM_LIKELY (y > 0)) - return scm_from_double (floor (x / y)); - else if (SCM_LIKELY (y < 0)) - return scm_from_double (ceil (x / y)); - else if (y == 0) - scm_num_overflow (s_scm_euclidean_quotient); /* or return a NaN? */ + if (SCM_UNLIKELY (y == 0)) + scm_num_overflow (s_scm_floor_quotient); /* or return a NaN? */ else - return scm_nan (); + return scm_from_double (floor (x / y)); } -/* Compute exact euclidean_quotient the slow way. - We use this only if both arguments are exact, - and at least one of them is a fraction */ static SCM -scm_i_slow_exact_euclidean_quotient (SCM x, SCM y) +scm_i_exact_rational_floor_quotient (SCM x, SCM y) { - if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x))) - SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient, x, y, SCM_ARG1, - s_scm_euclidean_quotient); - else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y))) - SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient, x, y, SCM_ARG2, - s_scm_euclidean_quotient); - else if (scm_is_true (scm_positive_p (y))) - return scm_floor (scm_divide (x, y)); - else if (scm_is_true (scm_negative_p (y))) - return scm_ceiling (scm_divide (x, y)); - else - scm_num_overflow (s_scm_euclidean_quotient); + return scm_floor_quotient + (scm_product (scm_numerator (x), scm_denominator (y)), + scm_product (scm_numerator (y), scm_denominator (x))); } -static SCM scm_i_inexact_euclidean_remainder (double x, double y); -static SCM scm_i_slow_exact_euclidean_remainder (SCM x, SCM y); +static SCM scm_i_inexact_floor_remainder (double x, double y); +static SCM scm_i_exact_rational_floor_remainder (SCM x, SCM y); -SCM_PRIMITIVE_GENERIC (scm_euclidean_remainder, "euclidean-remainder", 2, 0, 0, +SCM_PRIMITIVE_GENERIC (scm_floor_remainder, "floor-remainder", 2, 0, 0, (SCM x, SCM y), "Return the real number @var{r} such that\n" - "@math{0 <= @var{r} < abs(@var{y})} and\n" "@math{@var{x} = @var{q}*@var{y} + @var{r}}\n" - "for some integer @var{q}.\n" + "where @math{@var{q} = floor(@var{x} / @var{y})}.\n" "@lisp\n" - "(euclidean-remainder 123 10) @result{} 3\n" - "(euclidean-remainder 123 -10) @result{} 3\n" - "(euclidean-remainder -123 10) @result{} 7\n" - "(euclidean-remainder -123 -10) @result{} 7\n" - "(euclidean-remainder -123.2 -63.5) @result{} 3.8\n" - "(euclidean-remainder 16/3 -10/7) @result{} 22/21\n" + "(floor-remainder 123 10) @result{} 3\n" + "(floor-remainder 123 -10) @result{} -7\n" + "(floor-remainder -123 10) @result{} 7\n" + "(floor-remainder -123 -10) @result{} -3\n" + "(floor-remainder -123.2 -63.5) @result{} -59.7\n" + "(floor-remainder 16/3 -10/7) @result{} -8/21\n" "@end lisp") -#define FUNC_NAME s_scm_euclidean_remainder +#define FUNC_NAME s_scm_floor_remainder { 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); if (SCM_UNLIKELY (yy == 0)) - scm_num_overflow (s_scm_euclidean_remainder); + scm_num_overflow (s_scm_floor_remainder); else { - scm_t_inum rr = SCM_I_INUM (x) % yy; - if (rr >= 0) - return SCM_I_MAKINUM (rr); - else if (yy > 0) - return SCM_I_MAKINUM (rr + yy); + scm_t_inum rr = xx % yy; + int needs_adjustment; + + if (SCM_LIKELY (yy > 0)) + needs_adjustment = (rr < 0); else - return SCM_I_MAKINUM (rr - yy); + needs_adjustment = (rr > 0); + + if (needs_adjustment) + rr += yy; + return SCM_I_MAKINUM (rr); } } 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) + int sign = mpz_sgn (SCM_I_BIG_MPZ (y)); + scm_remember_upto_here_1 (y); + if (sign > 0) { - SCM r = scm_i_mkbig (); - mpz_sub_ui (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (y), -xx); - scm_remember_upto_here_1 (y); - return scm_i_normbig (r); + if (xx < 0) + { + SCM r = scm_i_mkbig (); + mpz_sub_ui (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (y), -xx); + scm_remember_upto_here_1 (y); + return scm_i_normbig (r); + } + else + return x; } + else if (xx <= 0) + return x; else { SCM r = scm_i_mkbig (); - mpz_add_ui (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (y), -xx); + mpz_add_ui (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (y), xx); scm_remember_upto_here_1 (y); - mpz_neg (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (r)); return scm_i_normbig (r); } } else if (SCM_REALP (y)) - return scm_i_inexact_euclidean_remainder - (SCM_I_INUM (x), SCM_REAL_VALUE (y)); + return scm_i_inexact_floor_remainder (xx, SCM_REAL_VALUE (y)); else if (SCM_FRACTIONP (y)) - return scm_i_slow_exact_euclidean_remainder (x, y); + return scm_i_exact_rational_floor_remainder (x, y); else - SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder, x, y, SCM_ARG2, - s_scm_euclidean_remainder); + return scm_wta_dispatch_2 (g_scm_floor_remainder, x, y, SCM_ARG2, + s_scm_floor_remainder); } else if (SCM_BIGP (x)) { @@ -1298,13 +1173,14 @@ SCM_PRIMITIVE_GENERIC (scm_euclidean_remainder, "euclidean-remainder", 2, 0, 0, { scm_t_inum yy = SCM_I_INUM (y); if (SCM_UNLIKELY (yy == 0)) - scm_num_overflow (s_scm_euclidean_remainder); + scm_num_overflow (s_scm_floor_remainder); else { scm_t_inum rr; - if (yy < 0) - yy = -yy; - rr = mpz_fdiv_ui (SCM_I_BIG_MPZ (x), yy); + if (yy > 0) + rr = mpz_fdiv_ui (SCM_I_BIG_MPZ (x), yy); + else + rr = -mpz_cdiv_ui (SCM_I_BIG_MPZ (x), -yy); scm_remember_upto_here_1 (x); return SCM_I_MAKINUM (rr); } @@ -1312,165 +1188,184 @@ SCM_PRIMITIVE_GENERIC (scm_euclidean_remainder, "euclidean-remainder", 2, 0, 0, else if (SCM_BIGP (y)) { SCM r = scm_i_mkbig (); - mpz_mod (SCM_I_BIG_MPZ (r), - SCM_I_BIG_MPZ (x), - SCM_I_BIG_MPZ (y)); + mpz_fdiv_r (SCM_I_BIG_MPZ (r), + SCM_I_BIG_MPZ (x), + SCM_I_BIG_MPZ (y)); scm_remember_upto_here_2 (x, y); return scm_i_normbig (r); } else if (SCM_REALP (y)) - return scm_i_inexact_euclidean_remainder + return scm_i_inexact_floor_remainder (scm_i_big2dbl (x), SCM_REAL_VALUE (y)); else if (SCM_FRACTIONP (y)) - return scm_i_slow_exact_euclidean_remainder (x, y); + return scm_i_exact_rational_floor_remainder (x, y); else - SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder, x, y, SCM_ARG2, - s_scm_euclidean_remainder); + return scm_wta_dispatch_2 (g_scm_floor_remainder, x, y, SCM_ARG2, + s_scm_floor_remainder); } else if (SCM_REALP (x)) { if (SCM_REALP (y) || SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) - return scm_i_inexact_euclidean_remainder + return scm_i_inexact_floor_remainder (SCM_REAL_VALUE (x), scm_to_double (y)); else - SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder, x, y, SCM_ARG2, - s_scm_euclidean_remainder); + return scm_wta_dispatch_2 (g_scm_floor_remainder, x, y, SCM_ARG2, + s_scm_floor_remainder); } else if (SCM_FRACTIONP (x)) { if (SCM_REALP (y)) - return scm_i_inexact_euclidean_remainder + return scm_i_inexact_floor_remainder (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_exact_rational_floor_remainder (x, y); else - return scm_i_slow_exact_euclidean_remainder (x, y); + return scm_wta_dispatch_2 (g_scm_floor_remainder, x, y, SCM_ARG2, + s_scm_floor_remainder); } else - SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder, x, y, SCM_ARG1, - s_scm_euclidean_remainder); + return scm_wta_dispatch_2 (g_scm_floor_remainder, x, y, SCM_ARG1, + s_scm_floor_remainder); } #undef FUNC_NAME static SCM -scm_i_inexact_euclidean_remainder (double x, double y) +scm_i_inexact_floor_remainder (double x, double y) { - double q; - /* Although it would be more efficient to use fmod here, we can't because it would in some cases produce results inconsistent with - scm_i_inexact_euclidean_quotient, such that x != q * y + r (not - even close). In particular, when x is very close to a multiple of - y, then r might be either 0.0 or abs(y)-epsilon, but those two - cases must correspond to different choices of q. If r = 0.0 then q - must be x/y, and if r = abs(y) then q must be (x-r)/y. If quotient - chooses one and remainder chooses the other, it would be bad. This - problem was observed with x = 130.0 and y = 10/7. */ - if (SCM_LIKELY (y > 0)) - q = floor (x / y); - else if (SCM_LIKELY (y < 0)) - q = ceil (x / y); - else if (y == 0) - scm_num_overflow (s_scm_euclidean_remainder); /* or return a NaN? */ + scm_i_inexact_floor_quotient, such that x != q * y + r (not even + close). In particular, when x is very close to a multiple of y, + then r might be either 0.0 or y, but those two cases must + correspond to different choices of q. If r = 0.0 then q must be + x/y, and if r = y then q must be x/y-1. If quotient chooses one + and remainder chooses the other, it would be bad. */ + if (SCM_UNLIKELY (y == 0)) + scm_num_overflow (s_scm_floor_remainder); /* or return a NaN? */ else - return scm_nan (); - return scm_from_double (x - q * y); + return scm_from_double (x - y * floor (x / y)); } -/* Compute exact euclidean_remainder the slow way. - We use this only if both arguments are exact, - and at least one of them is a fraction */ static SCM -scm_i_slow_exact_euclidean_remainder (SCM x, SCM y) +scm_i_exact_rational_floor_remainder (SCM x, SCM y) { - SCM q; - - if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x))) - SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder, x, y, SCM_ARG1, - s_scm_euclidean_remainder); - else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y))) - SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder, x, y, SCM_ARG2, - s_scm_euclidean_remainder); - else if (scm_is_true (scm_positive_p (y))) - q = scm_floor (scm_divide (x, y)); - else if (scm_is_true (scm_negative_p (y))) - q = scm_ceiling (scm_divide (x, y)); - else - scm_num_overflow (s_scm_euclidean_remainder); - return scm_difference (x, scm_product (y, q)); + SCM xd = scm_denominator (x); + SCM yd = scm_denominator (y); + SCM r1 = scm_floor_remainder (scm_product (scm_numerator (x), yd), + scm_product (scm_numerator (y), xd)); + return scm_divide (r1, scm_product (xd, yd)); } -static SCM scm_i_inexact_euclidean_divide (double x, double y); -static SCM scm_i_slow_exact_euclidean_divide (SCM x, SCM y); +static void scm_i_inexact_floor_divide (double x, double y, + SCM *qp, SCM *rp); +static void scm_i_exact_rational_floor_divide (SCM x, SCM y, + SCM *qp, SCM *rp); -SCM_PRIMITIVE_GENERIC (scm_euclidean_divide, "euclidean/", 2, 0, 0, +SCM_PRIMITIVE_GENERIC (scm_i_floor_divide, "floor/", 2, 0, 0, (SCM x, SCM y), "Return the integer @var{q} and the real number @var{r}\n" "such that @math{@var{x} = @var{q}*@var{y} + @var{r}}\n" - "and @math{0 <= @var{r} < abs(@var{y})}.\n" + "and @math{@var{q} = floor(@var{x} / @var{y})}.\n" "@lisp\n" - "(euclidean/ 123 10) @result{} 12 and 3\n" - "(euclidean/ 123 -10) @result{} -12 and 3\n" - "(euclidean/ -123 10) @result{} -13 and 7\n" - "(euclidean/ -123 -10) @result{} 13 and 7\n" - "(euclidean/ -123.2 -63.5) @result{} 2.0 and 3.8\n" - "(euclidean/ 16/3 -10/7) @result{} -3 and 22/21\n" + "(floor/ 123 10) @result{} 12 and 3\n" + "(floor/ 123 -10) @result{} -13 and -7\n" + "(floor/ -123 10) @result{} -13 and 7\n" + "(floor/ -123 -10) @result{} 12 and -3\n" + "(floor/ -123.2 -63.5) @result{} 1.0 and -59.7\n" + "(floor/ 16/3 -10/7) @result{} -4 and -8/21\n" "@end lisp") -#define FUNC_NAME s_scm_euclidean_divide +#define FUNC_NAME s_scm_i_floor_divide { - if (SCM_LIKELY (SCM_I_INUMP (x))) - { - if (SCM_LIKELY (SCM_I_INUMP (y))) - { - scm_t_inum yy = SCM_I_INUM (y); - if (SCM_UNLIKELY (yy == 0)) - scm_num_overflow (s_scm_euclidean_divide); - else + SCM q, r; + + scm_floor_divide(x, y, &q, &r); + return scm_values (scm_list_2 (q, r)); +} +#undef FUNC_NAME + +#define s_scm_floor_divide s_scm_i_floor_divide +#define g_scm_floor_divide g_scm_i_floor_divide + +void +scm_floor_divide (SCM x, SCM y, SCM *qp, SCM *rp) +{ + 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); + if (SCM_UNLIKELY (yy == 0)) + scm_num_overflow (s_scm_floor_divide); + else { - scm_t_inum xx = SCM_I_INUM (x); scm_t_inum qq = xx / yy; - scm_t_inum rr = xx - qq * yy; - if (rr < 0) + scm_t_inum rr = xx % yy; + int needs_adjustment; + + if (SCM_LIKELY (yy > 0)) + needs_adjustment = (rr < 0); + else + needs_adjustment = (rr > 0); + + if (needs_adjustment) { - if (yy > 0) - { rr += yy; qq--; } - else - { rr -= yy; qq++; } + rr += yy; + qq--; } - return scm_values (scm_list_2 (SCM_I_MAKINUM (qq), - SCM_I_MAKINUM (rr))); + + if (SCM_LIKELY (SCM_FIXABLE (qq))) + *qp = SCM_I_MAKINUM (qq); + else + *qp = scm_i_inum2big (qq); + *rp = SCM_I_MAKINUM (rr); } + return; } 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) + int sign = mpz_sgn (SCM_I_BIG_MPZ (y)); + scm_remember_upto_here_1 (y); + if (sign > 0) { - SCM r = scm_i_mkbig (); - mpz_sub_ui (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (y), -xx); - scm_remember_upto_here_1 (y); - return scm_values - (scm_list_2 (SCM_I_MAKINUM (-1), scm_i_normbig (r))); + if (xx < 0) + { + SCM r = scm_i_mkbig (); + mpz_sub_ui (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (y), -xx); + scm_remember_upto_here_1 (y); + *qp = SCM_I_MAKINUM (-1); + *rp = scm_i_normbig (r); + } + else + { + *qp = SCM_INUM0; + *rp = x; + } + } + else if (xx <= 0) + { + *qp = SCM_INUM0; + *rp = x; } else { SCM r = scm_i_mkbig (); - mpz_add_ui (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (y), -xx); + mpz_add_ui (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (y), xx); scm_remember_upto_here_1 (y); - mpz_neg (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (r)); - return scm_values (scm_list_2 (SCM_INUM1, scm_i_normbig (r))); + *qp = SCM_I_MAKINUM (-1); + *rp = scm_i_normbig (r); } + return; } else if (SCM_REALP (y)) - return scm_i_inexact_euclidean_divide - (SCM_I_INUM (x), SCM_REAL_VALUE (y)); + return scm_i_inexact_floor_divide (xx, SCM_REAL_VALUE (y), qp, rp); else if (SCM_FRACTIONP (y)) - return scm_i_slow_exact_euclidean_divide (x, y); + return scm_i_exact_rational_floor_divide (x, y, qp, rp); else - SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide, x, y, SCM_ARG2, - s_scm_euclidean_divide); + return two_valued_wta_dispatch_2 (g_scm_floor_divide, x, y, SCM_ARG2, + s_scm_floor_divide, qp, rp); } else if (SCM_BIGP (x)) { @@ -1478,7 +1373,7 @@ SCM_PRIMITIVE_GENERIC (scm_euclidean_divide, "euclidean/", 2, 0, 0, { scm_t_inum yy = SCM_I_INUM (y); if (SCM_UNLIKELY (yy == 0)) - scm_num_overflow (s_scm_euclidean_divide); + scm_num_overflow (s_scm_floor_divide); else { SCM q = scm_i_mkbig (); @@ -1488,179 +1383,1874 @@ SCM_PRIMITIVE_GENERIC (scm_euclidean_divide, "euclidean/", 2, 0, 0, SCM_I_BIG_MPZ (x), yy); else { - mpz_fdiv_qr_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r), + mpz_cdiv_qr_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r), 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))); + *qp = scm_i_normbig (q); + *rp = scm_i_normbig (r); } + return; } else if (SCM_BIGP (y)) { SCM q = scm_i_mkbig (); SCM r = scm_i_mkbig (); - if (mpz_sgn (SCM_I_BIG_MPZ (y)) > 0) - mpz_fdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r), - SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y)); + mpz_fdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r), + SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y)); + scm_remember_upto_here_2 (x, y); + *qp = scm_i_normbig (q); + *rp = scm_i_normbig (r); + return; + } + else if (SCM_REALP (y)) + return scm_i_inexact_floor_divide + (scm_i_big2dbl (x), SCM_REAL_VALUE (y), qp, rp); + else if (SCM_FRACTIONP (y)) + return scm_i_exact_rational_floor_divide (x, y, qp, rp); + else + return two_valued_wta_dispatch_2 (g_scm_floor_divide, x, y, SCM_ARG2, + s_scm_floor_divide, qp, rp); + } + else if (SCM_REALP (x)) + { + if (SCM_REALP (y) || SCM_I_INUMP (y) || + SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_inexact_floor_divide + (SCM_REAL_VALUE (x), scm_to_double (y), qp, rp); + else + return two_valued_wta_dispatch_2 (g_scm_floor_divide, x, y, SCM_ARG2, + s_scm_floor_divide, qp, rp); + } + else if (SCM_FRACTIONP (x)) + { + if (SCM_REALP (y)) + return scm_i_inexact_floor_divide + (scm_i_fraction2double (x), SCM_REAL_VALUE (y), qp, rp); + else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_exact_rational_floor_divide (x, y, qp, rp); + else + return two_valued_wta_dispatch_2 (g_scm_floor_divide, x, y, SCM_ARG2, + s_scm_floor_divide, qp, rp); + } + else + return two_valued_wta_dispatch_2 (g_scm_floor_divide, x, y, SCM_ARG1, + s_scm_floor_divide, qp, rp); +} + +static void +scm_i_inexact_floor_divide (double x, double y, SCM *qp, SCM *rp) +{ + if (SCM_UNLIKELY (y == 0)) + scm_num_overflow (s_scm_floor_divide); /* or return a NaN? */ + else + { + double q = floor (x / y); + double r = x - q * y; + *qp = scm_from_double (q); + *rp = scm_from_double (r); + } +} + +static void +scm_i_exact_rational_floor_divide (SCM x, SCM y, SCM *qp, SCM *rp) +{ + SCM r1; + SCM xd = scm_denominator (x); + SCM yd = scm_denominator (y); + + scm_floor_divide (scm_product (scm_numerator (x), yd), + scm_product (scm_numerator (y), xd), + qp, &r1); + *rp = scm_divide (r1, scm_product (xd, yd)); +} + +static SCM scm_i_inexact_ceiling_quotient (double x, double y); +static SCM scm_i_exact_rational_ceiling_quotient (SCM x, SCM y); + +SCM_PRIMITIVE_GENERIC (scm_ceiling_quotient, "ceiling-quotient", 2, 0, 0, + (SCM x, SCM y), + "Return the ceiling of @math{@var{x} / @var{y}}.\n" + "@lisp\n" + "(ceiling-quotient 123 10) @result{} 13\n" + "(ceiling-quotient 123 -10) @result{} -12\n" + "(ceiling-quotient -123 10) @result{} -12\n" + "(ceiling-quotient -123 -10) @result{} 13\n" + "(ceiling-quotient -123.2 -63.5) @result{} 2.0\n" + "(ceiling-quotient 16/3 -10/7) @result{} -3\n" + "@end lisp") +#define FUNC_NAME s_scm_ceiling_quotient +{ + 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); + if (SCM_UNLIKELY (yy == 0)) + scm_num_overflow (s_scm_ceiling_quotient); + else + { + scm_t_inum xx1 = xx; + scm_t_inum qq; + if (SCM_LIKELY (yy > 0)) + { + if (SCM_LIKELY (xx >= 0)) + xx1 = xx + yy - 1; + } + else if (SCM_UNLIKELY (yy == 0)) + scm_num_overflow (s_scm_ceiling_quotient); + else if (xx < 0) + xx1 = xx + yy + 1; + qq = xx1 / yy; + if (SCM_LIKELY (SCM_FIXABLE (qq))) + return SCM_I_MAKINUM (qq); + else + return scm_i_inum2big (qq); + } + } + else if (SCM_BIGP (y)) + { + int sign = mpz_sgn (SCM_I_BIG_MPZ (y)); + scm_remember_upto_here_1 (y); + if (SCM_LIKELY (sign > 0)) + { + if (SCM_LIKELY (xx > 0)) + return SCM_INUM1; + else if (SCM_UNLIKELY (xx == SCM_MOST_NEGATIVE_FIXNUM) + && SCM_UNLIKELY (mpz_cmp_ui (SCM_I_BIG_MPZ (y), + - SCM_MOST_NEGATIVE_FIXNUM) == 0)) + { + /* Special case: x == fixnum-min && y == abs (fixnum-min) */ + scm_remember_upto_here_1 (y); + return SCM_I_MAKINUM (-1); + } + else + return SCM_INUM0; + } + else if (xx >= 0) + return SCM_INUM0; + else + return SCM_INUM1; + } + else if (SCM_REALP (y)) + return scm_i_inexact_ceiling_quotient (xx, SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_exact_rational_ceiling_quotient (x, y); + else + return scm_wta_dispatch_2 (g_scm_ceiling_quotient, x, y, SCM_ARG2, + s_scm_ceiling_quotient); + } + else if (SCM_BIGP (x)) + { + if (SCM_LIKELY (SCM_I_INUMP (y))) + { + scm_t_inum yy = SCM_I_INUM (y); + if (SCM_UNLIKELY (yy == 0)) + scm_num_overflow (s_scm_ceiling_quotient); + else if (SCM_UNLIKELY (yy == 1)) + return x; + else + { + SCM q = scm_i_mkbig (); + if (yy > 0) + mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (x), yy); + else + { + 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_i_normbig (q); + } + } + else if (SCM_BIGP (y)) + { + SCM q = scm_i_mkbig (); + mpz_cdiv_q (SCM_I_BIG_MPZ (q), + SCM_I_BIG_MPZ (x), + SCM_I_BIG_MPZ (y)); + scm_remember_upto_here_2 (x, y); + return scm_i_normbig (q); + } + else if (SCM_REALP (y)) + return scm_i_inexact_ceiling_quotient + (scm_i_big2dbl (x), SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_exact_rational_ceiling_quotient (x, y); + else + return scm_wta_dispatch_2 (g_scm_ceiling_quotient, x, y, SCM_ARG2, + s_scm_ceiling_quotient); + } + else if (SCM_REALP (x)) + { + if (SCM_REALP (y) || SCM_I_INUMP (y) || + SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_inexact_ceiling_quotient + (SCM_REAL_VALUE (x), scm_to_double (y)); + else + return scm_wta_dispatch_2 (g_scm_ceiling_quotient, x, y, SCM_ARG2, + s_scm_ceiling_quotient); + } + else if (SCM_FRACTIONP (x)) + { + if (SCM_REALP (y)) + return scm_i_inexact_ceiling_quotient + (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_exact_rational_ceiling_quotient (x, y); + else + return scm_wta_dispatch_2 (g_scm_ceiling_quotient, x, y, SCM_ARG2, + s_scm_ceiling_quotient); + } + else + return scm_wta_dispatch_2 (g_scm_ceiling_quotient, x, y, SCM_ARG1, + s_scm_ceiling_quotient); +} +#undef FUNC_NAME + +static SCM +scm_i_inexact_ceiling_quotient (double x, double y) +{ + if (SCM_UNLIKELY (y == 0)) + scm_num_overflow (s_scm_ceiling_quotient); /* or return a NaN? */ + else + return scm_from_double (ceil (x / y)); +} + +static SCM +scm_i_exact_rational_ceiling_quotient (SCM x, SCM y) +{ + return scm_ceiling_quotient + (scm_product (scm_numerator (x), scm_denominator (y)), + scm_product (scm_numerator (y), scm_denominator (x))); +} + +static SCM scm_i_inexact_ceiling_remainder (double x, double y); +static SCM scm_i_exact_rational_ceiling_remainder (SCM x, SCM y); + +SCM_PRIMITIVE_GENERIC (scm_ceiling_remainder, "ceiling-remainder", 2, 0, 0, + (SCM x, SCM y), + "Return the real number @var{r} such that\n" + "@math{@var{x} = @var{q}*@var{y} + @var{r}}\n" + "where @math{@var{q} = ceiling(@var{x} / @var{y})}.\n" + "@lisp\n" + "(ceiling-remainder 123 10) @result{} -7\n" + "(ceiling-remainder 123 -10) @result{} 3\n" + "(ceiling-remainder -123 10) @result{} -3\n" + "(ceiling-remainder -123 -10) @result{} 7\n" + "(ceiling-remainder -123.2 -63.5) @result{} 3.8\n" + "(ceiling-remainder 16/3 -10/7) @result{} 22/21\n" + "@end lisp") +#define FUNC_NAME s_scm_ceiling_remainder +{ + 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); + if (SCM_UNLIKELY (yy == 0)) + scm_num_overflow (s_scm_ceiling_remainder); + else + { + scm_t_inum rr = xx % yy; + int needs_adjustment; + + if (SCM_LIKELY (yy > 0)) + needs_adjustment = (rr > 0); + else + needs_adjustment = (rr < 0); + + if (needs_adjustment) + rr -= yy; + return SCM_I_MAKINUM (rr); + } + } + else if (SCM_BIGP (y)) + { + int sign = mpz_sgn (SCM_I_BIG_MPZ (y)); + scm_remember_upto_here_1 (y); + if (SCM_LIKELY (sign > 0)) + { + if (SCM_LIKELY (xx > 0)) + { + SCM r = scm_i_mkbig (); + mpz_sub_ui (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (y), xx); + scm_remember_upto_here_1 (y); + mpz_neg (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (r)); + return scm_i_normbig (r); + } + else if (SCM_UNLIKELY (xx == SCM_MOST_NEGATIVE_FIXNUM) + && SCM_UNLIKELY (mpz_cmp_ui (SCM_I_BIG_MPZ (y), + - SCM_MOST_NEGATIVE_FIXNUM) == 0)) + { + /* Special case: x == fixnum-min && y == abs (fixnum-min) */ + scm_remember_upto_here_1 (y); + return SCM_INUM0; + } + else + return x; + } + else if (xx >= 0) + return x; + else + { + SCM r = scm_i_mkbig (); + mpz_add_ui (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (y), -xx); + scm_remember_upto_here_1 (y); + mpz_neg (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (r)); + return scm_i_normbig (r); + } + } + else if (SCM_REALP (y)) + return scm_i_inexact_ceiling_remainder (xx, SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_exact_rational_ceiling_remainder (x, y); + else + return scm_wta_dispatch_2 (g_scm_ceiling_remainder, x, y, SCM_ARG2, + s_scm_ceiling_remainder); + } + else if (SCM_BIGP (x)) + { + if (SCM_LIKELY (SCM_I_INUMP (y))) + { + scm_t_inum yy = SCM_I_INUM (y); + if (SCM_UNLIKELY (yy == 0)) + scm_num_overflow (s_scm_ceiling_remainder); else - mpz_cdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r), - SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y)); + { + scm_t_inum rr; + if (yy > 0) + rr = -mpz_cdiv_ui (SCM_I_BIG_MPZ (x), yy); + else + rr = mpz_fdiv_ui (SCM_I_BIG_MPZ (x), -yy); + scm_remember_upto_here_1 (x); + return SCM_I_MAKINUM (rr); + } + } + else if (SCM_BIGP (y)) + { + SCM r = scm_i_mkbig (); + mpz_cdiv_r (SCM_I_BIG_MPZ (r), + SCM_I_BIG_MPZ (x), + SCM_I_BIG_MPZ (y)); scm_remember_upto_here_2 (x, y); - return scm_values (scm_list_2 (scm_i_normbig (q), - scm_i_normbig (r))); + return scm_i_normbig (r); } else if (SCM_REALP (y)) - return scm_i_inexact_euclidean_divide + return scm_i_inexact_ceiling_remainder (scm_i_big2dbl (x), SCM_REAL_VALUE (y)); else if (SCM_FRACTIONP (y)) - return scm_i_slow_exact_euclidean_divide (x, y); + return scm_i_exact_rational_ceiling_remainder (x, y); else - SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide, x, y, SCM_ARG2, - s_scm_euclidean_divide); + return scm_wta_dispatch_2 (g_scm_ceiling_remainder, x, y, SCM_ARG2, + s_scm_ceiling_remainder); } else if (SCM_REALP (x)) { if (SCM_REALP (y) || SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) - return scm_i_inexact_euclidean_divide + return scm_i_inexact_ceiling_remainder (SCM_REAL_VALUE (x), scm_to_double (y)); else - SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide, x, y, SCM_ARG2, - s_scm_euclidean_divide); + return scm_wta_dispatch_2 (g_scm_ceiling_remainder, x, y, SCM_ARG2, + s_scm_ceiling_remainder); } else if (SCM_FRACTIONP (x)) { if (SCM_REALP (y)) - return scm_i_inexact_euclidean_divide + return scm_i_inexact_ceiling_remainder (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_exact_rational_ceiling_remainder (x, y); else - return scm_i_slow_exact_euclidean_divide (x, y); + return scm_wta_dispatch_2 (g_scm_ceiling_remainder, x, y, SCM_ARG2, + s_scm_ceiling_remainder); } else - SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide, x, y, SCM_ARG1, - s_scm_euclidean_divide); + return scm_wta_dispatch_2 (g_scm_ceiling_remainder, x, y, SCM_ARG1, + s_scm_ceiling_remainder); } #undef FUNC_NAME static SCM -scm_i_inexact_euclidean_divide (double x, double y) +scm_i_inexact_ceiling_remainder (double x, double y) { - double q, r; + /* Although it would be more efficient to use fmod here, we can't + because it would in some cases produce results inconsistent with + scm_i_inexact_ceiling_quotient, such that x != q * y + r (not even + close). In particular, when x is very close to a multiple of y, + then r might be either 0.0 or -y, but those two cases must + correspond to different choices of q. If r = 0.0 then q must be + x/y, and if r = -y then q must be x/y+1. If quotient chooses one + and remainder chooses the other, it would be bad. */ + if (SCM_UNLIKELY (y == 0)) + scm_num_overflow (s_scm_ceiling_remainder); /* or return a NaN? */ + else + return scm_from_double (x - y * ceil (x / y)); +} - if (SCM_LIKELY (y > 0)) - q = floor (x / y); - else if (SCM_LIKELY (y < 0)) - q = ceil (x / y); - else if (y == 0) - scm_num_overflow (s_scm_euclidean_divide); /* or return a NaN? */ +static SCM +scm_i_exact_rational_ceiling_remainder (SCM x, SCM y) +{ + SCM xd = scm_denominator (x); + SCM yd = scm_denominator (y); + SCM r1 = scm_ceiling_remainder (scm_product (scm_numerator (x), yd), + scm_product (scm_numerator (y), xd)); + return scm_divide (r1, scm_product (xd, yd)); +} + +static void scm_i_inexact_ceiling_divide (double x, double y, + SCM *qp, SCM *rp); +static void scm_i_exact_rational_ceiling_divide (SCM x, SCM y, + SCM *qp, SCM *rp); + +SCM_PRIMITIVE_GENERIC (scm_i_ceiling_divide, "ceiling/", 2, 0, 0, + (SCM x, SCM y), + "Return the integer @var{q} and the real number @var{r}\n" + "such that @math{@var{x} = @var{q}*@var{y} + @var{r}}\n" + "and @math{@var{q} = ceiling(@var{x} / @var{y})}.\n" + "@lisp\n" + "(ceiling/ 123 10) @result{} 13 and -7\n" + "(ceiling/ 123 -10) @result{} -12 and 3\n" + "(ceiling/ -123 10) @result{} -12 and -3\n" + "(ceiling/ -123 -10) @result{} 13 and 7\n" + "(ceiling/ -123.2 -63.5) @result{} 2.0 and 3.8\n" + "(ceiling/ 16/3 -10/7) @result{} -3 and 22/21\n" + "@end lisp") +#define FUNC_NAME s_scm_i_ceiling_divide +{ + SCM q, r; + + scm_ceiling_divide(x, y, &q, &r); + return scm_values (scm_list_2 (q, r)); +} +#undef FUNC_NAME + +#define s_scm_ceiling_divide s_scm_i_ceiling_divide +#define g_scm_ceiling_divide g_scm_i_ceiling_divide + +void +scm_ceiling_divide (SCM x, SCM y, SCM *qp, SCM *rp) +{ + 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); + if (SCM_UNLIKELY (yy == 0)) + scm_num_overflow (s_scm_ceiling_divide); + else + { + scm_t_inum qq = xx / yy; + scm_t_inum rr = xx % yy; + int needs_adjustment; + + if (SCM_LIKELY (yy > 0)) + needs_adjustment = (rr > 0); + else + needs_adjustment = (rr < 0); + + if (needs_adjustment) + { + rr -= yy; + qq++; + } + if (SCM_LIKELY (SCM_FIXABLE (qq))) + *qp = SCM_I_MAKINUM (qq); + else + *qp = scm_i_inum2big (qq); + *rp = SCM_I_MAKINUM (rr); + } + return; + } + else if (SCM_BIGP (y)) + { + int sign = mpz_sgn (SCM_I_BIG_MPZ (y)); + scm_remember_upto_here_1 (y); + if (SCM_LIKELY (sign > 0)) + { + if (SCM_LIKELY (xx > 0)) + { + SCM r = scm_i_mkbig (); + mpz_sub_ui (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (y), xx); + scm_remember_upto_here_1 (y); + mpz_neg (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (r)); + *qp = SCM_INUM1; + *rp = scm_i_normbig (r); + } + else if (SCM_UNLIKELY (xx == SCM_MOST_NEGATIVE_FIXNUM) + && SCM_UNLIKELY (mpz_cmp_ui (SCM_I_BIG_MPZ (y), + - SCM_MOST_NEGATIVE_FIXNUM) == 0)) + { + /* Special case: x == fixnum-min && y == abs (fixnum-min) */ + scm_remember_upto_here_1 (y); + *qp = SCM_I_MAKINUM (-1); + *rp = SCM_INUM0; + } + else + { + *qp = SCM_INUM0; + *rp = x; + } + } + else if (xx >= 0) + { + *qp = SCM_INUM0; + *rp = x; + } + else + { + SCM r = scm_i_mkbig (); + mpz_add_ui (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (y), -xx); + scm_remember_upto_here_1 (y); + mpz_neg (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (r)); + *qp = SCM_INUM1; + *rp = scm_i_normbig (r); + } + return; + } + else if (SCM_REALP (y)) + return scm_i_inexact_ceiling_divide (xx, SCM_REAL_VALUE (y), qp, rp); + else if (SCM_FRACTIONP (y)) + return scm_i_exact_rational_ceiling_divide (x, y, qp, rp); + else + return two_valued_wta_dispatch_2 (g_scm_ceiling_divide, x, y, SCM_ARG2, + s_scm_ceiling_divide, qp, rp); + } + else if (SCM_BIGP (x)) + { + if (SCM_LIKELY (SCM_I_INUMP (y))) + { + scm_t_inum yy = SCM_I_INUM (y); + if (SCM_UNLIKELY (yy == 0)) + scm_num_overflow (s_scm_ceiling_divide); + else + { + SCM q = scm_i_mkbig (); + SCM r = scm_i_mkbig (); + if (yy > 0) + mpz_cdiv_qr_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r), + 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); + mpz_neg (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (q)); + } + scm_remember_upto_here_1 (x); + *qp = scm_i_normbig (q); + *rp = scm_i_normbig (r); + } + return; + } + else if (SCM_BIGP (y)) + { + SCM q = scm_i_mkbig (); + SCM r = scm_i_mkbig (); + mpz_cdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r), + SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y)); + scm_remember_upto_here_2 (x, y); + *qp = scm_i_normbig (q); + *rp = scm_i_normbig (r); + return; + } + else if (SCM_REALP (y)) + return scm_i_inexact_ceiling_divide + (scm_i_big2dbl (x), SCM_REAL_VALUE (y), qp, rp); + else if (SCM_FRACTIONP (y)) + return scm_i_exact_rational_ceiling_divide (x, y, qp, rp); + else + return two_valued_wta_dispatch_2 (g_scm_ceiling_divide, x, y, SCM_ARG2, + s_scm_ceiling_divide, qp, rp); + } + else if (SCM_REALP (x)) + { + if (SCM_REALP (y) || SCM_I_INUMP (y) || + SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_inexact_ceiling_divide + (SCM_REAL_VALUE (x), scm_to_double (y), qp, rp); + else + return two_valued_wta_dispatch_2 (g_scm_ceiling_divide, x, y, SCM_ARG2, + s_scm_ceiling_divide, qp, rp); + } + else if (SCM_FRACTIONP (x)) + { + if (SCM_REALP (y)) + return scm_i_inexact_ceiling_divide + (scm_i_fraction2double (x), SCM_REAL_VALUE (y), qp, rp); + else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_exact_rational_ceiling_divide (x, y, qp, rp); + else + return two_valued_wta_dispatch_2 (g_scm_ceiling_divide, x, y, SCM_ARG2, + s_scm_ceiling_divide, qp, rp); + } + else + return two_valued_wta_dispatch_2 (g_scm_ceiling_divide, x, y, SCM_ARG1, + s_scm_ceiling_divide, qp, rp); +} + +static void +scm_i_inexact_ceiling_divide (double x, double y, SCM *qp, SCM *rp) +{ + if (SCM_UNLIKELY (y == 0)) + scm_num_overflow (s_scm_ceiling_divide); /* or return a NaN? */ + else + { + double q = ceil (x / y); + double r = x - q * y; + *qp = scm_from_double (q); + *rp = scm_from_double (r); + } +} + +static void +scm_i_exact_rational_ceiling_divide (SCM x, SCM y, SCM *qp, SCM *rp) +{ + SCM r1; + SCM xd = scm_denominator (x); + SCM yd = scm_denominator (y); + + scm_ceiling_divide (scm_product (scm_numerator (x), yd), + scm_product (scm_numerator (y), xd), + qp, &r1); + *rp = scm_divide (r1, scm_product (xd, yd)); +} + +static SCM scm_i_inexact_truncate_quotient (double x, double y); +static SCM scm_i_exact_rational_truncate_quotient (SCM x, SCM y); + +SCM_PRIMITIVE_GENERIC (scm_truncate_quotient, "truncate-quotient", 2, 0, 0, + (SCM x, SCM y), + "Return @math{@var{x} / @var{y}} rounded toward zero.\n" + "@lisp\n" + "(truncate-quotient 123 10) @result{} 12\n" + "(truncate-quotient 123 -10) @result{} -12\n" + "(truncate-quotient -123 10) @result{} -12\n" + "(truncate-quotient -123 -10) @result{} 12\n" + "(truncate-quotient -123.2 -63.5) @result{} 1.0\n" + "(truncate-quotient 16/3 -10/7) @result{} -3\n" + "@end lisp") +#define FUNC_NAME s_scm_truncate_quotient +{ + 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); + if (SCM_UNLIKELY (yy == 0)) + scm_num_overflow (s_scm_truncate_quotient); + else + { + scm_t_inum qq = xx / yy; + if (SCM_LIKELY (SCM_FIXABLE (qq))) + return SCM_I_MAKINUM (qq); + else + return scm_i_inum2big (qq); + } + } + else if (SCM_BIGP (y)) + { + if (SCM_UNLIKELY (xx == SCM_MOST_NEGATIVE_FIXNUM) + && SCM_UNLIKELY (mpz_cmp_ui (SCM_I_BIG_MPZ (y), + - SCM_MOST_NEGATIVE_FIXNUM) == 0)) + { + /* Special case: x == fixnum-min && y == abs (fixnum-min) */ + scm_remember_upto_here_1 (y); + return SCM_I_MAKINUM (-1); + } + else + return SCM_INUM0; + } + else if (SCM_REALP (y)) + return scm_i_inexact_truncate_quotient (xx, SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_exact_rational_truncate_quotient (x, y); + else + return scm_wta_dispatch_2 (g_scm_truncate_quotient, x, y, SCM_ARG2, + s_scm_truncate_quotient); + } + else if (SCM_BIGP (x)) + { + if (SCM_LIKELY (SCM_I_INUMP (y))) + { + scm_t_inum yy = SCM_I_INUM (y); + if (SCM_UNLIKELY (yy == 0)) + scm_num_overflow (s_scm_truncate_quotient); + else if (SCM_UNLIKELY (yy == 1)) + return x; + else + { + SCM q = scm_i_mkbig (); + if (yy > 0) + mpz_tdiv_q_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (x), yy); + else + { + mpz_tdiv_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_i_normbig (q); + } + } + else if (SCM_BIGP (y)) + { + SCM q = scm_i_mkbig (); + mpz_tdiv_q (SCM_I_BIG_MPZ (q), + SCM_I_BIG_MPZ (x), + SCM_I_BIG_MPZ (y)); + scm_remember_upto_here_2 (x, y); + return scm_i_normbig (q); + } + else if (SCM_REALP (y)) + return scm_i_inexact_truncate_quotient + (scm_i_big2dbl (x), SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_exact_rational_truncate_quotient (x, y); + else + return scm_wta_dispatch_2 (g_scm_truncate_quotient, x, y, SCM_ARG2, + s_scm_truncate_quotient); + } + else if (SCM_REALP (x)) + { + if (SCM_REALP (y) || SCM_I_INUMP (y) || + SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_inexact_truncate_quotient + (SCM_REAL_VALUE (x), scm_to_double (y)); + else + return scm_wta_dispatch_2 (g_scm_truncate_quotient, x, y, SCM_ARG2, + s_scm_truncate_quotient); + } + else if (SCM_FRACTIONP (x)) + { + if (SCM_REALP (y)) + return scm_i_inexact_truncate_quotient + (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_exact_rational_truncate_quotient (x, y); + else + return scm_wta_dispatch_2 (g_scm_truncate_quotient, x, y, SCM_ARG2, + s_scm_truncate_quotient); + } + else + return scm_wta_dispatch_2 (g_scm_truncate_quotient, x, y, SCM_ARG1, + s_scm_truncate_quotient); +} +#undef FUNC_NAME + +static SCM +scm_i_inexact_truncate_quotient (double x, double y) +{ + if (SCM_UNLIKELY (y == 0)) + scm_num_overflow (s_scm_truncate_quotient); /* or return a NaN? */ + else + return scm_from_double (trunc (x / y)); +} + +static SCM +scm_i_exact_rational_truncate_quotient (SCM x, SCM y) +{ + return scm_truncate_quotient + (scm_product (scm_numerator (x), scm_denominator (y)), + scm_product (scm_numerator (y), scm_denominator (x))); +} + +static SCM scm_i_inexact_truncate_remainder (double x, double y); +static SCM scm_i_exact_rational_truncate_remainder (SCM x, SCM y); + +SCM_PRIMITIVE_GENERIC (scm_truncate_remainder, "truncate-remainder", 2, 0, 0, + (SCM x, SCM y), + "Return the real number @var{r} such that\n" + "@math{@var{x} = @var{q}*@var{y} + @var{r}}\n" + "where @math{@var{q} = truncate(@var{x} / @var{y})}.\n" + "@lisp\n" + "(truncate-remainder 123 10) @result{} 3\n" + "(truncate-remainder 123 -10) @result{} 3\n" + "(truncate-remainder -123 10) @result{} -3\n" + "(truncate-remainder -123 -10) @result{} -3\n" + "(truncate-remainder -123.2 -63.5) @result{} -59.7\n" + "(truncate-remainder 16/3 -10/7) @result{} 22/21\n" + "@end lisp") +#define FUNC_NAME s_scm_truncate_remainder +{ + 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); + if (SCM_UNLIKELY (yy == 0)) + scm_num_overflow (s_scm_truncate_remainder); + else + return SCM_I_MAKINUM (xx % yy); + } + else if (SCM_BIGP (y)) + { + if (SCM_UNLIKELY (xx == SCM_MOST_NEGATIVE_FIXNUM) + && SCM_UNLIKELY (mpz_cmp_ui (SCM_I_BIG_MPZ (y), + - SCM_MOST_NEGATIVE_FIXNUM) == 0)) + { + /* Special case: x == fixnum-min && y == abs (fixnum-min) */ + scm_remember_upto_here_1 (y); + return SCM_INUM0; + } + else + return x; + } + else if (SCM_REALP (y)) + return scm_i_inexact_truncate_remainder (xx, SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_exact_rational_truncate_remainder (x, y); + else + return scm_wta_dispatch_2 (g_scm_truncate_remainder, x, y, SCM_ARG2, + s_scm_truncate_remainder); + } + else if (SCM_BIGP (x)) + { + if (SCM_LIKELY (SCM_I_INUMP (y))) + { + scm_t_inum yy = SCM_I_INUM (y); + if (SCM_UNLIKELY (yy == 0)) + scm_num_overflow (s_scm_truncate_remainder); + else + { + scm_t_inum rr = (mpz_tdiv_ui (SCM_I_BIG_MPZ (x), + (yy > 0) ? yy : -yy) + * mpz_sgn (SCM_I_BIG_MPZ (x))); + scm_remember_upto_here_1 (x); + return SCM_I_MAKINUM (rr); + } + } + else if (SCM_BIGP (y)) + { + SCM r = scm_i_mkbig (); + mpz_tdiv_r (SCM_I_BIG_MPZ (r), + SCM_I_BIG_MPZ (x), + SCM_I_BIG_MPZ (y)); + scm_remember_upto_here_2 (x, y); + return scm_i_normbig (r); + } + else if (SCM_REALP (y)) + return scm_i_inexact_truncate_remainder + (scm_i_big2dbl (x), SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_exact_rational_truncate_remainder (x, y); + else + return scm_wta_dispatch_2 (g_scm_truncate_remainder, x, y, SCM_ARG2, + s_scm_truncate_remainder); + } + else if (SCM_REALP (x)) + { + if (SCM_REALP (y) || SCM_I_INUMP (y) || + SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_inexact_truncate_remainder + (SCM_REAL_VALUE (x), scm_to_double (y)); + else + return scm_wta_dispatch_2 (g_scm_truncate_remainder, x, y, SCM_ARG2, + s_scm_truncate_remainder); + } + else if (SCM_FRACTIONP (x)) + { + if (SCM_REALP (y)) + return scm_i_inexact_truncate_remainder + (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_exact_rational_truncate_remainder (x, y); + else + return scm_wta_dispatch_2 (g_scm_truncate_remainder, x, y, SCM_ARG2, + s_scm_truncate_remainder); + } + else + return scm_wta_dispatch_2 (g_scm_truncate_remainder, x, y, SCM_ARG1, + s_scm_truncate_remainder); +} +#undef FUNC_NAME + +static SCM +scm_i_inexact_truncate_remainder (double x, double y) +{ + /* Although it would be more efficient to use fmod here, we can't + because it would in some cases produce results inconsistent with + scm_i_inexact_truncate_quotient, such that x != q * y + r (not even + close). In particular, when x is very close to a multiple of y, + then r might be either 0.0 or sgn(x)*|y|, but those two cases must + correspond to different choices of q. If quotient chooses one and + remainder chooses the other, it would be bad. */ + if (SCM_UNLIKELY (y == 0)) + scm_num_overflow (s_scm_truncate_remainder); /* or return a NaN? */ + else + return scm_from_double (x - y * trunc (x / y)); +} + +static SCM +scm_i_exact_rational_truncate_remainder (SCM x, SCM y) +{ + SCM xd = scm_denominator (x); + SCM yd = scm_denominator (y); + SCM r1 = scm_truncate_remainder (scm_product (scm_numerator (x), yd), + scm_product (scm_numerator (y), xd)); + return scm_divide (r1, scm_product (xd, yd)); +} + + +static void scm_i_inexact_truncate_divide (double x, double y, + SCM *qp, SCM *rp); +static void scm_i_exact_rational_truncate_divide (SCM x, SCM y, + SCM *qp, SCM *rp); + +SCM_PRIMITIVE_GENERIC (scm_i_truncate_divide, "truncate/", 2, 0, 0, + (SCM x, SCM y), + "Return the integer @var{q} and the real number @var{r}\n" + "such that @math{@var{x} = @var{q}*@var{y} + @var{r}}\n" + "and @math{@var{q} = truncate(@var{x} / @var{y})}.\n" + "@lisp\n" + "(truncate/ 123 10) @result{} 12 and 3\n" + "(truncate/ 123 -10) @result{} -12 and 3\n" + "(truncate/ -123 10) @result{} -12 and -3\n" + "(truncate/ -123 -10) @result{} 12 and -3\n" + "(truncate/ -123.2 -63.5) @result{} 1.0 and -59.7\n" + "(truncate/ 16/3 -10/7) @result{} -3 and 22/21\n" + "@end lisp") +#define FUNC_NAME s_scm_i_truncate_divide +{ + SCM q, r; + + scm_truncate_divide(x, y, &q, &r); + return scm_values (scm_list_2 (q, r)); +} +#undef FUNC_NAME + +#define s_scm_truncate_divide s_scm_i_truncate_divide +#define g_scm_truncate_divide g_scm_i_truncate_divide + +void +scm_truncate_divide (SCM x, SCM y, SCM *qp, SCM *rp) +{ + 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); + if (SCM_UNLIKELY (yy == 0)) + scm_num_overflow (s_scm_truncate_divide); + else + { + scm_t_inum qq = xx / yy; + scm_t_inum rr = xx % yy; + if (SCM_LIKELY (SCM_FIXABLE (qq))) + *qp = SCM_I_MAKINUM (qq); + else + *qp = scm_i_inum2big (qq); + *rp = SCM_I_MAKINUM (rr); + } + return; + } + else if (SCM_BIGP (y)) + { + if (SCM_UNLIKELY (xx == SCM_MOST_NEGATIVE_FIXNUM) + && SCM_UNLIKELY (mpz_cmp_ui (SCM_I_BIG_MPZ (y), + - SCM_MOST_NEGATIVE_FIXNUM) == 0)) + { + /* Special case: x == fixnum-min && y == abs (fixnum-min) */ + scm_remember_upto_here_1 (y); + *qp = SCM_I_MAKINUM (-1); + *rp = SCM_INUM0; + } + else + { + *qp = SCM_INUM0; + *rp = x; + } + return; + } + else if (SCM_REALP (y)) + return scm_i_inexact_truncate_divide (xx, SCM_REAL_VALUE (y), qp, rp); + else if (SCM_FRACTIONP (y)) + return scm_i_exact_rational_truncate_divide (x, y, qp, rp); + else + return two_valued_wta_dispatch_2 + (g_scm_truncate_divide, x, y, SCM_ARG2, + s_scm_truncate_divide, qp, rp); + } + else if (SCM_BIGP (x)) + { + if (SCM_LIKELY (SCM_I_INUMP (y))) + { + scm_t_inum yy = SCM_I_INUM (y); + if (SCM_UNLIKELY (yy == 0)) + scm_num_overflow (s_scm_truncate_divide); + else + { + SCM q = scm_i_mkbig (); + scm_t_inum rr; + if (yy > 0) + rr = mpz_tdiv_q_ui (SCM_I_BIG_MPZ (q), + SCM_I_BIG_MPZ (x), yy); + else + { + rr = mpz_tdiv_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)); + } + rr *= mpz_sgn (SCM_I_BIG_MPZ (x)); + scm_remember_upto_here_1 (x); + *qp = scm_i_normbig (q); + *rp = SCM_I_MAKINUM (rr); + } + return; + } + else if (SCM_BIGP (y)) + { + SCM q = scm_i_mkbig (); + SCM r = scm_i_mkbig (); + mpz_tdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r), + SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y)); + scm_remember_upto_here_2 (x, y); + *qp = scm_i_normbig (q); + *rp = scm_i_normbig (r); + } + else if (SCM_REALP (y)) + return scm_i_inexact_truncate_divide + (scm_i_big2dbl (x), SCM_REAL_VALUE (y), qp, rp); + else if (SCM_FRACTIONP (y)) + return scm_i_exact_rational_truncate_divide (x, y, qp, rp); + else + return two_valued_wta_dispatch_2 + (g_scm_truncate_divide, x, y, SCM_ARG2, + s_scm_truncate_divide, qp, rp); + } + else if (SCM_REALP (x)) + { + if (SCM_REALP (y) || SCM_I_INUMP (y) || + SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_inexact_truncate_divide + (SCM_REAL_VALUE (x), scm_to_double (y), qp, rp); + else + return two_valued_wta_dispatch_2 + (g_scm_truncate_divide, x, y, SCM_ARG2, + s_scm_truncate_divide, qp, rp); + } + else if (SCM_FRACTIONP (x)) + { + if (SCM_REALP (y)) + return scm_i_inexact_truncate_divide + (scm_i_fraction2double (x), SCM_REAL_VALUE (y), qp, rp); + else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_exact_rational_truncate_divide (x, y, qp, rp); + else + return two_valued_wta_dispatch_2 + (g_scm_truncate_divide, x, y, SCM_ARG2, + s_scm_truncate_divide, qp, rp); + } + else + return two_valued_wta_dispatch_2 (g_scm_truncate_divide, x, y, SCM_ARG1, + s_scm_truncate_divide, qp, rp); +} + +static void +scm_i_inexact_truncate_divide (double x, double y, SCM *qp, SCM *rp) +{ + if (SCM_UNLIKELY (y == 0)) + scm_num_overflow (s_scm_truncate_divide); /* or return a NaN? */ + else + { + double q = trunc (x / y); + double r = x - q * y; + *qp = scm_from_double (q); + *rp = scm_from_double (r); + } +} + +static void +scm_i_exact_rational_truncate_divide (SCM x, SCM y, SCM *qp, SCM *rp) +{ + SCM r1; + SCM xd = scm_denominator (x); + SCM yd = scm_denominator (y); + + scm_truncate_divide (scm_product (scm_numerator (x), yd), + scm_product (scm_numerator (y), xd), + qp, &r1); + *rp = scm_divide (r1, scm_product (xd, yd)); +} + +static SCM scm_i_inexact_centered_quotient (double x, double y); +static SCM scm_i_bigint_centered_quotient (SCM x, SCM y); +static SCM scm_i_exact_rational_centered_quotient (SCM x, SCM y); + +SCM_PRIMITIVE_GENERIC (scm_centered_quotient, "centered-quotient", 2, 0, 0, + (SCM x, SCM y), + "Return the integer @var{q} such that\n" + "@math{@var{x} = @var{q}*@var{y} + @var{r}} where\n" + "@math{-abs(@var{y}/2) <= @var{r} < abs(@var{y}/2)}.\n" + "@lisp\n" + "(centered-quotient 123 10) @result{} 12\n" + "(centered-quotient 123 -10) @result{} -12\n" + "(centered-quotient -123 10) @result{} -12\n" + "(centered-quotient -123 -10) @result{} 12\n" + "(centered-quotient -123.2 -63.5) @result{} 2.0\n" + "(centered-quotient 16/3 -10/7) @result{} -4\n" + "@end lisp") +#define FUNC_NAME s_scm_centered_quotient +{ + 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); + if (SCM_UNLIKELY (yy == 0)) + scm_num_overflow (s_scm_centered_quotient); + else + { + scm_t_inum qq = xx / yy; + scm_t_inum rr = xx % yy; + if (SCM_LIKELY (xx > 0)) + { + if (SCM_LIKELY (yy > 0)) + { + if (rr >= (yy + 1) / 2) + qq++; + } + else + { + if (rr >= (1 - yy) / 2) + qq--; + } + } + else + { + if (SCM_LIKELY (yy > 0)) + { + if (rr < -yy / 2) + qq--; + } + else + { + if (rr < yy / 2) + 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 (xx), y); + } + else if (SCM_REALP (y)) + return scm_i_inexact_centered_quotient (xx, SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_exact_rational_centered_quotient (x, y); + else + return scm_wta_dispatch_2 (g_scm_centered_quotient, x, y, SCM_ARG2, + s_scm_centered_quotient); + } + else if (SCM_BIGP (x)) + { + if (SCM_LIKELY (SCM_I_INUMP (y))) + { + 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 (); + scm_t_inum rr; + /* Arrange for rr to initially be non-positive, + because that simplifies the test to see + if it is within the needed bounds. */ + if (yy > 0) + { + rr = - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q), + SCM_I_BIG_MPZ (x), yy); + scm_remember_upto_here_1 (x); + if (rr < -yy / 2) + mpz_sub_ui (SCM_I_BIG_MPZ (q), + SCM_I_BIG_MPZ (q), 1); + } + else + { + rr = - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q), + SCM_I_BIG_MPZ (x), -yy); + scm_remember_upto_here_1 (x); + mpz_neg (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (q)); + if (rr < yy / 2) + mpz_add_ui (SCM_I_BIG_MPZ (q), + SCM_I_BIG_MPZ (q), 1); + } + return scm_i_normbig (q); + } + } + else if (SCM_BIGP (y)) + return scm_i_bigint_centered_quotient (x, y); + else if (SCM_REALP (y)) + return scm_i_inexact_centered_quotient + (scm_i_big2dbl (x), SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_exact_rational_centered_quotient (x, y); + else + return scm_wta_dispatch_2 (g_scm_centered_quotient, x, y, SCM_ARG2, + s_scm_centered_quotient); + } + else if (SCM_REALP (x)) + { + if (SCM_REALP (y) || SCM_I_INUMP (y) || + SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_inexact_centered_quotient + (SCM_REAL_VALUE (x), scm_to_double (y)); + else + return scm_wta_dispatch_2 (g_scm_centered_quotient, x, y, SCM_ARG2, + s_scm_centered_quotient); + } + else if (SCM_FRACTIONP (x)) + { + if (SCM_REALP (y)) + return scm_i_inexact_centered_quotient + (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_exact_rational_centered_quotient (x, y); + else + return scm_wta_dispatch_2 (g_scm_centered_quotient, x, y, SCM_ARG2, + s_scm_centered_quotient); + } + else + return scm_wta_dispatch_2 (g_scm_centered_quotient, x, y, SCM_ARG1, + s_scm_centered_quotient); +} +#undef FUNC_NAME + +static SCM +scm_i_inexact_centered_quotient (double x, double y) +{ + if (SCM_LIKELY (y > 0)) + return scm_from_double (floor (x/y + 0.5)); + else if (SCM_LIKELY (y < 0)) + return scm_from_double (ceil (x/y - 0.5)); + else if (y == 0) + scm_num_overflow (s_scm_centered_quotient); /* or return a NaN? */ + else + return scm_nan (); +} + +/* Assumes that both x and y are bigints, though + x might be able to fit into a fixnum. */ +static SCM +scm_i_bigint_centered_quotient (SCM x, SCM y) +{ + SCM q, r, min_r; + + /* Note that x might be small enough to fit into a + fixnum, so we must not let it escape into the wild */ + q = scm_i_mkbig (); + r = scm_i_mkbig (); + + /* min_r will eventually become -abs(y)/2 */ + min_r = scm_i_mkbig (); + mpz_tdiv_q_2exp (SCM_I_BIG_MPZ (min_r), + SCM_I_BIG_MPZ (y), 1); + + /* Arrange for rr to initially be non-positive, + because that simplifies the test to see + if it is within the needed bounds. */ + if (mpz_sgn (SCM_I_BIG_MPZ (y)) > 0) + { + mpz_cdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r), + SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y)); + scm_remember_upto_here_2 (x, y); + mpz_neg (SCM_I_BIG_MPZ (min_r), SCM_I_BIG_MPZ (min_r)); + if (mpz_cmp (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (min_r)) < 0) + mpz_sub_ui (SCM_I_BIG_MPZ (q), + SCM_I_BIG_MPZ (q), 1); + } + else + { + mpz_fdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r), + SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y)); + scm_remember_upto_here_2 (x, y); + if (mpz_cmp (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (min_r)) < 0) + mpz_add_ui (SCM_I_BIG_MPZ (q), + SCM_I_BIG_MPZ (q), 1); + } + scm_remember_upto_here_2 (r, min_r); + return scm_i_normbig (q); +} + +static SCM +scm_i_exact_rational_centered_quotient (SCM x, SCM y) +{ + return scm_centered_quotient + (scm_product (scm_numerator (x), scm_denominator (y)), + scm_product (scm_numerator (y), scm_denominator (x))); +} + +static SCM scm_i_inexact_centered_remainder (double x, double y); +static SCM scm_i_bigint_centered_remainder (SCM x, SCM y); +static SCM scm_i_exact_rational_centered_remainder (SCM x, SCM y); + +SCM_PRIMITIVE_GENERIC (scm_centered_remainder, "centered-remainder", 2, 0, 0, + (SCM x, SCM y), + "Return the real number @var{r} such that\n" + "@math{-abs(@var{y}/2) <= @var{r} < abs(@var{y}/2)}\n" + "and @math{@var{x} = @var{q}*@var{y} + @var{r}}\n" + "for some integer @var{q}.\n" + "@lisp\n" + "(centered-remainder 123 10) @result{} 3\n" + "(centered-remainder 123 -10) @result{} 3\n" + "(centered-remainder -123 10) @result{} -3\n" + "(centered-remainder -123 -10) @result{} -3\n" + "(centered-remainder -123.2 -63.5) @result{} 3.8\n" + "(centered-remainder 16/3 -10/7) @result{} -8/21\n" + "@end lisp") +#define FUNC_NAME s_scm_centered_remainder +{ + 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); + if (SCM_UNLIKELY (yy == 0)) + scm_num_overflow (s_scm_centered_remainder); + else + { + scm_t_inum rr = xx % yy; + if (SCM_LIKELY (xx > 0)) + { + if (SCM_LIKELY (yy > 0)) + { + if (rr >= (yy + 1) / 2) + rr -= yy; + } + else + { + if (rr >= (1 - yy) / 2) + rr += yy; + } + } + else + { + if (SCM_LIKELY (yy > 0)) + { + if (rr < -yy / 2) + rr += yy; + } + else + { + if (rr < yy / 2) + rr -= yy; + } + } + return 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_remainder */ + return scm_i_bigint_centered_remainder (scm_i_long2big (xx), y); + } + else if (SCM_REALP (y)) + return scm_i_inexact_centered_remainder (xx, SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_exact_rational_centered_remainder (x, y); + else + return scm_wta_dispatch_2 (g_scm_centered_remainder, x, y, SCM_ARG2, + s_scm_centered_remainder); + } + else if (SCM_BIGP (x)) + { + if (SCM_LIKELY (SCM_I_INUMP (y))) + { + scm_t_inum yy = SCM_I_INUM (y); + if (SCM_UNLIKELY (yy == 0)) + scm_num_overflow (s_scm_centered_remainder); + else + { + scm_t_inum rr; + /* Arrange for rr to initially be non-positive, + because that simplifies the test to see + if it is within the needed bounds. */ + if (yy > 0) + { + rr = - mpz_cdiv_ui (SCM_I_BIG_MPZ (x), yy); + scm_remember_upto_here_1 (x); + if (rr < -yy / 2) + rr += yy; + } + else + { + rr = - mpz_cdiv_ui (SCM_I_BIG_MPZ (x), -yy); + scm_remember_upto_here_1 (x); + if (rr < yy / 2) + rr -= yy; + } + return SCM_I_MAKINUM (rr); + } + } + else if (SCM_BIGP (y)) + return scm_i_bigint_centered_remainder (x, y); + else if (SCM_REALP (y)) + return scm_i_inexact_centered_remainder + (scm_i_big2dbl (x), SCM_REAL_VALUE (y)); + else if (SCM_FRACTIONP (y)) + return scm_i_exact_rational_centered_remainder (x, y); + else + return scm_wta_dispatch_2 (g_scm_centered_remainder, x, y, SCM_ARG2, + s_scm_centered_remainder); + } + else if (SCM_REALP (x)) + { + if (SCM_REALP (y) || SCM_I_INUMP (y) || + SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_inexact_centered_remainder + (SCM_REAL_VALUE (x), scm_to_double (y)); + else + return scm_wta_dispatch_2 (g_scm_centered_remainder, x, y, SCM_ARG2, + s_scm_centered_remainder); + } + else if (SCM_FRACTIONP (x)) + { + if (SCM_REALP (y)) + return scm_i_inexact_centered_remainder + (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_exact_rational_centered_remainder (x, y); + else + return scm_wta_dispatch_2 (g_scm_centered_remainder, x, y, SCM_ARG2, + s_scm_centered_remainder); + } + else + return scm_wta_dispatch_2 (g_scm_centered_remainder, x, y, SCM_ARG1, + s_scm_centered_remainder); +} +#undef FUNC_NAME + +static SCM +scm_i_inexact_centered_remainder (double x, double y) +{ + double q; + + /* Although it would be more efficient to use fmod here, we can't + because it would in some cases produce results inconsistent with + scm_i_inexact_centered_quotient, such that x != r + q * y (not even + close). In particular, when x-y/2 is very close to a multiple of + y, then r might be either -abs(y/2) or abs(y/2)-epsilon, but those + two cases must correspond to different choices of q. If quotient + chooses one and remainder chooses the other, it would be bad. */ + if (SCM_LIKELY (y > 0)) + q = floor (x/y + 0.5); + else if (SCM_LIKELY (y < 0)) + q = ceil (x/y - 0.5); + else if (y == 0) + scm_num_overflow (s_scm_centered_remainder); /* or return a NaN? */ + else + return scm_nan (); + return scm_from_double (x - q * y); +} + +/* Assumes that both x and y are bigints, though + x might be able to fit into a fixnum. */ +static SCM +scm_i_bigint_centered_remainder (SCM x, SCM y) +{ + SCM r, min_r; + + /* Note that x might be small enough to fit into a + fixnum, so we must not let it escape into the wild */ + r = scm_i_mkbig (); + + /* min_r will eventually become -abs(y)/2 */ + min_r = scm_i_mkbig (); + mpz_tdiv_q_2exp (SCM_I_BIG_MPZ (min_r), + SCM_I_BIG_MPZ (y), 1); + + /* Arrange for rr to initially be non-positive, + because that simplifies the test to see + if it is within the needed bounds. */ + if (mpz_sgn (SCM_I_BIG_MPZ (y)) > 0) + { + mpz_cdiv_r (SCM_I_BIG_MPZ (r), + SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y)); + mpz_neg (SCM_I_BIG_MPZ (min_r), SCM_I_BIG_MPZ (min_r)); + if (mpz_cmp (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (min_r)) < 0) + mpz_add (SCM_I_BIG_MPZ (r), + SCM_I_BIG_MPZ (r), + SCM_I_BIG_MPZ (y)); + } + else + { + mpz_fdiv_r (SCM_I_BIG_MPZ (r), + SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y)); + if (mpz_cmp (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (min_r)) < 0) + mpz_sub (SCM_I_BIG_MPZ (r), + SCM_I_BIG_MPZ (r), + SCM_I_BIG_MPZ (y)); + } + scm_remember_upto_here_2 (x, y); + return scm_i_normbig (r); +} + +static SCM +scm_i_exact_rational_centered_remainder (SCM x, SCM y) +{ + SCM xd = scm_denominator (x); + SCM yd = scm_denominator (y); + SCM r1 = scm_centered_remainder (scm_product (scm_numerator (x), yd), + scm_product (scm_numerator (y), xd)); + return scm_divide (r1, scm_product (xd, yd)); +} + + +static void scm_i_inexact_centered_divide (double x, double y, + SCM *qp, SCM *rp); +static void scm_i_bigint_centered_divide (SCM x, SCM y, SCM *qp, SCM *rp); +static void scm_i_exact_rational_centered_divide (SCM x, SCM y, + SCM *qp, SCM *rp); + +SCM_PRIMITIVE_GENERIC (scm_i_centered_divide, "centered/", 2, 0, 0, + (SCM x, SCM y), + "Return the integer @var{q} and the real number @var{r}\n" + "such that @math{@var{x} = @var{q}*@var{y} + @var{r}}\n" + "and @math{-abs(@var{y}/2) <= @var{r} < abs(@var{y}/2)}.\n" + "@lisp\n" + "(centered/ 123 10) @result{} 12 and 3\n" + "(centered/ 123 -10) @result{} -12 and 3\n" + "(centered/ -123 10) @result{} -12 and -3\n" + "(centered/ -123 -10) @result{} 12 and -3\n" + "(centered/ -123.2 -63.5) @result{} 2.0 and 3.8\n" + "(centered/ 16/3 -10/7) @result{} -4 and -8/21\n" + "@end lisp") +#define FUNC_NAME s_scm_i_centered_divide +{ + SCM q, r; + + scm_centered_divide(x, y, &q, &r); + return scm_values (scm_list_2 (q, r)); +} +#undef FUNC_NAME + +#define s_scm_centered_divide s_scm_i_centered_divide +#define g_scm_centered_divide g_scm_i_centered_divide + +void +scm_centered_divide (SCM x, SCM y, SCM *qp, SCM *rp) +{ + 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); + if (SCM_UNLIKELY (yy == 0)) + scm_num_overflow (s_scm_centered_divide); + else + { + scm_t_inum qq = xx / yy; + scm_t_inum rr = xx % yy; + if (SCM_LIKELY (xx > 0)) + { + if (SCM_LIKELY (yy > 0)) + { + if (rr >= (yy + 1) / 2) + { qq++; rr -= yy; } + } + else + { + if (rr >= (1 - yy) / 2) + { qq--; rr += yy; } + } + } + else + { + if (SCM_LIKELY (yy > 0)) + { + if (rr < -yy / 2) + { qq--; rr += yy; } + } + else + { + if (rr < yy / 2) + { qq++; rr -= yy; } + } + } + if (SCM_LIKELY (SCM_FIXABLE (qq))) + *qp = SCM_I_MAKINUM (qq); + else + *qp = scm_i_inum2big (qq); + *rp = SCM_I_MAKINUM (rr); + } + return; + } + 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 (xx), y, qp, rp); + } + else if (SCM_REALP (y)) + return scm_i_inexact_centered_divide (xx, SCM_REAL_VALUE (y), qp, rp); + else if (SCM_FRACTIONP (y)) + return scm_i_exact_rational_centered_divide (x, y, qp, rp); + else + return two_valued_wta_dispatch_2 + (g_scm_centered_divide, x, y, SCM_ARG2, + s_scm_centered_divide, qp, rp); + } + else if (SCM_BIGP (x)) + { + if (SCM_LIKELY (SCM_I_INUMP (y))) + { + scm_t_inum yy = SCM_I_INUM (y); + if (SCM_UNLIKELY (yy == 0)) + scm_num_overflow (s_scm_centered_divide); + else + { + SCM q = scm_i_mkbig (); + scm_t_inum rr; + /* Arrange for rr to initially be non-positive, + because that simplifies the test to see + if it is within the needed bounds. */ + if (yy > 0) + { + rr = - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q), + SCM_I_BIG_MPZ (x), yy); + scm_remember_upto_here_1 (x); + if (rr < -yy / 2) + { + mpz_sub_ui (SCM_I_BIG_MPZ (q), + SCM_I_BIG_MPZ (q), 1); + rr += yy; + } + } + else + { + rr = - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q), + SCM_I_BIG_MPZ (x), -yy); + scm_remember_upto_here_1 (x); + mpz_neg (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (q)); + if (rr < yy / 2) + { + mpz_add_ui (SCM_I_BIG_MPZ (q), + SCM_I_BIG_MPZ (q), 1); + rr -= yy; + } + } + *qp = scm_i_normbig (q); + *rp = SCM_I_MAKINUM (rr); + } + return; + } + else if (SCM_BIGP (y)) + return scm_i_bigint_centered_divide (x, y, qp, rp); + else if (SCM_REALP (y)) + return scm_i_inexact_centered_divide + (scm_i_big2dbl (x), SCM_REAL_VALUE (y), qp, rp); + else if (SCM_FRACTIONP (y)) + return scm_i_exact_rational_centered_divide (x, y, qp, rp); + else + return two_valued_wta_dispatch_2 + (g_scm_centered_divide, x, y, SCM_ARG2, + s_scm_centered_divide, qp, rp); + } + else if (SCM_REALP (x)) + { + if (SCM_REALP (y) || SCM_I_INUMP (y) || + SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_inexact_centered_divide + (SCM_REAL_VALUE (x), scm_to_double (y), qp, rp); + else + return two_valued_wta_dispatch_2 + (g_scm_centered_divide, x, y, SCM_ARG2, + s_scm_centered_divide, qp, rp); + } + else if (SCM_FRACTIONP (x)) + { + if (SCM_REALP (y)) + return scm_i_inexact_centered_divide + (scm_i_fraction2double (x), SCM_REAL_VALUE (y), qp, rp); + else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_exact_rational_centered_divide (x, y, qp, rp); + else + return two_valued_wta_dispatch_2 + (g_scm_centered_divide, x, y, SCM_ARG2, + s_scm_centered_divide, qp, rp); + } + else + return two_valued_wta_dispatch_2 (g_scm_centered_divide, x, y, SCM_ARG1, + s_scm_centered_divide, qp, rp); +} + +static void +scm_i_inexact_centered_divide (double x, double y, SCM *qp, SCM *rp) +{ + double q, r; + + if (SCM_LIKELY (y > 0)) + q = floor (x/y + 0.5); + else if (SCM_LIKELY (y < 0)) + q = ceil (x/y - 0.5); + else if (y == 0) + scm_num_overflow (s_scm_centered_divide); /* or return a NaN? */ + else + q = guile_NaN; + r = x - q * y; + *qp = scm_from_double (q); + *rp = scm_from_double (r); +} + +/* Assumes that both x and y are bigints, though + x might be able to fit into a fixnum. */ +static void +scm_i_bigint_centered_divide (SCM x, SCM y, SCM *qp, SCM *rp) +{ + SCM q, r, min_r; + + /* Note that x might be small enough to fit into a + fixnum, so we must not let it escape into the wild */ + q = scm_i_mkbig (); + r = scm_i_mkbig (); + + /* min_r will eventually become -abs(y/2) */ + min_r = scm_i_mkbig (); + mpz_tdiv_q_2exp (SCM_I_BIG_MPZ (min_r), + SCM_I_BIG_MPZ (y), 1); + + /* Arrange for rr to initially be non-positive, + because that simplifies the test to see + if it is within the needed bounds. */ + if (mpz_sgn (SCM_I_BIG_MPZ (y)) > 0) + { + mpz_cdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r), + SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y)); + mpz_neg (SCM_I_BIG_MPZ (min_r), SCM_I_BIG_MPZ (min_r)); + if (mpz_cmp (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (min_r)) < 0) + { + mpz_sub_ui (SCM_I_BIG_MPZ (q), + SCM_I_BIG_MPZ (q), 1); + mpz_add (SCM_I_BIG_MPZ (r), + SCM_I_BIG_MPZ (r), + SCM_I_BIG_MPZ (y)); + } + } else - q = guile_NaN; - r = x - q * y; - return scm_values (scm_list_2 (scm_from_double (q), - scm_from_double (r))); + { + mpz_fdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r), + SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y)); + if (mpz_cmp (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (min_r)) < 0) + { + mpz_add_ui (SCM_I_BIG_MPZ (q), + SCM_I_BIG_MPZ (q), 1); + mpz_sub (SCM_I_BIG_MPZ (r), + SCM_I_BIG_MPZ (r), + SCM_I_BIG_MPZ (y)); + } + } + scm_remember_upto_here_2 (x, y); + *qp = scm_i_normbig (q); + *rp = scm_i_normbig (r); } -/* Compute exact euclidean quotient and remainder the slow way. - We use this only if both arguments are exact, - and at least one of them is a fraction */ -static SCM -scm_i_slow_exact_euclidean_divide (SCM x, SCM y) +static void +scm_i_exact_rational_centered_divide (SCM x, SCM y, SCM *qp, SCM *rp) { - SCM q, r; + SCM r1; + SCM xd = scm_denominator (x); + SCM yd = scm_denominator (y); - if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x))) - SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide, x, y, SCM_ARG1, - s_scm_euclidean_divide); - else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y))) - SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide, x, y, SCM_ARG2, - s_scm_euclidean_divide); - else if (scm_is_true (scm_positive_p (y))) - q = scm_floor (scm_divide (x, y)); - else if (scm_is_true (scm_negative_p (y))) - q = scm_ceiling (scm_divide (x, y)); - else - scm_num_overflow (s_scm_euclidean_divide); - r = scm_difference (x, scm_product (q, y)); - return scm_values (scm_list_2 (q, r)); + scm_centered_divide (scm_product (scm_numerator (x), yd), + scm_product (scm_numerator (y), xd), + qp, &r1); + *rp = scm_divide (r1, scm_product (xd, yd)); } -static SCM scm_i_inexact_centered_quotient (double x, double y); -static SCM scm_i_bigint_centered_quotient (SCM x, SCM y); -static SCM scm_i_slow_exact_centered_quotient (SCM x, SCM y); +static SCM scm_i_inexact_round_quotient (double x, double y); +static SCM scm_i_bigint_round_quotient (SCM x, SCM y); +static SCM scm_i_exact_rational_round_quotient (SCM x, SCM y); -SCM_PRIMITIVE_GENERIC (scm_centered_quotient, "centered-quotient", 2, 0, 0, +SCM_PRIMITIVE_GENERIC (scm_round_quotient, "round-quotient", 2, 0, 0, (SCM x, SCM y), - "Return the integer @var{q} such that\n" - "@math{@var{x} = @var{q}*@var{y} + @var{r}} where\n" - "@math{-abs(@var{y}/2) <= @var{r} < abs(@var{y}/2)}.\n" + "Return @math{@var{x} / @var{y}} to the nearest integer,\n" + "with ties going to the nearest even integer.\n" "@lisp\n" - "(centered-quotient 123 10) @result{} 12\n" - "(centered-quotient 123 -10) @result{} -12\n" - "(centered-quotient -123 10) @result{} -12\n" - "(centered-quotient -123 -10) @result{} 12\n" - "(centered-quotient -123.2 -63.5) @result{} 2.0\n" - "(centered-quotient 16/3 -10/7) @result{} -4\n" + "(round-quotient 123 10) @result{} 12\n" + "(round-quotient 123 -10) @result{} -12\n" + "(round-quotient -123 10) @result{} -12\n" + "(round-quotient -123 -10) @result{} 12\n" + "(round-quotient 125 10) @result{} 12\n" + "(round-quotient 127 10) @result{} 13\n" + "(round-quotient 135 10) @result{} 14\n" + "(round-quotient -123.2 -63.5) @result{} 2.0\n" + "(round-quotient 16/3 -10/7) @result{} -4\n" "@end lisp") -#define FUNC_NAME s_scm_centered_quotient +#define FUNC_NAME s_scm_round_quotient { 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); if (SCM_UNLIKELY (yy == 0)) - scm_num_overflow (s_scm_centered_quotient); + scm_num_overflow (s_scm_round_quotient); else { - scm_t_inum xx = SCM_I_INUM (x); scm_t_inum qq = xx / yy; - scm_t_inum rr = xx - qq * yy; - if (SCM_LIKELY (xx > 0)) + scm_t_inum rr = xx % yy; + scm_t_inum ay = yy; + scm_t_inum r2 = 2 * rr; + + if (SCM_LIKELY (yy < 0)) { - if (SCM_LIKELY (yy > 0)) - { - if (rr >= (yy + 1) / 2) - qq++; - } - else - { - if (rr >= (1 - yy) / 2) - qq--; - } + ay = -ay; + r2 = -r2; + } + + if (qq & 1L) + { + if (r2 >= ay) + qq++; + else if (r2 <= -ay) + qq--; } else { - if (SCM_LIKELY (yy > 0)) - { - if (rr < -yy / 2) - qq--; - } - else - { - if (rr < yy / 2) - qq++; - } + if (r2 > ay) + qq++; + else if (r2 < -ay) + 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); + can fit in a fixnum) to scm_i_bigint_round_quotient */ + return scm_i_bigint_round_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_round_quotient (xx, SCM_REAL_VALUE (y)); else if (SCM_FRACTIONP (y)) - return scm_i_slow_exact_centered_quotient (x, y); + return scm_i_exact_rational_round_quotient (x, y); else - SCM_WTA_DISPATCH_2 (g_scm_centered_quotient, x, y, SCM_ARG2, - s_scm_centered_quotient); + return scm_wta_dispatch_2 (g_scm_round_quotient, x, y, SCM_ARG2, + s_scm_round_quotient); } else if (SCM_BIGP (x)) { @@ -1668,205 +3258,185 @@ 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); + scm_num_overflow (s_scm_round_quotient); + else if (SCM_UNLIKELY (yy == 1)) + return x; else { SCM q = scm_i_mkbig (); scm_t_inum rr; - /* Arrange for rr to initially be non-positive, - because that simplifies the test to see - if it is within the needed bounds. */ + int needs_adjustment; + if (yy > 0) { - rr = - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q), - SCM_I_BIG_MPZ (x), yy); - scm_remember_upto_here_1 (x); - if (rr < -yy / 2) - mpz_sub_ui (SCM_I_BIG_MPZ (q), - SCM_I_BIG_MPZ (q), 1); + rr = mpz_fdiv_q_ui (SCM_I_BIG_MPZ (q), + SCM_I_BIG_MPZ (x), yy); + if (mpz_odd_p (SCM_I_BIG_MPZ (q))) + needs_adjustment = (2*rr >= yy); + else + needs_adjustment = (2*rr > yy); } else { rr = - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (x), -yy); - scm_remember_upto_here_1 (x); mpz_neg (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (q)); - if (rr < yy / 2) - mpz_add_ui (SCM_I_BIG_MPZ (q), - SCM_I_BIG_MPZ (q), 1); + if (mpz_odd_p (SCM_I_BIG_MPZ (q))) + needs_adjustment = (2*rr <= yy); + else + needs_adjustment = (2*rr < yy); } + scm_remember_upto_here_1 (x); + if (needs_adjustment) + mpz_add_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (q), 1); return scm_i_normbig (q); } } else if (SCM_BIGP (y)) - return scm_i_bigint_centered_quotient (x, y); + return scm_i_bigint_round_quotient (x, y); else if (SCM_REALP (y)) - return scm_i_inexact_centered_quotient + return scm_i_inexact_round_quotient (scm_i_big2dbl (x), SCM_REAL_VALUE (y)); else if (SCM_FRACTIONP (y)) - return scm_i_slow_exact_centered_quotient (x, y); + return scm_i_exact_rational_round_quotient (x, y); else - SCM_WTA_DISPATCH_2 (g_scm_centered_quotient, x, y, SCM_ARG2, - s_scm_centered_quotient); + return scm_wta_dispatch_2 (g_scm_round_quotient, x, y, SCM_ARG2, + s_scm_round_quotient); } else if (SCM_REALP (x)) { if (SCM_REALP (y) || SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) - return scm_i_inexact_centered_quotient + return scm_i_inexact_round_quotient (SCM_REAL_VALUE (x), scm_to_double (y)); else - SCM_WTA_DISPATCH_2 (g_scm_centered_quotient, x, y, SCM_ARG2, - s_scm_centered_quotient); + return scm_wta_dispatch_2 (g_scm_round_quotient, x, y, SCM_ARG2, + s_scm_round_quotient); } else if (SCM_FRACTIONP (x)) { if (SCM_REALP (y)) - return scm_i_inexact_centered_quotient + return scm_i_inexact_round_quotient (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_exact_rational_round_quotient (x, y); else - return scm_i_slow_exact_centered_quotient (x, y); + return scm_wta_dispatch_2 (g_scm_round_quotient, x, y, SCM_ARG2, + s_scm_round_quotient); } else - SCM_WTA_DISPATCH_2 (g_scm_centered_quotient, x, y, SCM_ARG1, - s_scm_centered_quotient); + return scm_wta_dispatch_2 (g_scm_round_quotient, x, y, SCM_ARG1, + s_scm_round_quotient); } #undef FUNC_NAME static SCM -scm_i_inexact_centered_quotient (double x, double y) +scm_i_inexact_round_quotient (double x, double y) { - if (SCM_LIKELY (y > 0)) - return scm_from_double (floor (x/y + 0.5)); - else if (SCM_LIKELY (y < 0)) - return scm_from_double (ceil (x/y - 0.5)); - else if (y == 0) - scm_num_overflow (s_scm_centered_quotient); /* or return a NaN? */ + if (SCM_UNLIKELY (y == 0)) + scm_num_overflow (s_scm_round_quotient); /* or return a NaN? */ else - return scm_nan (); + return scm_from_double (scm_c_round (x / y)); } /* Assumes that both x and y are bigints, though x might be able to fit into a fixnum. */ static SCM -scm_i_bigint_centered_quotient (SCM x, SCM y) +scm_i_bigint_round_quotient (SCM x, SCM y) { - SCM q, r, min_r; + SCM q, r, r2; + int cmp, needs_adjustment; /* Note that x might be small enough to fit into a fixnum, so we must not let it escape into the wild */ q = scm_i_mkbig (); r = scm_i_mkbig (); + r2 = scm_i_mkbig (); - /* min_r will eventually become -abs(y)/2 */ - min_r = scm_i_mkbig (); - mpz_tdiv_q_2exp (SCM_I_BIG_MPZ (min_r), - SCM_I_BIG_MPZ (y), 1); + mpz_fdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r), + SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y)); + mpz_mul_2exp (SCM_I_BIG_MPZ (r2), SCM_I_BIG_MPZ (r), 1); /* r2 = 2*r */ + scm_remember_upto_here_2 (x, r); - /* Arrange for rr to initially be non-positive, - because that simplifies the test to see - if it is within the needed bounds. */ - if (mpz_sgn (SCM_I_BIG_MPZ (y)) > 0) - { - mpz_cdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r), - SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y)); - scm_remember_upto_here_2 (x, y); - mpz_neg (SCM_I_BIG_MPZ (min_r), SCM_I_BIG_MPZ (min_r)); - if (mpz_cmp (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (min_r)) < 0) - mpz_sub_ui (SCM_I_BIG_MPZ (q), - SCM_I_BIG_MPZ (q), 1); - } + cmp = mpz_cmpabs (SCM_I_BIG_MPZ (r2), SCM_I_BIG_MPZ (y)); + if (mpz_odd_p (SCM_I_BIG_MPZ (q))) + needs_adjustment = (cmp >= 0); else - { - mpz_fdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r), - SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y)); - scm_remember_upto_here_2 (x, y); - if (mpz_cmp (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (min_r)) < 0) - mpz_add_ui (SCM_I_BIG_MPZ (q), - SCM_I_BIG_MPZ (q), 1); - } - scm_remember_upto_here_2 (r, min_r); + needs_adjustment = (cmp > 0); + scm_remember_upto_here_2 (r2, y); + + if (needs_adjustment) + mpz_add_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (q), 1); + return scm_i_normbig (q); } -/* Compute exact centered quotient the slow way. - We use this only if both arguments are exact, - and at least one of them is a fraction */ static SCM -scm_i_slow_exact_centered_quotient (SCM x, SCM y) +scm_i_exact_rational_round_quotient (SCM x, SCM y) { - if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x))) - SCM_WTA_DISPATCH_2 (g_scm_centered_quotient, x, y, SCM_ARG1, - s_scm_centered_quotient); - else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y))) - SCM_WTA_DISPATCH_2 (g_scm_centered_quotient, x, y, SCM_ARG2, - s_scm_centered_quotient); - else if (scm_is_true (scm_positive_p (y))) - return scm_floor (scm_sum (scm_divide (x, y), - exactly_one_half)); - else if (scm_is_true (scm_negative_p (y))) - return scm_ceiling (scm_difference (scm_divide (x, y), - exactly_one_half)); - else - scm_num_overflow (s_scm_centered_quotient); + return scm_round_quotient + (scm_product (scm_numerator (x), scm_denominator (y)), + scm_product (scm_numerator (y), scm_denominator (x))); } -static SCM scm_i_inexact_centered_remainder (double x, double y); -static SCM scm_i_bigint_centered_remainder (SCM x, SCM y); -static SCM scm_i_slow_exact_centered_remainder (SCM x, SCM y); +static SCM scm_i_inexact_round_remainder (double x, double y); +static SCM scm_i_bigint_round_remainder (SCM x, SCM y); +static SCM scm_i_exact_rational_round_remainder (SCM x, SCM y); -SCM_PRIMITIVE_GENERIC (scm_centered_remainder, "centered-remainder", 2, 0, 0, +SCM_PRIMITIVE_GENERIC (scm_round_remainder, "round-remainder", 2, 0, 0, (SCM x, SCM y), "Return the real number @var{r} such that\n" - "@math{-abs(@var{y}/2) <= @var{r} < abs(@var{y}/2)}\n" - "and @math{@var{x} = @var{q}*@var{y} + @var{r}}\n" - "for some integer @var{q}.\n" + "@math{@var{x} = @var{q}*@var{y} + @var{r}}, where\n" + "@var{q} is @math{@var{x} / @var{y}} rounded to the\n" + "nearest integer, with ties going to the nearest\n" + "even integer.\n" "@lisp\n" - "(centered-remainder 123 10) @result{} 3\n" - "(centered-remainder 123 -10) @result{} 3\n" - "(centered-remainder -123 10) @result{} -3\n" - "(centered-remainder -123 -10) @result{} -3\n" - "(centered-remainder -123.2 -63.5) @result{} 3.8\n" - "(centered-remainder 16/3 -10/7) @result{} -8/21\n" + "(round-remainder 123 10) @result{} 3\n" + "(round-remainder 123 -10) @result{} 3\n" + "(round-remainder -123 10) @result{} -3\n" + "(round-remainder -123 -10) @result{} -3\n" + "(round-remainder 125 10) @result{} 5\n" + "(round-remainder 127 10) @result{} -3\n" + "(round-remainder 135 10) @result{} -5\n" + "(round-remainder -123.2 -63.5) @result{} 3.8\n" + "(round-remainder 16/3 -10/7) @result{} -8/21\n" "@end lisp") -#define FUNC_NAME s_scm_centered_remainder +#define FUNC_NAME s_scm_round_remainder { 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); if (SCM_UNLIKELY (yy == 0)) - scm_num_overflow (s_scm_centered_remainder); + scm_num_overflow (s_scm_round_remainder); else { - scm_t_inum xx = SCM_I_INUM (x); + scm_t_inum qq = xx / yy; scm_t_inum rr = xx % yy; - if (SCM_LIKELY (xx > 0)) + scm_t_inum ay = yy; + scm_t_inum r2 = 2 * rr; + + if (SCM_LIKELY (yy < 0)) { - if (SCM_LIKELY (yy > 0)) - { - if (rr >= (yy + 1) / 2) - rr -= yy; - } - else - { - if (rr >= (1 - yy) / 2) - rr += yy; - } + ay = -ay; + r2 = -r2; + } + + if (qq & 1L) + { + if (r2 >= ay) + rr -= yy; + else if (r2 <= -ay) + rr += yy; } else { - if (SCM_LIKELY (yy > 0)) - { - if (rr < -yy / 2) - rr += yy; - } - else - { - if (rr < yy / 2) - rr -= yy; - } + if (r2 > ay) + rr -= yy; + else if (r2 < -ay) + rr += yy; } return SCM_I_MAKINUM (rr); } @@ -1874,18 +3444,17 @@ SCM_PRIMITIVE_GENERIC (scm_centered_remainder, "centered-remainder", 2, 0, 0, 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_remainder */ - return scm_i_bigint_centered_remainder - (scm_i_long2big (SCM_I_INUM (x)), y); + can fit in a fixnum) to scm_i_bigint_round_remainder */ + return scm_i_bigint_round_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_round_remainder (xx, SCM_REAL_VALUE (y)); else if (SCM_FRACTIONP (y)) - return scm_i_slow_exact_centered_remainder (x, y); + return scm_i_exact_rational_round_remainder (x, y); else - SCM_WTA_DISPATCH_2 (g_scm_centered_remainder, x, y, SCM_ARG2, - s_scm_centered_remainder); + return scm_wta_dispatch_2 (g_scm_round_remainder, x, y, SCM_ARG2, + s_scm_round_remainder); } else if (SCM_BIGP (x)) { @@ -1893,230 +3462,232 @@ SCM_PRIMITIVE_GENERIC (scm_centered_remainder, "centered-remainder", 2, 0, 0, { scm_t_inum yy = SCM_I_INUM (y); if (SCM_UNLIKELY (yy == 0)) - scm_num_overflow (s_scm_centered_remainder); + scm_num_overflow (s_scm_round_remainder); else { + SCM q = scm_i_mkbig (); scm_t_inum rr; - /* Arrange for rr to initially be non-positive, - because that simplifies the test to see - if it is within the needed bounds. */ + int needs_adjustment; + if (yy > 0) { - rr = - mpz_cdiv_ui (SCM_I_BIG_MPZ (x), yy); - scm_remember_upto_here_1 (x); - if (rr < -yy / 2) - rr += yy; + rr = mpz_fdiv_q_ui (SCM_I_BIG_MPZ (q), + SCM_I_BIG_MPZ (x), yy); + if (mpz_odd_p (SCM_I_BIG_MPZ (q))) + needs_adjustment = (2*rr >= yy); + else + needs_adjustment = (2*rr > yy); } else { - rr = - mpz_cdiv_ui (SCM_I_BIG_MPZ (x), -yy); - scm_remember_upto_here_1 (x); - if (rr < yy / 2) - rr -= yy; + rr = - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q), + SCM_I_BIG_MPZ (x), -yy); + if (mpz_odd_p (SCM_I_BIG_MPZ (q))) + needs_adjustment = (2*rr <= yy); + else + needs_adjustment = (2*rr < yy); } + scm_remember_upto_here_2 (x, q); + if (needs_adjustment) + rr -= yy; return SCM_I_MAKINUM (rr); } } else if (SCM_BIGP (y)) - return scm_i_bigint_centered_remainder (x, y); + return scm_i_bigint_round_remainder (x, y); else if (SCM_REALP (y)) - return scm_i_inexact_centered_remainder + return scm_i_inexact_round_remainder (scm_i_big2dbl (x), SCM_REAL_VALUE (y)); else if (SCM_FRACTIONP (y)) - return scm_i_slow_exact_centered_remainder (x, y); + return scm_i_exact_rational_round_remainder (x, y); else - SCM_WTA_DISPATCH_2 (g_scm_centered_remainder, x, y, SCM_ARG2, - s_scm_centered_remainder); + return scm_wta_dispatch_2 (g_scm_round_remainder, x, y, SCM_ARG2, + s_scm_round_remainder); } else if (SCM_REALP (x)) { if (SCM_REALP (y) || SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) - return scm_i_inexact_centered_remainder + return scm_i_inexact_round_remainder (SCM_REAL_VALUE (x), scm_to_double (y)); else - SCM_WTA_DISPATCH_2 (g_scm_centered_remainder, x, y, SCM_ARG2, - s_scm_centered_remainder); + return scm_wta_dispatch_2 (g_scm_round_remainder, x, y, SCM_ARG2, + s_scm_round_remainder); } else if (SCM_FRACTIONP (x)) { if (SCM_REALP (y)) - return scm_i_inexact_centered_remainder + return scm_i_inexact_round_remainder (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_exact_rational_round_remainder (x, y); else - return scm_i_slow_exact_centered_remainder (x, y); + return scm_wta_dispatch_2 (g_scm_round_remainder, x, y, SCM_ARG2, + s_scm_round_remainder); } else - SCM_WTA_DISPATCH_2 (g_scm_centered_remainder, x, y, SCM_ARG1, - s_scm_centered_remainder); + return scm_wta_dispatch_2 (g_scm_round_remainder, x, y, SCM_ARG1, + s_scm_round_remainder); } #undef FUNC_NAME static SCM -scm_i_inexact_centered_remainder (double x, double y) +scm_i_inexact_round_remainder (double x, double y) { - double q; - /* Although it would be more efficient to use fmod here, we can't because it would in some cases produce results inconsistent with - scm_i_inexact_centered_quotient, such that x != r + q * y (not even + scm_i_inexact_round_quotient, such that x != r + q * y (not even close). In particular, when x-y/2 is very close to a multiple of - y, then r might be either -abs(y/2) or abs(y/2)-epsilon, but those - two cases must correspond to different choices of q. If quotient + y, then r might be either -abs(y/2) or abs(y/2), but those two + cases must correspond to different choices of q. If quotient chooses one and remainder chooses the other, it would be bad. */ - if (SCM_LIKELY (y > 0)) - q = floor (x/y + 0.5); - else if (SCM_LIKELY (y < 0)) - q = ceil (x/y - 0.5); - else if (y == 0) - scm_num_overflow (s_scm_centered_remainder); /* or return a NaN? */ + + if (SCM_UNLIKELY (y == 0)) + scm_num_overflow (s_scm_round_remainder); /* or return a NaN? */ else - return scm_nan (); - return scm_from_double (x - q * y); + { + double q = scm_c_round (x / y); + return scm_from_double (x - q * y); + } } /* Assumes that both x and y are bigints, though x might be able to fit into a fixnum. */ static SCM -scm_i_bigint_centered_remainder (SCM x, SCM y) +scm_i_bigint_round_remainder (SCM x, SCM y) { - SCM r, min_r; + SCM q, r, r2; + int cmp, needs_adjustment; /* Note that x might be small enough to fit into a fixnum, so we must not let it escape into the wild */ + q = scm_i_mkbig (); r = scm_i_mkbig (); + r2 = scm_i_mkbig (); - /* min_r will eventually become -abs(y)/2 */ - min_r = scm_i_mkbig (); - mpz_tdiv_q_2exp (SCM_I_BIG_MPZ (min_r), - SCM_I_BIG_MPZ (y), 1); + mpz_fdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r), + SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y)); + scm_remember_upto_here_1 (x); + mpz_mul_2exp (SCM_I_BIG_MPZ (r2), SCM_I_BIG_MPZ (r), 1); /* r2 = 2*r */ - /* Arrange for rr to initially be non-positive, - because that simplifies the test to see - if it is within the needed bounds. */ - if (mpz_sgn (SCM_I_BIG_MPZ (y)) > 0) - { - mpz_cdiv_r (SCM_I_BIG_MPZ (r), - SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y)); - mpz_neg (SCM_I_BIG_MPZ (min_r), SCM_I_BIG_MPZ (min_r)); - if (mpz_cmp (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (min_r)) < 0) - mpz_add (SCM_I_BIG_MPZ (r), - SCM_I_BIG_MPZ (r), - SCM_I_BIG_MPZ (y)); - } + cmp = mpz_cmpabs (SCM_I_BIG_MPZ (r2), SCM_I_BIG_MPZ (y)); + if (mpz_odd_p (SCM_I_BIG_MPZ (q))) + needs_adjustment = (cmp >= 0); else - { - mpz_fdiv_r (SCM_I_BIG_MPZ (r), - SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y)); - if (mpz_cmp (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (min_r)) < 0) - mpz_sub (SCM_I_BIG_MPZ (r), - SCM_I_BIG_MPZ (r), - SCM_I_BIG_MPZ (y)); - } - scm_remember_upto_here_2 (x, y); + needs_adjustment = (cmp > 0); + scm_remember_upto_here_2 (q, r2); + + if (needs_adjustment) + mpz_sub (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (y)); + + scm_remember_upto_here_1 (y); return scm_i_normbig (r); } -/* Compute exact centered_remainder the slow way. - We use this only if both arguments are exact, - and at least one of them is a fraction */ static SCM -scm_i_slow_exact_centered_remainder (SCM x, SCM y) +scm_i_exact_rational_round_remainder (SCM x, SCM y) { - SCM q; - - if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x))) - SCM_WTA_DISPATCH_2 (g_scm_centered_remainder, x, y, SCM_ARG1, - s_scm_centered_remainder); - else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y))) - SCM_WTA_DISPATCH_2 (g_scm_centered_remainder, x, y, SCM_ARG2, - s_scm_centered_remainder); - else if (scm_is_true (scm_positive_p (y))) - q = scm_floor (scm_sum (scm_divide (x, y), exactly_one_half)); - else if (scm_is_true (scm_negative_p (y))) - q = scm_ceiling (scm_difference (scm_divide (x, y), exactly_one_half)); - else - scm_num_overflow (s_scm_centered_remainder); - return scm_difference (x, scm_product (y, q)); + SCM xd = scm_denominator (x); + SCM yd = scm_denominator (y); + SCM r1 = scm_round_remainder (scm_product (scm_numerator (x), yd), + scm_product (scm_numerator (y), xd)); + return scm_divide (r1, scm_product (xd, yd)); } -static SCM scm_i_inexact_centered_divide (double x, double y); -static SCM scm_i_bigint_centered_divide (SCM x, SCM y); -static SCM scm_i_slow_exact_centered_divide (SCM x, SCM y); +static void scm_i_inexact_round_divide (double x, double y, SCM *qp, SCM *rp); +static void scm_i_bigint_round_divide (SCM x, SCM y, SCM *qp, SCM *rp); +static void scm_i_exact_rational_round_divide (SCM x, SCM y, SCM *qp, SCM *rp); -SCM_PRIMITIVE_GENERIC (scm_centered_divide, "centered/", 2, 0, 0, +SCM_PRIMITIVE_GENERIC (scm_i_round_divide, "round/", 2, 0, 0, (SCM x, SCM y), "Return the integer @var{q} and the real number @var{r}\n" "such that @math{@var{x} = @var{q}*@var{y} + @var{r}}\n" - "and @math{-abs(@var{y}/2) <= @var{r} < abs(@var{y}/2)}.\n" + "and @var{q} is @math{@var{x} / @var{y}} rounded to the\n" + "nearest integer, with ties going to the nearest even integer.\n" "@lisp\n" - "(centered/ 123 10) @result{} 12 and 3\n" - "(centered/ 123 -10) @result{} -12 and 3\n" - "(centered/ -123 10) @result{} -12 and -3\n" - "(centered/ -123 -10) @result{} 12 and -3\n" - "(centered/ -123.2 -63.5) @result{} 2.0 and 3.8\n" - "(centered/ 16/3 -10/7) @result{} -4 and -8/21\n" + "(round/ 123 10) @result{} 12 and 3\n" + "(round/ 123 -10) @result{} -12 and 3\n" + "(round/ -123 10) @result{} -12 and -3\n" + "(round/ -123 -10) @result{} 12 and -3\n" + "(round/ 125 10) @result{} 12 and 5\n" + "(round/ 127 10) @result{} 13 and -3\n" + "(round/ 135 10) @result{} 14 and -5\n" + "(round/ -123.2 -63.5) @result{} 2.0 and 3.8\n" + "(round/ 16/3 -10/7) @result{} -4 and -8/21\n" "@end lisp") -#define FUNC_NAME s_scm_centered_divide +#define FUNC_NAME s_scm_i_round_divide +{ + SCM q, r; + + scm_round_divide(x, y, &q, &r); + return scm_values (scm_list_2 (q, r)); +} +#undef FUNC_NAME + +#define s_scm_round_divide s_scm_i_round_divide +#define g_scm_round_divide g_scm_i_round_divide + +void +scm_round_divide (SCM x, SCM y, SCM *qp, SCM *rp) { 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); if (SCM_UNLIKELY (yy == 0)) - scm_num_overflow (s_scm_centered_divide); + scm_num_overflow (s_scm_round_divide); else { - scm_t_inum xx = SCM_I_INUM (x); scm_t_inum qq = xx / yy; - scm_t_inum rr = xx - qq * yy; - if (SCM_LIKELY (xx > 0)) + scm_t_inum rr = xx % yy; + scm_t_inum ay = yy; + scm_t_inum r2 = 2 * rr; + + if (SCM_LIKELY (yy < 0)) { - if (SCM_LIKELY (yy > 0)) - { - if (rr >= (yy + 1) / 2) - { qq++; rr -= yy; } - } - else - { - if (rr >= (1 - yy) / 2) - { qq--; rr += yy; } - } + ay = -ay; + r2 = -r2; + } + + if (qq & 1L) + { + if (r2 >= ay) + { qq++; rr -= yy; } + else if (r2 <= -ay) + { qq--; rr += yy; } } else { - if (SCM_LIKELY (yy > 0)) - { - if (rr < -yy / 2) - { qq--; rr += yy; } - } - else - { - if (rr < yy / 2) - { qq++; rr -= yy; } - } + if (r2 > ay) + { qq++; rr -= yy; } + else if (r2 < -ay) + { qq--; rr += yy; } } - return scm_values (scm_list_2 (SCM_I_MAKINUM (qq), - SCM_I_MAKINUM (rr))); + if (SCM_LIKELY (SCM_FIXABLE (qq))) + *qp = SCM_I_MAKINUM (qq); + else + *qp = scm_i_inum2big (qq); + *rp = SCM_I_MAKINUM (rr); } + return; } 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); + can fit in a fixnum) to scm_i_bigint_round_divide */ + return scm_i_bigint_round_divide + (scm_i_long2big (SCM_I_INUM (x)), y, qp, rp); } else if (SCM_REALP (y)) - return scm_i_inexact_centered_divide - (SCM_I_INUM (x), SCM_REAL_VALUE (y)); + return scm_i_inexact_round_divide (xx, SCM_REAL_VALUE (y), qp, rp); else if (SCM_FRACTIONP (y)) - return scm_i_slow_exact_centered_divide (x, y); + return scm_i_exact_rational_round_divide (x, y, qp, rp); else - SCM_WTA_DISPATCH_2 (g_scm_centered_divide, x, y, SCM_ARG2, - s_scm_centered_divide); + return two_valued_wta_dispatch_2 (g_scm_round_divide, x, y, SCM_ARG2, + s_scm_round_divide, qp, rp); } else if (SCM_BIGP (x)) { @@ -2124,172 +3695,141 @@ SCM_PRIMITIVE_GENERIC (scm_centered_divide, "centered/", 2, 0, 0, { scm_t_inum yy = SCM_I_INUM (y); if (SCM_UNLIKELY (yy == 0)) - scm_num_overflow (s_scm_centered_divide); + scm_num_overflow (s_scm_round_divide); else { SCM q = scm_i_mkbig (); scm_t_inum rr; - /* Arrange for rr to initially be non-positive, - because that simplifies the test to see - if it is within the needed bounds. */ + int needs_adjustment; + if (yy > 0) { - rr = - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q), - SCM_I_BIG_MPZ (x), yy); - scm_remember_upto_here_1 (x); - if (rr < -yy / 2) - { - mpz_sub_ui (SCM_I_BIG_MPZ (q), - SCM_I_BIG_MPZ (q), 1); - rr += yy; - } + rr = mpz_fdiv_q_ui (SCM_I_BIG_MPZ (q), + SCM_I_BIG_MPZ (x), yy); + if (mpz_odd_p (SCM_I_BIG_MPZ (q))) + needs_adjustment = (2*rr >= yy); + else + needs_adjustment = (2*rr > yy); } else { rr = - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (x), -yy); - scm_remember_upto_here_1 (x); mpz_neg (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (q)); - if (rr < yy / 2) - { - mpz_add_ui (SCM_I_BIG_MPZ (q), - SCM_I_BIG_MPZ (q), 1); - rr -= yy; - } + if (mpz_odd_p (SCM_I_BIG_MPZ (q))) + needs_adjustment = (2*rr <= yy); + else + needs_adjustment = (2*rr < yy); + } + scm_remember_upto_here_1 (x); + if (needs_adjustment) + { + mpz_add_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (q), 1); + rr -= yy; } - return scm_values (scm_list_2 (scm_i_normbig (q), - SCM_I_MAKINUM (rr))); + *qp = scm_i_normbig (q); + *rp = SCM_I_MAKINUM (rr); } + return; } else if (SCM_BIGP (y)) - return scm_i_bigint_centered_divide (x, y); + return scm_i_bigint_round_divide (x, y, qp, rp); else if (SCM_REALP (y)) - return scm_i_inexact_centered_divide - (scm_i_big2dbl (x), SCM_REAL_VALUE (y)); + return scm_i_inexact_round_divide + (scm_i_big2dbl (x), SCM_REAL_VALUE (y), qp, rp); else if (SCM_FRACTIONP (y)) - return scm_i_slow_exact_centered_divide (x, y); + return scm_i_exact_rational_round_divide (x, y, qp, rp); else - SCM_WTA_DISPATCH_2 (g_scm_centered_divide, x, y, SCM_ARG2, - s_scm_centered_divide); + return two_valued_wta_dispatch_2 (g_scm_round_divide, x, y, SCM_ARG2, + s_scm_round_divide, qp, rp); } else if (SCM_REALP (x)) { if (SCM_REALP (y) || SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) - return scm_i_inexact_centered_divide - (SCM_REAL_VALUE (x), scm_to_double (y)); - else - SCM_WTA_DISPATCH_2 (g_scm_centered_divide, x, y, SCM_ARG2, - s_scm_centered_divide); + return scm_i_inexact_round_divide + (SCM_REAL_VALUE (x), scm_to_double (y), qp, rp); + else + return two_valued_wta_dispatch_2 (g_scm_round_divide, x, y, SCM_ARG2, + s_scm_round_divide, qp, rp); } else if (SCM_FRACTIONP (x)) { if (SCM_REALP (y)) - return scm_i_inexact_centered_divide - (scm_i_fraction2double (x), SCM_REAL_VALUE (y)); + return scm_i_inexact_round_divide + (scm_i_fraction2double (x), SCM_REAL_VALUE (y), qp, rp); + else if (SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y)) + return scm_i_exact_rational_round_divide (x, y, qp, rp); else - return scm_i_slow_exact_centered_divide (x, y); + return two_valued_wta_dispatch_2 (g_scm_round_divide, x, y, SCM_ARG2, + s_scm_round_divide, qp, rp); } else - SCM_WTA_DISPATCH_2 (g_scm_centered_divide, x, y, SCM_ARG1, - s_scm_centered_divide); + return two_valued_wta_dispatch_2 (g_scm_round_divide, x, y, SCM_ARG1, + s_scm_round_divide, qp, rp); } -#undef FUNC_NAME -static SCM -scm_i_inexact_centered_divide (double x, double y) +static void +scm_i_inexact_round_divide (double x, double y, SCM *qp, SCM *rp) { - double q, r; - - if (SCM_LIKELY (y > 0)) - q = floor (x/y + 0.5); - else if (SCM_LIKELY (y < 0)) - q = ceil (x/y - 0.5); - else if (y == 0) - scm_num_overflow (s_scm_centered_divide); /* or return a NaN? */ + if (SCM_UNLIKELY (y == 0)) + scm_num_overflow (s_scm_round_divide); /* or return a NaN? */ else - q = guile_NaN; - r = x - q * y; - return scm_values (scm_list_2 (scm_from_double (q), - scm_from_double (r))); + { + double q = scm_c_round (x / y); + double r = x - q * y; + *qp = scm_from_double (q); + *rp = scm_from_double (r); + } } /* Assumes that both x and y are bigints, though x might be able to fit into a fixnum. */ -static SCM -scm_i_bigint_centered_divide (SCM x, SCM y) +static void +scm_i_bigint_round_divide (SCM x, SCM y, SCM *qp, SCM *rp) { - SCM q, r, min_r; + SCM q, r, r2; + int cmp, needs_adjustment; /* Note that x might be small enough to fit into a fixnum, so we must not let it escape into the wild */ q = scm_i_mkbig (); r = scm_i_mkbig (); + r2 = scm_i_mkbig (); - /* min_r will eventually become -abs(y/2) */ - min_r = scm_i_mkbig (); - mpz_tdiv_q_2exp (SCM_I_BIG_MPZ (min_r), - SCM_I_BIG_MPZ (y), 1); + mpz_fdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r), + SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y)); + scm_remember_upto_here_1 (x); + mpz_mul_2exp (SCM_I_BIG_MPZ (r2), SCM_I_BIG_MPZ (r), 1); /* r2 = 2*r */ - /* Arrange for rr to initially be non-positive, - because that simplifies the test to see - if it is within the needed bounds. */ - if (mpz_sgn (SCM_I_BIG_MPZ (y)) > 0) - { - mpz_cdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r), - SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y)); - mpz_neg (SCM_I_BIG_MPZ (min_r), SCM_I_BIG_MPZ (min_r)); - if (mpz_cmp (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (min_r)) < 0) - { - mpz_sub_ui (SCM_I_BIG_MPZ (q), - SCM_I_BIG_MPZ (q), 1); - mpz_add (SCM_I_BIG_MPZ (r), - SCM_I_BIG_MPZ (r), - SCM_I_BIG_MPZ (y)); - } - } + cmp = mpz_cmpabs (SCM_I_BIG_MPZ (r2), SCM_I_BIG_MPZ (y)); + if (mpz_odd_p (SCM_I_BIG_MPZ (q))) + needs_adjustment = (cmp >= 0); else + needs_adjustment = (cmp > 0); + + if (needs_adjustment) { - mpz_fdiv_qr (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (r), - SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y)); - if (mpz_cmp (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (min_r)) < 0) - { - mpz_add_ui (SCM_I_BIG_MPZ (q), - SCM_I_BIG_MPZ (q), 1); - mpz_sub (SCM_I_BIG_MPZ (r), - SCM_I_BIG_MPZ (r), - SCM_I_BIG_MPZ (y)); - } + mpz_add_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (q), 1); + mpz_sub (SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (y)); } - scm_remember_upto_here_2 (x, y); - return scm_values (scm_list_2 (scm_i_normbig (q), - scm_i_normbig (r))); + + scm_remember_upto_here_2 (r2, y); + *qp = scm_i_normbig (q); + *rp = scm_i_normbig (r); } -/* Compute exact centered quotient and remainder the slow way. - We use this only if both arguments are exact, - and at least one of them is a fraction */ -static SCM -scm_i_slow_exact_centered_divide (SCM x, SCM y) +static void +scm_i_exact_rational_round_divide (SCM x, SCM y, SCM *qp, SCM *rp) { - SCM q, r; + SCM r1; + SCM xd = scm_denominator (x); + SCM yd = scm_denominator (y); - if (!(SCM_I_INUMP (x) || SCM_BIGP (x) || SCM_FRACTIONP (x))) - SCM_WTA_DISPATCH_2 (g_scm_centered_divide, x, y, SCM_ARG1, - s_scm_centered_divide); - else if (!(SCM_I_INUMP (y) || SCM_BIGP (y) || SCM_FRACTIONP (y))) - SCM_WTA_DISPATCH_2 (g_scm_centered_divide, x, y, SCM_ARG2, - s_scm_centered_divide); - else if (scm_is_true (scm_positive_p (y))) - q = scm_floor (scm_sum (scm_divide (x, y), - exactly_one_half)); - else if (scm_is_true (scm_negative_p (y))) - q = scm_ceiling (scm_difference (scm_divide (x, y), - exactly_one_half)); - else - scm_num_overflow (s_scm_centered_divide); - r = scm_difference (x, scm_product (q, y)); - return scm_values (scm_list_2 (q, r)); + scm_round_divide (scm_product (scm_numerator (x), yd), + scm_product (scm_numerator (y), xd), + qp, &r1); + *rp = scm_divide (r1, scm_product (xd, yd)); } @@ -2371,7 +3911,7 @@ scm_gcd (SCM x, SCM y) goto big_inum; } else - SCM_WTA_DISPATCH_2 (g_gcd, x, y, SCM_ARG2, s_gcd); + return scm_wta_dispatch_2 (g_gcd, x, y, SCM_ARG2, s_gcd); } else if (SCM_BIGP (x)) { @@ -2401,10 +3941,10 @@ scm_gcd (SCM x, SCM y) return scm_i_normbig (result); } else - SCM_WTA_DISPATCH_2 (g_gcd, x, y, SCM_ARG2, s_gcd); + return scm_wta_dispatch_2 (g_gcd, x, y, SCM_ARG2, s_gcd); } else - SCM_WTA_DISPATCH_2 (g_gcd, x, y, SCM_ARG1, s_gcd); + return scm_wta_dispatch_2 (g_gcd, x, y, SCM_ARG1, s_gcd); } SCM_PRIMITIVE_GENERIC (scm_i_lcm, "lcm", 0, 2, 1, @@ -2435,10 +3975,11 @@ scm_lcm (SCM n1, SCM n2) n2 = SCM_I_MAKINUM (1L); } - SCM_GASSERT2 (SCM_I_INUMP (n1) || SCM_BIGP (n1), - g_lcm, n1, n2, SCM_ARG1, s_lcm); - SCM_GASSERT2 (SCM_I_INUMP (n2) || SCM_BIGP (n2), - g_lcm, n1, n2, SCM_ARGn, s_lcm); + if (SCM_UNLIKELY (!(SCM_I_INUMP (n1) || SCM_BIGP (n1)))) + return scm_wta_dispatch_2 (g_lcm, n1, n2, SCM_ARG1, s_lcm); + + if (SCM_UNLIKELY (!(SCM_I_INUMP (n2) || SCM_BIGP (n2)))) + return scm_wta_dispatch_2 (g_lcm, n1, n2, SCM_ARG2, s_lcm); if (SCM_I_INUMP (n1)) { @@ -2571,7 +4112,7 @@ SCM scm_logand (SCM n1, SCM n2) else if SCM_BIGP (n2) { intbig: - if (n1 == 0) + if (nn1 == 0) return SCM_INUM0; { SCM result_z = scm_i_mkbig (); @@ -3568,12 +5109,6 @@ idbl2str (double f, char *a, int radix) exp++; } zero: -#ifdef ENGNOT - /* adding 9999 makes this equivalent to abs(x) % 3 */ - dpt = (exp + 9999) % 3; - exp -= dpt++; - efmt = 1; -#else efmt = (exp < -3) || (exp > wp + 2); if (!efmt) { @@ -3590,7 +5125,6 @@ idbl2str (double f, char *a, int radix) } else dpt = 1; -#endif do { @@ -3612,7 +5146,6 @@ idbl2str (double f, char *a, int radix) if (dpt > 0) { -#ifndef ENGNOT if ((dpt > 4) && (exp > 6)) { d = (a[0] == '-' ? 2 : 1); @@ -3622,7 +5155,6 @@ idbl2str (double f, char *a, int radix) efmt = 1; } else -#endif { while (--dpt) a[ch++] = '0'; @@ -3834,14 +5366,15 @@ scm_bigprint (SCM exp, SCM port, scm_print_state *pstate SCM_UNUSED) * in R5RS. Thus, the functions resemble syntactic units (, * , ...) that are used to build up numbers in the grammar. Some * points should be noted about the implementation: + * * * Each function keeps a local index variable 'idx' that points at the * current position within the parsed string. The global index is only * updated if the function could parse the corresponding syntactic unit * successfully. + * * * Similarly, the functions keep track of indicators of inexactness ('#', - * '.' or exponents) using local variables ('hash_seen', 'x'). Again, the - * global exactness information is only updated after each part has been - * successfully parsed. + * '.' or exponents) using local variables ('hash_seen', 'x'). + * * * Sequences of digits are parsed into temporary variables holding fixnums. * Only if these fixnums would overflow, the result variables are updated * using the standard functions scm_add, scm_product, scm_divide etc. Then, @@ -3850,6 +5383,34 @@ scm_bigprint (SCM exp, SCM port, scm_print_state *pstate SCM_UNUSED) * digits, a number 1234567890 would be parsed in two parts 12345 and 67890, * and the result was computed as 12345 * 100000 + 67890. In other words, * only every five digits two bignum operations were performed. + * + * Notes on the handling of exactness specifiers: + * + * When parsing non-real complex numbers, we apply exactness specifiers on + * per-component basis, as is done in PLT Scheme. For complex numbers + * written in rectangular form, exactness specifiers are applied to the + * real and imaginary parts before calling scm_make_rectangular. For + * complex numbers written in polar form, exactness specifiers are applied + * to the magnitude and angle before calling scm_make_polar. + * + * There are two kinds of exactness specifiers: forced and implicit. A + * forced exactness specifier is a "#e" or "#i" prefix at the beginning of + * the entire number, and applies to both components of a complex number. + * "#e" causes each component to be made exact, and "#i" causes each + * component to be made inexact. If no forced exactness specifier is + * present, then the exactness of each component is determined + * independently by the presence or absence of a decimal point or hash mark + * within that component. If a decimal point or hash mark is present, the + * component is made inexact, otherwise it is made exact. + * + * After the exactness specifiers have been applied to each component, they + * are passed to either scm_make_rectangular or scm_make_polar to produce + * the final result. Note that this will result in a real number if the + * imaginary part, magnitude, or angle is an exact 0. + * + * For example, (string->number "#i5.0+0i") does the equivalent of: + * + * (make-rectangular (exact->inexact 5) (exact->inexact 0)) */ enum t_exactness {NO_EXACTNESS, INEXACT, EXACT}; @@ -4100,7 +5661,7 @@ mem2decimal_from_point (SCM result, SCM mem, if (sign == 1) result = scm_product (result, e); else - result = scm_divide2real (result, e); + result = scm_divide (result, e); /* We've seen an exponent, thus the value is implicitly inexact. */ x = INEXACT; @@ -4616,7 +6177,8 @@ scm_num_eq_p (SCM x, SCM y) else if (SCM_FRACTIONP (y)) return SCM_BOOL_F; else - SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p, x, y, SCM_ARGn, s_scm_i_num_eq_p); + return scm_wta_dispatch_2 (g_scm_i_num_eq_p, x, y, SCM_ARGn, + s_scm_i_num_eq_p); } else if (SCM_BIGP (x)) { @@ -4651,7 +6213,8 @@ scm_num_eq_p (SCM x, SCM y) else if (SCM_FRACTIONP (y)) return SCM_BOOL_F; else - SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p, x, y, SCM_ARGn, s_scm_i_num_eq_p); + return scm_wta_dispatch_2 (g_scm_i_num_eq_p, x, y, SCM_ARGn, + s_scm_i_num_eq_p); } else if (SCM_REALP (x)) { @@ -4689,7 +6252,8 @@ scm_num_eq_p (SCM x, SCM y) goto again; } else - SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p, x, y, SCM_ARGn, s_scm_i_num_eq_p); + return scm_wta_dispatch_2 (g_scm_i_num_eq_p, x, y, SCM_ARGn, + s_scm_i_num_eq_p); } else if (SCM_COMPLEXP (x)) { @@ -4727,7 +6291,8 @@ scm_num_eq_p (SCM x, SCM y) goto again; } else - SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p, x, y, SCM_ARGn, s_scm_i_num_eq_p); + return scm_wta_dispatch_2 (g_scm_i_num_eq_p, x, y, SCM_ARGn, + s_scm_i_num_eq_p); } else if (SCM_FRACTIONP (x)) { @@ -4761,10 +6326,12 @@ scm_num_eq_p (SCM x, SCM y) else if (SCM_FRACTIONP (y)) return scm_i_fraction_equalp (x, y); else - SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p, x, y, SCM_ARGn, s_scm_i_num_eq_p); + return scm_wta_dispatch_2 (g_scm_i_num_eq_p, x, y, SCM_ARGn, + s_scm_i_num_eq_p); } else - SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p, x, y, SCM_ARG1, s_scm_i_num_eq_p); + return scm_wta_dispatch_2 (g_scm_i_num_eq_p, x, y, SCM_ARG1, + s_scm_i_num_eq_p); } @@ -4823,7 +6390,8 @@ scm_less_p (SCM x, SCM y) goto again; } else - SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p, x, y, SCM_ARGn, s_scm_i_num_less_p); + return scm_wta_dispatch_2 (g_scm_i_num_less_p, x, y, SCM_ARGn, + s_scm_i_num_less_p); } else if (SCM_BIGP (x)) { @@ -4851,7 +6419,8 @@ scm_less_p (SCM x, SCM y) else if (SCM_FRACTIONP (y)) goto int_frac; else - SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p, x, y, SCM_ARGn, s_scm_i_num_less_p); + return scm_wta_dispatch_2 (g_scm_i_num_less_p, x, y, SCM_ARGn, + s_scm_i_num_less_p); } else if (SCM_REALP (x)) { @@ -4879,7 +6448,8 @@ scm_less_p (SCM x, SCM y) goto again; } else - SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p, x, y, SCM_ARGn, s_scm_i_num_less_p); + return scm_wta_dispatch_2 (g_scm_i_num_less_p, x, y, SCM_ARGn, + s_scm_i_num_less_p); } else if (SCM_FRACTIONP (x)) { @@ -4912,10 +6482,12 @@ scm_less_p (SCM x, SCM y) goto again; } else - SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p, x, y, SCM_ARGn, s_scm_i_num_less_p); + return scm_wta_dispatch_2 (g_scm_i_num_less_p, x, y, SCM_ARGn, + s_scm_i_num_less_p); } else - SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p, x, y, SCM_ARG1, s_scm_i_num_less_p); + return scm_wta_dispatch_2 (g_scm_i_num_less_p, x, y, SCM_ARG1, + s_scm_i_num_less_p); } @@ -4944,9 +6516,9 @@ SCM scm_gr_p (SCM x, SCM y) { if (!SCM_NUMBERP (x)) - SCM_WTA_DISPATCH_2 (g_scm_i_num_gr_p, x, y, SCM_ARG1, FUNC_NAME); + return scm_wta_dispatch_2 (g_scm_i_num_gr_p, x, y, SCM_ARG1, FUNC_NAME); else if (!SCM_NUMBERP (y)) - SCM_WTA_DISPATCH_2 (g_scm_i_num_gr_p, x, y, SCM_ARG2, FUNC_NAME); + return scm_wta_dispatch_2 (g_scm_i_num_gr_p, x, y, SCM_ARG2, FUNC_NAME); else return scm_less_p (y, x); } @@ -4978,9 +6550,9 @@ SCM scm_leq_p (SCM x, SCM y) { if (!SCM_NUMBERP (x)) - SCM_WTA_DISPATCH_2 (g_scm_i_num_leq_p, x, y, SCM_ARG1, FUNC_NAME); + return scm_wta_dispatch_2 (g_scm_i_num_leq_p, x, y, SCM_ARG1, FUNC_NAME); else if (!SCM_NUMBERP (y)) - SCM_WTA_DISPATCH_2 (g_scm_i_num_leq_p, x, y, SCM_ARG2, FUNC_NAME); + return scm_wta_dispatch_2 (g_scm_i_num_leq_p, x, y, SCM_ARG2, FUNC_NAME); else if (scm_is_true (scm_nan_p (x)) || scm_is_true (scm_nan_p (y))) return SCM_BOOL_F; else @@ -5014,9 +6586,9 @@ SCM scm_geq_p (SCM x, SCM y) { if (!SCM_NUMBERP (x)) - SCM_WTA_DISPATCH_2 (g_scm_i_num_geq_p, x, y, SCM_ARG1, FUNC_NAME); + return scm_wta_dispatch_2 (g_scm_i_num_geq_p, x, y, SCM_ARG1, FUNC_NAME); else if (!SCM_NUMBERP (y)) - SCM_WTA_DISPATCH_2 (g_scm_i_num_geq_p, x, y, SCM_ARG2, FUNC_NAME); + return scm_wta_dispatch_2 (g_scm_i_num_geq_p, x, y, SCM_ARG2, FUNC_NAME); else if (scm_is_true (scm_nan_p (x)) || scm_is_true (scm_nan_p (y))) return SCM_BOOL_F; else @@ -5043,7 +6615,7 @@ SCM_PRIMITIVE_GENERIC (scm_zero_p, "zero?", 1, 0, 0, else if (SCM_FRACTIONP (z)) return SCM_BOOL_F; else - SCM_WTA_DISPATCH_1 (g_scm_zero_p, z, SCM_ARG1, s_scm_zero_p); + return scm_wta_dispatch_1 (g_scm_zero_p, z, SCM_ARG1, s_scm_zero_p); } #undef FUNC_NAME @@ -5067,7 +6639,7 @@ SCM_PRIMITIVE_GENERIC (scm_positive_p, "positive?", 1, 0, 0, else if (SCM_FRACTIONP (x)) return scm_positive_p (SCM_FRACTION_NUMERATOR (x)); else - SCM_WTA_DISPATCH_1 (g_scm_positive_p, x, SCM_ARG1, s_scm_positive_p); + return scm_wta_dispatch_1 (g_scm_positive_p, x, SCM_ARG1, s_scm_positive_p); } #undef FUNC_NAME @@ -5091,7 +6663,7 @@ SCM_PRIMITIVE_GENERIC (scm_negative_p, "negative?", 1, 0, 0, else if (SCM_FRACTIONP (x)) return scm_negative_p (SCM_FRACTION_NUMERATOR (x)); else - SCM_WTA_DISPATCH_1 (g_scm_negative_p, x, SCM_ARG1, s_scm_negative_p); + return scm_wta_dispatch_1 (g_scm_negative_p, x, SCM_ARG1, s_scm_negative_p); } #undef FUNC_NAME @@ -5125,11 +6697,11 @@ scm_max (SCM x, SCM y) if (SCM_UNBNDP (y)) { if (SCM_UNBNDP (x)) - SCM_WTA_DISPATCH_0 (g_max, s_max); + return scm_wta_dispatch_0 (g_max, s_max); else if (SCM_I_INUMP(x) || SCM_BIGP(x) || SCM_REALP(x) || SCM_FRACTIONP(x)) return x; else - SCM_WTA_DISPATCH_1 (g_max, x, SCM_ARG1, s_max); + return scm_wta_dispatch_1 (g_max, x, SCM_ARG1, s_max); } if (SCM_I_INUMP (x)) @@ -5168,7 +6740,7 @@ scm_max (SCM x, SCM y) return (scm_is_false (scm_less_p (x, y)) ? x : y); } else - SCM_WTA_DISPATCH_2 (g_max, x, y, SCM_ARGn, s_max); + return scm_wta_dispatch_2 (g_max, x, y, SCM_ARGn, s_max); } else if (SCM_BIGP (x)) { @@ -5198,7 +6770,7 @@ scm_max (SCM x, SCM y) goto use_less; } else - SCM_WTA_DISPATCH_2 (g_max, x, y, SCM_ARGn, s_max); + return scm_wta_dispatch_2 (g_max, x, y, SCM_ARGn, s_max); } else if (SCM_REALP (x)) { @@ -5253,7 +6825,7 @@ scm_max (SCM x, SCM y) return (xx < yy) ? scm_from_double (yy) : x; } else - SCM_WTA_DISPATCH_2 (g_max, x, y, SCM_ARGn, s_max); + return scm_wta_dispatch_2 (g_max, x, y, SCM_ARGn, s_max); } else if (SCM_FRACTIONP (x)) { @@ -5276,10 +6848,10 @@ scm_max (SCM x, SCM y) goto use_less; } else - SCM_WTA_DISPATCH_2 (g_max, x, y, SCM_ARGn, s_max); + return scm_wta_dispatch_2 (g_max, x, y, SCM_ARGn, s_max); } else - SCM_WTA_DISPATCH_2 (g_max, x, y, SCM_ARG1, s_max); + return scm_wta_dispatch_2 (g_max, x, y, SCM_ARG1, s_max); } @@ -5306,11 +6878,11 @@ scm_min (SCM x, SCM y) if (SCM_UNBNDP (y)) { if (SCM_UNBNDP (x)) - SCM_WTA_DISPATCH_0 (g_min, s_min); + return scm_wta_dispatch_0 (g_min, s_min); else if (SCM_I_INUMP(x) || SCM_BIGP(x) || SCM_REALP(x) || SCM_FRACTIONP(x)) return x; else - SCM_WTA_DISPATCH_1 (g_min, x, SCM_ARG1, s_min); + return scm_wta_dispatch_1 (g_min, x, SCM_ARG1, s_min); } if (SCM_I_INUMP (x)) @@ -5339,7 +6911,7 @@ scm_min (SCM x, SCM y) return (scm_is_false (scm_less_p (x, y)) ? y : x); } else - SCM_WTA_DISPATCH_2 (g_min, x, y, SCM_ARGn, s_min); + return scm_wta_dispatch_2 (g_min, x, y, SCM_ARGn, s_min); } else if (SCM_BIGP (x)) { @@ -5369,7 +6941,7 @@ scm_min (SCM x, SCM y) goto use_less; } else - SCM_WTA_DISPATCH_2 (g_min, x, y, SCM_ARGn, s_min); + return scm_wta_dispatch_2 (g_min, x, y, SCM_ARGn, s_min); } else if (SCM_REALP (x)) { @@ -5413,7 +6985,7 @@ scm_min (SCM x, SCM y) return (yy < xx) ? scm_from_double (yy) : x; } else - SCM_WTA_DISPATCH_2 (g_min, x, y, SCM_ARGn, s_min); + return scm_wta_dispatch_2 (g_min, x, y, SCM_ARGn, s_min); } else if (SCM_FRACTIONP (x)) { @@ -5436,10 +7008,10 @@ scm_min (SCM x, SCM y) goto use_less; } else - SCM_WTA_DISPATCH_2 (g_min, x, y, SCM_ARGn, s_min); + return scm_wta_dispatch_2 (g_min, x, y, SCM_ARGn, s_min); } else - SCM_WTA_DISPATCH_2 (g_min, x, y, SCM_ARG1, s_min); + return scm_wta_dispatch_2 (g_min, x, y, SCM_ARG1, s_min); } @@ -5468,7 +7040,7 @@ scm_sum (SCM x, SCM y) { if (SCM_NUMBERP (x)) return x; if (SCM_UNBNDP (x)) return SCM_INUM0; - SCM_WTA_DISPATCH_1 (g_sum, x, SCM_ARG1, s_sum); + return scm_wta_dispatch_1 (g_sum, x, SCM_ARG1, s_sum); } if (SCM_LIKELY (SCM_I_INUMP (x))) @@ -5501,7 +7073,7 @@ scm_sum (SCM x, SCM y) scm_product (x, SCM_FRACTION_DENOMINATOR (y))), SCM_FRACTION_DENOMINATOR (y)); else - SCM_WTA_DISPATCH_2 (g_sum, x, y, SCM_ARGn, s_sum); + return scm_wta_dispatch_2 (g_sum, x, y, SCM_ARGn, s_sum); } else if (SCM_BIGP (x)) { if (SCM_I_INUMP (y)) @@ -5566,7 +7138,7 @@ scm_sum (SCM x, SCM y) scm_product (x, SCM_FRACTION_DENOMINATOR (y))), SCM_FRACTION_DENOMINATOR (y)); else - SCM_WTA_DISPATCH_2 (g_sum, x, y, SCM_ARGn, s_sum); + return scm_wta_dispatch_2 (g_sum, x, y, SCM_ARGn, s_sum); } else if (SCM_REALP (x)) { @@ -5586,7 +7158,7 @@ scm_sum (SCM x, SCM y) else if (SCM_FRACTIONP (y)) return scm_from_double (SCM_REAL_VALUE (x) + scm_i_fraction2double (y)); else - SCM_WTA_DISPATCH_2 (g_sum, x, y, SCM_ARGn, s_sum); + return scm_wta_dispatch_2 (g_sum, x, y, SCM_ARGn, s_sum); } else if (SCM_COMPLEXP (x)) { @@ -5610,7 +7182,7 @@ scm_sum (SCM x, SCM y) return scm_c_make_rectangular (SCM_COMPLEX_REAL (x) + scm_i_fraction2double (y), SCM_COMPLEX_IMAG (x)); else - SCM_WTA_DISPATCH_2 (g_sum, x, y, SCM_ARGn, s_sum); + return scm_wta_dispatch_2 (g_sum, x, y, SCM_ARGn, s_sum); } else if (SCM_FRACTIONP (x)) { @@ -5633,10 +7205,10 @@ scm_sum (SCM x, SCM y) scm_product (SCM_FRACTION_NUMERATOR (y), SCM_FRACTION_DENOMINATOR (x))), scm_product (SCM_FRACTION_DENOMINATOR (x), SCM_FRACTION_DENOMINATOR (y))); else - SCM_WTA_DISPATCH_2 (g_sum, x, y, SCM_ARGn, s_sum); + return scm_wta_dispatch_2 (g_sum, x, y, SCM_ARGn, s_sum); } else - SCM_WTA_DISPATCH_2 (g_sum, x, y, SCM_ARG1, s_sum); + return scm_wta_dispatch_2 (g_sum, x, y, SCM_ARG1, s_sum); } @@ -5676,7 +7248,7 @@ scm_difference (SCM x, SCM y) if (SCM_UNLIKELY (SCM_UNBNDP (y))) { if (SCM_UNBNDP (x)) - SCM_WTA_DISPATCH_0 (g_difference, s_difference); + return scm_wta_dispatch_0 (g_difference, s_difference); else if (SCM_I_INUMP (x)) { @@ -5699,7 +7271,7 @@ scm_difference (SCM x, SCM y) return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x), SCM_UNDEFINED), SCM_FRACTION_DENOMINATOR (x)); else - SCM_WTA_DISPATCH_1 (g_difference, x, SCM_ARG1, s_difference); + return scm_wta_dispatch_1 (g_difference, x, SCM_ARG1, s_difference); } if (SCM_LIKELY (SCM_I_INUMP (x))) @@ -5786,7 +7358,7 @@ scm_difference (SCM x, SCM y) SCM_FRACTION_NUMERATOR (y)), SCM_FRACTION_DENOMINATOR (y)); else - SCM_WTA_DISPATCH_2 (g_difference, x, y, SCM_ARGn, s_difference); + return scm_wta_dispatch_2 (g_difference, x, y, SCM_ARGn, s_difference); } else if (SCM_BIGP (x)) { @@ -5850,7 +7422,8 @@ scm_difference (SCM x, SCM y) return scm_i_make_ratio (scm_difference (scm_product (x, SCM_FRACTION_DENOMINATOR (y)), SCM_FRACTION_NUMERATOR (y)), SCM_FRACTION_DENOMINATOR (y)); - else SCM_WTA_DISPATCH_2 (g_difference, x, y, SCM_ARGn, s_difference); + else + return scm_wta_dispatch_2 (g_difference, x, y, SCM_ARGn, s_difference); } else if (SCM_REALP (x)) { @@ -5870,7 +7443,7 @@ scm_difference (SCM x, SCM y) else if (SCM_FRACTIONP (y)) return scm_from_double (SCM_REAL_VALUE (x) - scm_i_fraction2double (y)); else - SCM_WTA_DISPATCH_2 (g_difference, x, y, SCM_ARGn, s_difference); + return scm_wta_dispatch_2 (g_difference, x, y, SCM_ARGn, s_difference); } else if (SCM_COMPLEXP (x)) { @@ -5894,7 +7467,7 @@ scm_difference (SCM x, SCM y) return scm_c_make_rectangular (SCM_COMPLEX_REAL (x) - scm_i_fraction2double (y), SCM_COMPLEX_IMAG (x)); else - SCM_WTA_DISPATCH_2 (g_difference, x, y, SCM_ARGn, s_difference); + return scm_wta_dispatch_2 (g_difference, x, y, SCM_ARGn, s_difference); } else if (SCM_FRACTIONP (x)) { @@ -5918,10 +7491,10 @@ scm_difference (SCM x, SCM y) scm_product (SCM_FRACTION_NUMERATOR (y), SCM_FRACTION_DENOMINATOR (x))), scm_product (SCM_FRACTION_DENOMINATOR (x), SCM_FRACTION_DENOMINATOR (y))); else - SCM_WTA_DISPATCH_2 (g_difference, x, y, SCM_ARGn, s_difference); + return scm_wta_dispatch_2 (g_difference, x, y, SCM_ARGn, s_difference); } else - SCM_WTA_DISPATCH_2 (g_difference, x, y, SCM_ARG1, s_difference); + return scm_wta_dispatch_2 (g_difference, x, y, SCM_ARG1, s_difference); } #undef FUNC_NAME @@ -5964,7 +7537,7 @@ scm_product (SCM x, SCM y) else if (SCM_NUMBERP (x)) return x; else - SCM_WTA_DISPATCH_1 (g_product, x, SCM_ARG1, s_product); + return scm_wta_dispatch_1 (g_product, x, SCM_ARG1, s_product); } if (SCM_LIKELY (SCM_I_INUMP (x))) @@ -5997,7 +7570,7 @@ scm_product (SCM x, SCM y) else if (SCM_NUMP (y)) return SCM_INUM0; else - SCM_WTA_DISPATCH_2 (g_product, x, y, SCM_ARGn, s_product); + return scm_wta_dispatch_2 (g_product, x, y, SCM_ARGn, s_product); break; case -1: /* @@ -6042,7 +7615,7 @@ scm_product (SCM x, SCM y) return scm_i_make_ratio (scm_product (x, SCM_FRACTION_NUMERATOR (y)), SCM_FRACTION_DENOMINATOR (y)); else - SCM_WTA_DISPATCH_2 (g_product, x, y, SCM_ARGn, s_product); + return scm_wta_dispatch_2 (g_product, x, y, SCM_ARGn, s_product); } else if (SCM_BIGP (x)) { @@ -6077,7 +7650,7 @@ scm_product (SCM x, SCM y) return scm_i_make_ratio (scm_product (x, SCM_FRACTION_NUMERATOR (y)), SCM_FRACTION_DENOMINATOR (y)); else - SCM_WTA_DISPATCH_2 (g_product, x, y, SCM_ARGn, s_product); + return scm_wta_dispatch_2 (g_product, x, y, SCM_ARGn, s_product); } else if (SCM_REALP (x)) { @@ -6100,7 +7673,7 @@ scm_product (SCM x, SCM y) else if (SCM_FRACTIONP (y)) return scm_from_double (SCM_REAL_VALUE (x) * scm_i_fraction2double (y)); else - SCM_WTA_DISPATCH_2 (g_product, x, y, SCM_ARGn, s_product); + return scm_wta_dispatch_2 (g_product, x, y, SCM_ARGn, s_product); } else if (SCM_COMPLEXP (x)) { @@ -6133,7 +7706,7 @@ scm_product (SCM x, SCM y) yy * SCM_COMPLEX_IMAG (x)); } else - SCM_WTA_DISPATCH_2 (g_product, x, y, SCM_ARGn, s_product); + return scm_wta_dispatch_2 (g_product, x, y, SCM_ARGn, s_product); } else if (SCM_FRACTIONP (x)) { @@ -6158,10 +7731,10 @@ scm_product (SCM x, SCM y) scm_product (SCM_FRACTION_DENOMINATOR (x), SCM_FRACTION_DENOMINATOR (y))); else - SCM_WTA_DISPATCH_2 (g_product, x, y, SCM_ARGn, s_product); + return scm_wta_dispatch_2 (g_product, x, y, SCM_ARGn, s_product); } else - SCM_WTA_DISPATCH_2 (g_product, x, y, SCM_ARG1, s_product); + return scm_wta_dispatch_2 (g_product, x, y, SCM_ARG1, s_product); } #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \ @@ -6225,7 +7798,7 @@ do_divide (SCM x, SCM y, int inexact) if (SCM_UNLIKELY (SCM_UNBNDP (y))) { if (SCM_UNBNDP (x)) - SCM_WTA_DISPATCH_0 (g_divide, s_divide); + return scm_wta_dispatch_0 (g_divide, s_divide); else if (SCM_I_INUMP (x)) { scm_t_inum xx = SCM_I_INUM (x); @@ -6279,7 +7852,7 @@ do_divide (SCM x, SCM y, int inexact) return scm_i_make_ratio (SCM_FRACTION_DENOMINATOR (x), SCM_FRACTION_NUMERATOR (x)); else - SCM_WTA_DISPATCH_1 (g_divide, x, SCM_ARG1, s_divide); + return scm_wta_dispatch_1 (g_divide, x, SCM_ARG1, s_divide); } if (SCM_LIKELY (SCM_I_INUMP (x))) @@ -6353,7 +7926,7 @@ do_divide (SCM x, SCM y, int inexact) return scm_i_make_ratio (scm_product (x, SCM_FRACTION_DENOMINATOR (y)), SCM_FRACTION_NUMERATOR (y)); else - SCM_WTA_DISPATCH_2 (g_divide, x, y, SCM_ARGn, s_divide); + return scm_wta_dispatch_2 (g_divide, x, y, SCM_ARGn, s_divide); } else if (SCM_BIGP (x)) { @@ -6452,7 +8025,7 @@ do_divide (SCM x, SCM y, int inexact) return scm_i_make_ratio (scm_product (x, SCM_FRACTION_DENOMINATOR (y)), SCM_FRACTION_NUMERATOR (y)); else - SCM_WTA_DISPATCH_2 (g_divide, x, y, SCM_ARGn, s_divide); + return scm_wta_dispatch_2 (g_divide, x, y, SCM_ARGn, s_divide); } else if (SCM_REALP (x)) { @@ -6491,7 +8064,7 @@ do_divide (SCM x, SCM y, int inexact) else if (SCM_FRACTIONP (y)) return scm_from_double (rx / scm_i_fraction2double (y)); else - SCM_WTA_DISPATCH_2 (g_divide, x, y, SCM_ARGn, s_divide); + return scm_wta_dispatch_2 (g_divide, x, y, SCM_ARGn, s_divide); } else if (SCM_COMPLEXP (x)) { @@ -6549,7 +8122,7 @@ do_divide (SCM x, SCM y, int inexact) return scm_c_make_rectangular (rx / yy, ix / yy); } else - SCM_WTA_DISPATCH_2 (g_divide, x, y, SCM_ARGn, s_divide); + return scm_wta_dispatch_2 (g_divide, x, y, SCM_ARGn, s_divide); } else if (SCM_FRACTIONP (x)) { @@ -6588,10 +8161,10 @@ do_divide (SCM x, SCM y, int inexact) return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x), SCM_FRACTION_DENOMINATOR (y)), scm_product (SCM_FRACTION_NUMERATOR (y), SCM_FRACTION_DENOMINATOR (x))); else - SCM_WTA_DISPATCH_2 (g_divide, x, y, SCM_ARGn, s_divide); + return scm_wta_dispatch_2 (g_divide, x, y, SCM_ARGn, s_divide); } else - SCM_WTA_DISPATCH_2 (g_divide, x, y, SCM_ARG1, s_divide); + return scm_wta_dispatch_2 (g_divide, x, y, SCM_ARG1, s_divide); } SCM @@ -6610,13 +8183,7 @@ static SCM scm_divide2real (SCM x, SCM y) double scm_c_truncate (double x) { -#if HAVE_TRUNC return trunc (x); -#else - if (x < 0.0) - return -floor (-x); - return floor (x); -#endif } /* scm_c_round is done using floor(x+0.5) to round to nearest and with @@ -6661,43 +8228,41 @@ scm_c_round (double x) : result); } -SCM_DEFINE (scm_truncate_number, "truncate", 1, 0, 0, - (SCM x), - "Round the number @var{x} towards zero.") +SCM_PRIMITIVE_GENERIC (scm_truncate_number, "truncate", 1, 0, 0, + (SCM x), + "Round the number @var{x} towards zero.") #define FUNC_NAME s_scm_truncate_number { - if (scm_is_false (scm_negative_p (x))) - return scm_floor (x); + if (SCM_I_INUMP (x) || SCM_BIGP (x)) + return x; + else if (SCM_REALP (x)) + return scm_from_double (trunc (SCM_REAL_VALUE (x))); + else if (SCM_FRACTIONP (x)) + return scm_truncate_quotient (SCM_FRACTION_NUMERATOR (x), + SCM_FRACTION_DENOMINATOR (x)); else - return scm_ceiling (x); + return scm_wta_dispatch_1 (g_scm_truncate_number, x, SCM_ARG1, + s_scm_truncate_number); } #undef FUNC_NAME -SCM_DEFINE (scm_round_number, "round", 1, 0, 0, - (SCM x), - "Round the number @var{x} towards the nearest integer. " - "When it is exactly halfway between two integers, " - "round towards the even one.") +SCM_PRIMITIVE_GENERIC (scm_round_number, "round", 1, 0, 0, + (SCM x), + "Round the number @var{x} towards the nearest integer. " + "When it is exactly halfway between two integers, " + "round towards the even one.") #define FUNC_NAME s_scm_round_number { if (SCM_I_INUMP (x) || SCM_BIGP (x)) return x; else if (SCM_REALP (x)) return scm_from_double (scm_c_round (SCM_REAL_VALUE (x))); + else if (SCM_FRACTIONP (x)) + return scm_round_quotient (SCM_FRACTION_NUMERATOR (x), + SCM_FRACTION_DENOMINATOR (x)); else - { - /* OPTIMIZE-ME: Fraction case could be done more efficiently by a - single quotient+remainder division then examining to see which way - the rounding should go. */ - SCM plus_half = scm_sum (x, exactly_one_half); - SCM result = scm_floor (plus_half); - /* Adjust so that the rounding is towards even. */ - if (scm_is_true (scm_num_eq_p (plus_half, result)) - && scm_is_true (scm_odd_p (result))) - return scm_difference (result, SCM_INUM1); - else - return result; - } + return scm_wta_dispatch_1 (g_scm_round_number, x, SCM_ARG1, + s_scm_round_number); } #undef FUNC_NAME @@ -6711,24 +8276,10 @@ SCM_PRIMITIVE_GENERIC (scm_floor, "floor", 1, 0, 0, else if (SCM_REALP (x)) return scm_from_double (floor (SCM_REAL_VALUE (x))); else if (SCM_FRACTIONP (x)) - { - SCM q = scm_quotient (SCM_FRACTION_NUMERATOR (x), - SCM_FRACTION_DENOMINATOR (x)); - if (scm_is_false (scm_negative_p (x))) - { - /* For positive x, rounding towards zero is correct. */ - return q; - } - else - { - /* For negative x, we need to return q-1 unless x is an - integer. But fractions are never integer, per our - assumptions. */ - return scm_difference (q, SCM_INUM1); - } - } + return scm_floor_quotient (SCM_FRACTION_NUMERATOR (x), + SCM_FRACTION_DENOMINATOR (x)); else - SCM_WTA_DISPATCH_1 (g_scm_floor, x, 1, s_scm_floor); + return scm_wta_dispatch_1 (g_scm_floor, x, 1, s_scm_floor); } #undef FUNC_NAME @@ -6742,24 +8293,10 @@ SCM_PRIMITIVE_GENERIC (scm_ceiling, "ceiling", 1, 0, 0, else if (SCM_REALP (x)) return scm_from_double (ceil (SCM_REAL_VALUE (x))); else if (SCM_FRACTIONP (x)) - { - SCM q = scm_quotient (SCM_FRACTION_NUMERATOR (x), - SCM_FRACTION_DENOMINATOR (x)); - if (scm_is_false (scm_positive_p (x))) - { - /* For negative x, rounding towards zero is correct. */ - return q; - } - else - { - /* For positive x, we need to return q+1 unless x is an - integer. But fractions are never integer, per our - assumptions. */ - return scm_sum (q, SCM_INUM1); - } - } + return scm_ceiling_quotient (SCM_FRACTION_NUMERATOR (x), + SCM_FRACTION_DENOMINATOR (x)); else - SCM_WTA_DISPATCH_1 (g_scm_ceiling, x, 1, s_scm_ceiling); + return scm_wta_dispatch_1 (g_scm_ceiling, x, 1, s_scm_ceiling); } #undef FUNC_NAME @@ -6798,9 +8335,9 @@ SCM_PRIMITIVE_GENERIC (scm_expt, "expt", 2, 0, 0, else if (scm_is_complex (x) && scm_is_complex (y)) return scm_exp (scm_product (scm_log (x), y)); else if (scm_is_complex (x)) - SCM_WTA_DISPATCH_2 (g_scm_expt, x, y, SCM_ARG2, s_scm_expt); + return scm_wta_dispatch_2 (g_scm_expt, x, y, SCM_ARG2, s_scm_expt); else - SCM_WTA_DISPATCH_2 (g_scm_expt, x, y, SCM_ARG1, s_scm_expt); + return scm_wta_dispatch_2 (g_scm_expt, x, y, SCM_ARG1, s_scm_expt); } #undef FUNC_NAME @@ -6827,7 +8364,7 @@ SCM_PRIMITIVE_GENERIC (scm_sin, "sin", 1, 0, 0, cos (x) * sinh (y)); } else - SCM_WTA_DISPATCH_1 (g_scm_sin, z, 1, s_scm_sin); + return scm_wta_dispatch_1 (g_scm_sin, z, 1, s_scm_sin); } #undef FUNC_NAME @@ -6848,7 +8385,7 @@ SCM_PRIMITIVE_GENERIC (scm_cos, "cos", 1, 0, 0, -sin (x) * sinh (y)); } else - SCM_WTA_DISPATCH_1 (g_scm_cos, z, 1, s_scm_cos); + return scm_wta_dispatch_1 (g_scm_cos, z, 1, s_scm_cos); } #undef FUNC_NAME @@ -6873,7 +8410,7 @@ SCM_PRIMITIVE_GENERIC (scm_tan, "tan", 1, 0, 0, return scm_c_make_rectangular (sin (x) / w, sinh (y) / w); } else - SCM_WTA_DISPATCH_1 (g_scm_tan, z, 1, s_scm_tan); + return scm_wta_dispatch_1 (g_scm_tan, z, 1, s_scm_tan); } #undef FUNC_NAME @@ -6894,7 +8431,7 @@ SCM_PRIMITIVE_GENERIC (scm_sinh, "sinh", 1, 0, 0, cosh (x) * sin (y)); } else - SCM_WTA_DISPATCH_1 (g_scm_sinh, z, 1, s_scm_sinh); + return scm_wta_dispatch_1 (g_scm_sinh, z, 1, s_scm_sinh); } #undef FUNC_NAME @@ -6915,7 +8452,7 @@ SCM_PRIMITIVE_GENERIC (scm_cosh, "cosh", 1, 0, 0, sinh (x) * sin (y)); } else - SCM_WTA_DISPATCH_1 (g_scm_cosh, z, 1, s_scm_cosh); + return scm_wta_dispatch_1 (g_scm_cosh, z, 1, s_scm_cosh); } #undef FUNC_NAME @@ -6940,7 +8477,7 @@ SCM_PRIMITIVE_GENERIC (scm_tanh, "tanh", 1, 0, 0, return scm_c_make_rectangular (sinh (x) / w, sin (y) / w); } else - SCM_WTA_DISPATCH_1 (g_scm_tanh, z, 1, s_scm_tanh); + return scm_wta_dispatch_1 (g_scm_tanh, z, 1, s_scm_tanh); } #undef FUNC_NAME @@ -6968,7 +8505,7 @@ SCM_PRIMITIVE_GENERIC (scm_asin, "asin", 1, 0, 0, scm_sys_asinh (scm_c_make_rectangular (-y, x))); } else - SCM_WTA_DISPATCH_1 (g_scm_asin, z, 1, s_scm_asin); + return scm_wta_dispatch_1 (g_scm_asin, z, 1, s_scm_asin); } #undef FUNC_NAME @@ -6998,7 +8535,7 @@ SCM_PRIMITIVE_GENERIC (scm_acos, "acos", 1, 0, 0, scm_sys_asinh (scm_c_make_rectangular (-y, x)))); } else - SCM_WTA_DISPATCH_1 (g_scm_acos, z, 1, s_scm_acos); + return scm_wta_dispatch_1 (g_scm_acos, z, 1, s_scm_acos); } #undef FUNC_NAME @@ -7025,17 +8562,17 @@ SCM_PRIMITIVE_GENERIC (scm_atan, "atan", 1, 1, 0, scm_c_make_rectangular (0, 2)); } else - SCM_WTA_DISPATCH_2 (g_scm_atan, z, y, SCM_ARG1, s_scm_atan); + return scm_wta_dispatch_1 (g_scm_atan, z, SCM_ARG1, s_scm_atan); } else if (scm_is_real (z)) { if (scm_is_real (y)) return scm_from_double (atan2 (scm_to_double (z), scm_to_double (y))); else - SCM_WTA_DISPATCH_2 (g_scm_atan, z, y, SCM_ARG2, s_scm_atan); + return scm_wta_dispatch_2 (g_scm_atan, z, y, SCM_ARG2, s_scm_atan); } else - SCM_WTA_DISPATCH_2 (g_scm_atan, z, y, SCM_ARG1, s_scm_atan); + return scm_wta_dispatch_2 (g_scm_atan, z, y, SCM_ARG1, s_scm_atan); } #undef FUNC_NAME @@ -7053,7 +8590,7 @@ SCM_PRIMITIVE_GENERIC (scm_sys_asinh, "asinh", 1, 0, 0, scm_sqrt (scm_sum (scm_product (z, z), SCM_INUM1)))); else - SCM_WTA_DISPATCH_1 (g_scm_sys_asinh, z, 1, s_scm_sys_asinh); + return scm_wta_dispatch_1 (g_scm_sys_asinh, z, 1, s_scm_sys_asinh); } #undef FUNC_NAME @@ -7071,7 +8608,7 @@ SCM_PRIMITIVE_GENERIC (scm_sys_acosh, "acosh", 1, 0, 0, scm_sqrt (scm_difference (scm_product (z, z), SCM_INUM1)))); else - SCM_WTA_DISPATCH_1 (g_scm_sys_acosh, z, 1, s_scm_sys_acosh); + return scm_wta_dispatch_1 (g_scm_sys_acosh, z, 1, s_scm_sys_acosh); } #undef FUNC_NAME @@ -7089,7 +8626,7 @@ SCM_PRIMITIVE_GENERIC (scm_sys_atanh, "atanh", 1, 0, 0, scm_difference (SCM_INUM1, z))), SCM_I_MAKINUM (2)); else - SCM_WTA_DISPATCH_1 (g_scm_sys_atanh, z, 1, s_scm_sys_atanh); + return scm_wta_dispatch_1 (g_scm_sys_atanh, z, 1, s_scm_sys_atanh); } #undef FUNC_NAME @@ -7190,7 +8727,7 @@ SCM_PRIMITIVE_GENERIC (scm_real_part, "real-part", 1, 0, 0, else if (SCM_I_INUMP (z) || SCM_BIGP (z) || SCM_REALP (z) || SCM_FRACTIONP (z)) return z; else - SCM_WTA_DISPATCH_1 (g_scm_real_part, z, SCM_ARG1, s_scm_real_part); + return scm_wta_dispatch_1 (g_scm_real_part, z, SCM_ARG1, s_scm_real_part); } #undef FUNC_NAME @@ -7205,7 +8742,7 @@ SCM_PRIMITIVE_GENERIC (scm_imag_part, "imag-part", 1, 0, 0, else if (SCM_I_INUMP (z) || SCM_REALP (z) || SCM_BIGP (z) || SCM_FRACTIONP (z)) return SCM_INUM0; else - SCM_WTA_DISPATCH_1 (g_scm_imag_part, z, SCM_ARG1, s_scm_imag_part); + return scm_wta_dispatch_1 (g_scm_imag_part, z, SCM_ARG1, s_scm_imag_part); } #undef FUNC_NAME @@ -7221,7 +8758,7 @@ SCM_PRIMITIVE_GENERIC (scm_numerator, "numerator", 1, 0, 0, else if (SCM_REALP (z)) return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z))); else - SCM_WTA_DISPATCH_1 (g_scm_numerator, z, SCM_ARG1, s_scm_numerator); + return scm_wta_dispatch_1 (g_scm_numerator, z, SCM_ARG1, s_scm_numerator); } #undef FUNC_NAME @@ -7238,7 +8775,8 @@ SCM_PRIMITIVE_GENERIC (scm_denominator, "denominator", 1, 0, 0, else if (SCM_REALP (z)) return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z))); else - SCM_WTA_DISPATCH_1 (g_scm_denominator, z, SCM_ARG1, s_scm_denominator); + return scm_wta_dispatch_1 (g_scm_denominator, z, SCM_ARG1, + s_scm_denominator); } #undef FUNC_NAME @@ -7280,7 +8818,8 @@ SCM_PRIMITIVE_GENERIC (scm_magnitude, "magnitude", 1, 0, 0, SCM_FRACTION_DENOMINATOR (z)); } else - SCM_WTA_DISPATCH_1 (g_scm_magnitude, z, SCM_ARG1, s_scm_magnitude); + return scm_wta_dispatch_1 (g_scm_magnitude, z, SCM_ARG1, + s_scm_magnitude); } #undef FUNC_NAME @@ -7326,7 +8865,7 @@ SCM_PRIMITIVE_GENERIC (scm_angle, "angle", 1, 0, 0, else return scm_from_double (atan2 (0.0, -1.0)); } else - SCM_WTA_DISPATCH_1 (g_scm_angle, z, SCM_ARG1, s_scm_angle); + return scm_wta_dispatch_1 (g_scm_angle, z, SCM_ARG1, s_scm_angle); } #undef FUNC_NAME @@ -7345,7 +8884,8 @@ SCM_PRIMITIVE_GENERIC (scm_exact_to_inexact, "exact->inexact", 1, 0, 0, else if (SCM_INEXACTP (z)) return z; else - SCM_WTA_DISPATCH_1 (g_scm_exact_to_inexact, z, 1, s_scm_exact_to_inexact); + return scm_wta_dispatch_1 (g_scm_exact_to_inexact, z, 1, + s_scm_exact_to_inexact); } #undef FUNC_NAME @@ -7366,7 +8906,8 @@ SCM_PRIMITIVE_GENERIC (scm_inexact_to_exact, "inexact->exact", 1, 0, 0, else if (SCM_COMPLEXP (z) && SCM_COMPLEX_IMAG (z) == 0.0) val = SCM_COMPLEX_REAL (z); else - SCM_WTA_DISPATCH_1 (g_scm_inexact_to_exact, z, 1, s_scm_inexact_to_exact); + return scm_wta_dispatch_1 (g_scm_inexact_to_exact, z, 1, + s_scm_inexact_to_exact); if (!SCM_LIKELY (DOUBLE_IS_FINITE (val))) SCM_OUT_OF_RANGE (1, z); @@ -7754,46 +9295,6 @@ scm_from_double (double val) return z; } -#if SCM_ENABLE_DEPRECATED == 1 - -float -scm_num2float (SCM num, unsigned long pos, const char *s_caller) -{ - scm_c_issue_deprecation_warning - ("`scm_num2float' is deprecated. Use scm_to_double instead."); - - if (SCM_BIGP (num)) - { - float res = mpz_get_d (SCM_I_BIG_MPZ (num)); - if (!isinf (res)) - return res; - else - scm_out_of_range (NULL, num); - } - else - return scm_to_double (num); -} - -double -scm_num2double (SCM num, unsigned long pos, const char *s_caller) -{ - scm_c_issue_deprecation_warning - ("`scm_num2double' is deprecated. Use scm_to_double instead."); - - if (SCM_BIGP (num)) - { - double res = mpz_get_d (SCM_I_BIG_MPZ (num)); - if (!isinf (res)) - return res; - else - scm_out_of_range (NULL, num); - } - else - return scm_to_double (num); -} - -#endif - int scm_is_complex (SCM val) { @@ -7848,6 +9349,62 @@ scm_is_number (SCM z) } +/* Returns log(x * 2^shift) */ +static SCM +log_of_shifted_double (double x, long shift) +{ + double ans = log (fabs (x)) + shift * M_LN2; + + if (x > 0.0 || double_is_non_negative_zero (x)) + return scm_from_double (ans); + else + return scm_c_make_rectangular (ans, M_PI); +} + +/* Returns log(n), for exact integer n of integer-length size */ +static SCM +log_of_exact_integer_with_size (SCM n, long size) +{ + long shift = size - 2 * scm_dblprec[0]; + + if (shift > 0) + return log_of_shifted_double + (scm_to_double (scm_ash (n, scm_from_long(-shift))), + shift); + else + return log_of_shifted_double (scm_to_double (n), 0); +} + +/* Returns log(n), for exact integer n */ +static SCM +log_of_exact_integer (SCM n) +{ + return log_of_exact_integer_with_size + (n, scm_to_long (scm_integer_length (n))); +} + +/* Returns log(n/d), for exact non-zero integers n and d */ +static SCM +log_of_fraction (SCM n, SCM d) +{ + long n_size = scm_to_long (scm_integer_length (n)); + long d_size = scm_to_long (scm_integer_length (d)); + + if (abs (n_size - d_size) > 1) + return (scm_difference (log_of_exact_integer_with_size (n, n_size), + log_of_exact_integer_with_size (d, d_size))); + else if (scm_is_false (scm_negative_p (n))) + return scm_from_double + (log1p (scm_to_double (scm_divide2real (scm_difference (n, d), d)))); + else + return scm_c_make_rectangular + (log1p (scm_to_double (scm_divide2real + (scm_difference (scm_abs (n), d), + d))), + M_PI); +} + + /* In the following functions we dispatch to the real-arg funcs like log() when we know the arg is real, instead of just handing everything to clog() for instance. This is in case clog() doesn't optimize for a @@ -7861,7 +9418,8 @@ SCM_PRIMITIVE_GENERIC (scm_log, "log", 1, 0, 0, { if (SCM_COMPLEXP (z)) { -#if HAVE_COMPLEX_DOUBLE && HAVE_CLOG && defined (SCM_COMPLEX_VALUE) +#if defined HAVE_COMPLEX_DOUBLE && defined HAVE_CLOG \ + && defined (SCM_COMPLEX_VALUE) return scm_from_complex_double (clog (SCM_COMPLEX_VALUE (z))); #else double re = SCM_COMPLEX_REAL (z); @@ -7870,19 +9428,23 @@ SCM_PRIMITIVE_GENERIC (scm_log, "log", 1, 0, 0, atan2 (im, re)); #endif } - else if (SCM_NUMBERP (z)) + else if (SCM_REALP (z)) + return log_of_shifted_double (SCM_REAL_VALUE (z), 0); + else if (SCM_I_INUMP (z)) { - /* ENHANCE-ME: When z is a bignum the logarithm will fit a double - although the value itself overflows. */ - double re = scm_to_double (z); - double l = log (fabs (re)); - if (re >= 0.0) - return scm_from_double (l); - else - return scm_c_make_rectangular (l, M_PI); +#ifndef ALLOW_DIVIDE_BY_EXACT_ZERO + if (scm_is_eq (z, SCM_INUM0)) + scm_num_overflow (s_scm_log); +#endif + return log_of_shifted_double (SCM_I_INUM (z), 0); } + else if (SCM_BIGP (z)) + return log_of_exact_integer (z); + else if (SCM_FRACTIONP (z)) + return log_of_fraction (SCM_FRACTION_NUMERATOR (z), + SCM_FRACTION_DENOMINATOR (z)); else - SCM_WTA_DISPATCH_1 (g_scm_log, z, 1, s_scm_log); + return scm_wta_dispatch_1 (g_scm_log, z, 1, s_scm_log); } #undef FUNC_NAME @@ -7907,19 +9469,29 @@ SCM_PRIMITIVE_GENERIC (scm_log10, "log10", 1, 0, 0, M_LOG10E * atan2 (im, re)); #endif } - else if (SCM_NUMBERP (z)) + else if (SCM_REALP (z) || SCM_I_INUMP (z)) { - /* ENHANCE-ME: When z is a bignum the logarithm will fit a double - although the value itself overflows. */ - double re = scm_to_double (z); - double l = log10 (fabs (re)); - if (re >= 0.0) - return scm_from_double (l); - else - return scm_c_make_rectangular (l, M_LOG10E * M_PI); +#ifndef ALLOW_DIVIDE_BY_EXACT_ZERO + if (scm_is_eq (z, SCM_INUM0)) + scm_num_overflow (s_scm_log10); +#endif + { + double re = scm_to_double (z); + double l = log10 (fabs (re)); + if (re > 0.0 || double_is_non_negative_zero (re)) + return scm_from_double (l); + else + return scm_c_make_rectangular (l, M_LOG10E * M_PI); + } } + else if (SCM_BIGP (z)) + return scm_product (flo_log10e, log_of_exact_integer (z)); + else if (SCM_FRACTIONP (z)) + return scm_product (flo_log10e, + log_of_fraction (SCM_FRACTION_NUMERATOR (z), + SCM_FRACTION_DENOMINATOR (z))); else - SCM_WTA_DISPATCH_1 (g_scm_log10, z, 1, s_scm_log10); + return scm_wta_dispatch_1 (g_scm_log10, z, 1, s_scm_log10); } #undef FUNC_NAME @@ -7932,7 +9504,8 @@ SCM_PRIMITIVE_GENERIC (scm_exp, "exp", 1, 0, 0, { if (SCM_COMPLEXP (z)) { -#if HAVE_COMPLEX_DOUBLE && HAVE_CEXP && defined (SCM_COMPLEX_VALUE) +#if defined HAVE_COMPLEX_DOUBLE && defined HAVE_CEXP \ + && defined (SCM_COMPLEX_VALUE) return scm_from_complex_double (cexp (SCM_COMPLEX_VALUE (z))); #else return scm_c_make_polar (exp (SCM_COMPLEX_REAL (z)), @@ -7946,15 +9519,79 @@ SCM_PRIMITIVE_GENERIC (scm_exp, "exp", 1, 0, 0, return scm_from_double (exp (scm_to_double (z))); } else - SCM_WTA_DISPATCH_1 (g_scm_exp, z, 1, s_scm_exp); + return scm_wta_dispatch_1 (g_scm_exp, z, 1, s_scm_exp); +} +#undef FUNC_NAME + + +SCM_DEFINE (scm_i_exact_integer_sqrt, "exact-integer-sqrt", 1, 0, 0, + (SCM k), + "Return two exact non-negative integers @var{s} and @var{r}\n" + "such that @math{@var{k} = @var{s}^2 + @var{r}} and\n" + "@math{@var{s}^2 <= @var{k} < (@var{s} + 1)^2}.\n" + "An error is raised if @var{k} is not an exact non-negative integer.\n" + "\n" + "@lisp\n" + "(exact-integer-sqrt 10) @result{} 3 and 1\n" + "@end lisp") +#define FUNC_NAME s_scm_i_exact_integer_sqrt +{ + SCM s, r; + + scm_exact_integer_sqrt (k, &s, &r); + return scm_values (scm_list_2 (s, r)); } #undef FUNC_NAME +void +scm_exact_integer_sqrt (SCM k, SCM *sp, SCM *rp) +{ + if (SCM_LIKELY (SCM_I_INUMP (k))) + { + scm_t_inum kk = SCM_I_INUM (k); + scm_t_inum uu = kk; + scm_t_inum ss; + + if (SCM_LIKELY (kk > 0)) + { + do + { + ss = uu; + uu = (ss + kk/ss) / 2; + } while (uu < ss); + *sp = SCM_I_MAKINUM (ss); + *rp = SCM_I_MAKINUM (kk - ss*ss); + } + else if (SCM_LIKELY (kk == 0)) + *sp = *rp = SCM_INUM0; + else + scm_wrong_type_arg_msg ("exact-integer-sqrt", SCM_ARG1, k, + "exact non-negative integer"); + } + else if (SCM_LIKELY (SCM_BIGP (k))) + { + SCM s, r; + + if (mpz_sgn (SCM_I_BIG_MPZ (k)) < 0) + scm_wrong_type_arg_msg ("exact-integer-sqrt", SCM_ARG1, k, + "exact non-negative integer"); + s = scm_i_mkbig (); + r = scm_i_mkbig (); + mpz_sqrtrem (SCM_I_BIG_MPZ (s), SCM_I_BIG_MPZ (r), SCM_I_BIG_MPZ (k)); + scm_remember_upto_here_1 (k); + *sp = scm_i_normbig (s); + *rp = scm_i_normbig (r); + } + else + scm_wrong_type_arg_msg ("exact-integer-sqrt", SCM_ARG1, k, + "exact non-negative integer"); +} + SCM_PRIMITIVE_GENERIC (scm_sqrt, "sqrt", 1, 0, 0, (SCM z), "Return the square root of @var{z}. Of the two possible roots\n" - "(positive and negative), the one with the a positive real part\n" + "(positive and negative), the one with positive real part\n" "is returned, or if that's zero then a positive imaginary part.\n" "Thus,\n" "\n" @@ -7987,7 +9624,7 @@ SCM_PRIMITIVE_GENERIC (scm_sqrt, "sqrt", 1, 0, 0, return scm_from_double (sqrt (xx)); } else - SCM_WTA_DISPATCH_1 (g_scm_sqrt, z, 1, s_scm_sqrt); + return scm_wta_dispatch_1 (g_scm_sqrt, z, 1, s_scm_sqrt); } #undef FUNC_NAME @@ -8012,6 +9649,7 @@ scm_init_numbers () scm_add_feature ("complex"); scm_add_feature ("inexact"); flo0 = scm_from_double (0.0); + flo_log10e = scm_from_double (M_LOG10E); /* determine floating point precision */ for (i=2; i <= SCM_MAX_DBL_RADIX; ++i)