[Guile-commits] 28/69: Implement scm_modulo_expt with new integer library

[Andy Wingo]

wingo pushed a commit to branch wip-inline-digits
in repository guile.

commit f9825d986350548c3356b535372dce310b6cd3b6
Author: Andy Wingo <wingo@pobox.com>
AuthorDate: Wed Dec 29 20:39:11 2021 +0100

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