summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorAndy Wingo <wingo@pobox.com>2021-12-29 20:39:11 +0100
committerAndy Wingo <wingo@pobox.com>2022-01-13 09:37:16 +0100
commit2d5dc6a14c4e503105c5805bafc6699ad202ac17 (patch)
treef297614a82899bb5cc906fe504b30a4b0b463677
parentb41714d277b8d6629e358344b81d55b7ca4ef377 (diff)
downloadguile-2d5dc6a14c4e503105c5805bafc6699ad202ac17.tar.gz
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.
-rw-r--r--libguile/integers.c57
-rw-r--r--libguile/integers.h2
-rw-r--r--libguile/numbers.c110
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;
-
- 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 (!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);
- 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