mirror/guile
mirror/guile
1
Fork 0
mirror of https://git.savannah.gnu.org/git/guile.git synced 2025-05-20 11:40:18 +02:00
Code Activity

Implement scm_modulo_expt with new integer library

Browse source 
* 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.
This commit is contained in:
Andy Wingo 2021-12-29 20:39:11 +01:00
parent b41714d277
commit 2d5dc6a14c
3 changed files with 66 additions and 103 deletions
Download patch file Download diff file Expand all files Collapse all files

57
libguile/integers.c
View file

@ -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);
}

2
libguile/integers.h
View file

@ -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 */

110
libguile/numbers.c
View file

@ -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