From 2d5dc6a14c4e503105c5805bafc6699ad202ac17 Mon Sep 17 00:00:00 2001 From: Andy Wingo Date: Wed, 29 Dec 2021 20:39:11 +0100 Subject: [PATCH] Implement scm_modulo_expt with new integer library * libguile/integers.c (scm_integer_modulo_expt_nnn): (integer_init_mpz): New helper. * libguile/integers.h: Declare the new internal function. * libguile/numbers.c (scm_modulo_expt): Use new internal function. --- libguile/integers.c | 57 +++++++++++++++++++++++ libguile/integers.h | 2 + libguile/numbers.c | 110 +++----------------------------------------- 3 files changed, 66 insertions(+), 103 deletions(-) diff --git a/libguile/integers.c b/libguile/integers.c index 2ae2c30d5..a20fdf7db 100644 --- a/libguile/integers.c +++ b/libguile/integers.c @@ -2045,3 +2045,60 @@ scm_integer_lognot_z (SCM n) scm_remember_upto_here_1 (n); return take_mpz (result); } + +static void +integer_init_mpz (mpz_ptr z, SCM n) +{ + if (SCM_I_INUMP (n)) + mpz_init_set_si (z, SCM_I_INUM (n)); + else + { + ASSERT (SCM_BIGP (n)); + mpz_t zn; + alias_bignum_to_mpz (scm_bignum (n), zn); + mpz_init_set (z, zn); + scm_remember_upto_here_1 (n); + } +} + +SCM +scm_integer_modulo_expt_nnn (SCM n, SCM k, SCM m) +{ + if (scm_is_eq (m, SCM_INUM0)) + scm_num_overflow ("modulo-expt"); + + mpz_t n_tmp, k_tmp, m_tmp; + + integer_init_mpz (n_tmp, n); + integer_init_mpz (k_tmp, k); + integer_init_mpz (m_tmp, m); + + /* if the exponent K is negative, and we simply call mpz_powm, we + will get a divide-by-zero exception when an inverse 1/n mod m + doesn't exist (or is not unique). Since exceptions are hard to + handle, we'll attempt the inversion "by hand" -- that way, we get + a simple failure code, which is easy to handle. */ + + if (-1 == mpz_sgn (k_tmp)) + { + if (!mpz_invert (n_tmp, n_tmp, m_tmp)) + { + mpz_clear (n_tmp); + mpz_clear (k_tmp); + mpz_clear (m_tmp); + + scm_num_overflow ("modulo-expt"); + } + mpz_neg (k_tmp, k_tmp); + } + + mpz_powm (n_tmp, n_tmp, k_tmp, m_tmp); + + if (mpz_sgn (m_tmp) < 0 && mpz_sgn (n_tmp) != 0) + mpz_add (n_tmp, n_tmp, m_tmp); + + mpz_clear (m_tmp); + mpz_clear (k_tmp); + + return take_mpz (n_tmp); +} diff --git a/libguile/integers.h b/libguile/integers.h index 105b86b63..74624f143 100644 --- a/libguile/integers.h +++ b/libguile/integers.h @@ -154,6 +154,8 @@ SCM_INTERNAL int scm_integer_logbit_uz (unsigned long bit, SCM n); SCM_INTERNAL SCM scm_integer_lognot_i (scm_t_inum n); SCM_INTERNAL SCM scm_integer_lognot_z (SCM n); +SCM_INTERNAL SCM scm_integer_modulo_expt_nnn (SCM n, SCM k, SCM m); + #endif /* SCM_INTEGERS_H */ diff --git a/libguile/numbers.c b/libguile/numbers.c index 3e8431757..0afa83b79 100644 --- a/libguile/numbers.c +++ b/libguile/numbers.c @@ -3186,21 +3186,6 @@ SCM_DEFINE (scm_lognot, "lognot", 1, 0, 0, } #undef FUNC_NAME -/* returns 0 if IN is not an integer. OUT must already be - initialized. */ -static int -coerce_to_big (SCM in, mpz_t out) -{ - if (SCM_BIGP (in)) - mpz_set (out, SCM_I_BIG_MPZ (in)); - else if (SCM_I_INUMP (in)) - mpz_set_si (out, SCM_I_INUM (in)); - else - return 0; - - return 1; -} - SCM_DEFINE (scm_modulo_expt, "modulo-expt", 3, 0, 0, (SCM n, SCM k, SCM m), "Return @var{n} raised to the integer exponent\n" @@ -3212,95 +3197,14 @@ SCM_DEFINE (scm_modulo_expt, "modulo-expt", 3, 0, 0, "@end lisp") #define FUNC_NAME s_scm_modulo_expt { - mpz_t n_tmp; - mpz_t k_tmp; - mpz_t m_tmp; - - /* There are two classes of error we might encounter -- - 1) Math errors, which we'll report by calling scm_num_overflow, - and - 2) wrong-type errors, which of course we'll report by calling - SCM_WRONG_TYPE_ARG. - We don't report those errors immediately, however; instead we do - some cleanup first. These variables tell us which error (if - any) we should report after cleaning up. - */ - int report_overflow = 0; + if (!scm_is_exact_integer (n)) + SCM_WRONG_TYPE_ARG (SCM_ARG1, n); + if (!scm_is_exact_integer (k)) + SCM_WRONG_TYPE_ARG (SCM_ARG2, k); + if (!scm_is_exact_integer (m)) + SCM_WRONG_TYPE_ARG (SCM_ARG3, m); - int position_of_wrong_type = 0; - SCM value_of_wrong_type = SCM_INUM0; - - SCM result = SCM_UNDEFINED; - - mpz_init (n_tmp); - mpz_init (k_tmp); - mpz_init (m_tmp); - - if (scm_is_eq (m, SCM_INUM0)) - { - report_overflow = 1; - goto cleanup; - } - - if (!coerce_to_big (n, n_tmp)) - { - value_of_wrong_type = n; - position_of_wrong_type = 1; - goto cleanup; - } - - if (!coerce_to_big (k, k_tmp)) - { - value_of_wrong_type = k; - position_of_wrong_type = 2; - goto cleanup; - } - - if (!coerce_to_big (m, m_tmp)) - { - value_of_wrong_type = m; - position_of_wrong_type = 3; - goto cleanup; - } - - /* if the exponent K is negative, and we simply call mpz_powm, we - will get a divide-by-zero exception when an inverse 1/n mod m - doesn't exist (or is not unique). Since exceptions are hard to - handle, we'll attempt the inversion "by hand" -- that way, we get - a simple failure code, which is easy to handle. */ - - if (-1 == mpz_sgn (k_tmp)) - { - if (!mpz_invert (n_tmp, n_tmp, m_tmp)) - { - report_overflow = 1; - goto cleanup; - } - mpz_neg (k_tmp, k_tmp); - } - - result = scm_i_mkbig (); - mpz_powm (SCM_I_BIG_MPZ (result), - n_tmp, - k_tmp, - m_tmp); - - if (mpz_sgn (m_tmp) < 0 && mpz_sgn (SCM_I_BIG_MPZ (result)) != 0) - mpz_add (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (result), m_tmp); - - cleanup: - mpz_clear (m_tmp); - mpz_clear (k_tmp); - mpz_clear (n_tmp); - - if (report_overflow) - scm_num_overflow (FUNC_NAME); - - if (position_of_wrong_type) - SCM_WRONG_TYPE_ARG (position_of_wrong_type, - value_of_wrong_type); - - return scm_i_normbig (result); + return scm_integer_modulo_expt_nnn (n, k, m); } #undef FUNC_NAME