Re: Forwarded patch for modular exponentiation support (GMP powm)

[Eric Hanchrow <offby1@blarg.net>]

[-- Attachment #1: Type: text/plain, Size: 115 bytes --]

Latest version.  I think I've incorporated every suggestion that
Kevin has made (thanks for your patience, Kevin).

[-- Attachment #2: modulo-expt path --]
[-- Type: application/octet-stream, Size: 4903 bytes --]

? libguile/guile.info
? libguile/little-test.scm
? test-suite/tests/hmm.pl
? test-suite/tests/mexp.test
Index: libguile/numbers.c
===================================================================
RCS file: /cvsroot/guile/guile/guile-core/libguile/numbers.c,v
retrieving revision 1.221
diff -w -u -r1.221 numbers.c
--- libguile/numbers.c	6 Jan 2004 21:55:29 -0000	1.221
+++ libguile/numbers.c	12 Feb 2004 04:52:26 -0000
@@ -1519,6 +1519,76 @@
 }
 #undef FUNC_NAME
 
+static void
+coerce_to_big (SCM in, mpz_t out)
+{
+  if (SCM_BIGP (in))
+    mpz_init_set (out, SCM_I_BIG_MPZ (in));
+  else if (SCM_INUMP (in))
+    mpz_init_set_si (out, SCM_INUM (in));
+  else
+    scm_wrong_type_arg ("modulo-expt", 1, in);
+}
+
+SCM_DEFINE (scm_modular_expt, "modulo-expt", 3, 0, 0,
+            (SCM n, SCM k, SCM m),
+            "Return @var{n} raised to the integer exponent\n"
+	    "@var{k}, modulo @var{m}.\n"
+	    "\n"
+	    "@lisp\n"
+	    "(modulo-expt 2 3 5)\n"
+	    "   @result{} 3\n"
+	    "@end lisp")
+#define FUNC_NAME s_scm_modular_expt
+{
+  mpz_t n_tmp; 
+  mpz_t k_tmp; 
+  mpz_t m_tmp; 
+
+  int report_overflow = 0;
+
+  SCM result = SCM_UNDEFINED;
+  
+  if (SCM_NFALSEP (scm_zero_p (m)))
+    scm_num_overflow (FUNC_NAME);
+  
+  coerce_to_big (n, n_tmp);
+  coerce_to_big (k, k_tmp);
+  coerce_to_big (m, m_tmp);
+  
+  /* 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);
+ cleanup:
+  mpz_clear (m_tmp);
+  mpz_clear (k_tmp);
+  mpz_clear (n_tmp);
+
+  if (report_overflow)
+    scm_num_overflow (FUNC_NAME);
+
+  return scm_i_normbig (result);
+}
+#undef FUNC_NAME
+
 SCM_DEFINE (scm_integer_expt, "integer-expt", 2, 0, 0,
             (SCM n, SCM k),
 	    "Return @var{n} raised to the non-negative integer exponent\n"
Index: libguile/numbers.h
===================================================================
RCS file: /cvsroot/guile/guile/guile-core/libguile/numbers.h,v
retrieving revision 1.77
diff -w -u -r1.77 numbers.h
--- libguile/numbers.h	18 Nov 2003 19:59:51 -0000	1.77
+++ libguile/numbers.h	12 Feb 2004 04:52:26 -0000
@@ -201,6 +201,7 @@
 SCM_API SCM scm_logtest (SCM n1, SCM n2);
 SCM_API SCM scm_logbit_p (SCM n1, SCM n2);
 SCM_API SCM scm_lognot (SCM n);
+SCM_API SCM scm_modular_expt (SCM n, SCM k, SCM m);
 SCM_API SCM scm_integer_expt (SCM z1, SCM z2);
 SCM_API SCM scm_ash (SCM n, SCM cnt);
 SCM_API SCM scm_bit_extract (SCM n, SCM start, SCM end);
--- /dev/null	1969-12-31 16:00:00.000000000 -0800
+++ test-suite/tests/mexp.test	2004-02-10 16:36:58.000000000 -0800
@@ -0,0 +1,42 @@
+;;;; mexp.test --- test suite for Guile's modular exponentiation functions    -*- scheme -*-
+;;;; Eric Hanchrow <offby1@blarg.net> --- January 2004
+;;;;
+;;;; Copyright (C) 2004 Free Software Foundation, Inc.
+;;;; 
+;;;; This program is free software; you can redistribute it and/or modify
+;;;; it under the terms of the GNU General Public License as published by
+;;;; the Free Software Foundation; either version 2, or (at your option)
+;;;; any later version.
+;;;; 
+;;;; This program is distributed in the hope that it will be useful,
+;;;; but WITHOUT ANY WARRANTY; without even the implied warranty of
+;;;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+;;;; GNU General Public License for more details.
+;;;; 
+;;;; You should have received a copy of the GNU General Public License
+;;;; along with this software; see the file COPYING.  If not, write to
+;;;; the Free Software Foundation, Inc., 59 Temple Place, Suite 330,
+;;;; Boston, MA 02111-1307 USA
+
+;; You can easily double-check these numbers by using perl with the Math::BigInt::GMP library.
+(with-test-prefix "modulo-expt"
+   (pass-if (= 1 (modulo-expt 17 23 47)))
+   
+   (pass-if (= 1 (modulo-expt 17 -23 47)))
+   
+   (pass-if (= 17 (modulo-expt 17 -22 47)))
+   
+   (pass-if (= 36 (modulo-expt 17 22 47)))
+   
+   (pass-if (= 183658794479969134816674175082294846241553725240 (modulo-expt 111122223333444455556666 111122223333444455556666 1153478690012629968439432872520758982731022934717)))
+
+   (pass-if-exception
+    "Proper exception with 0 modulus"
+    (cons 'numerical-overflow "")
+    (modulo-expt 17 23 0))
+
+   (pass-if-exception
+    "Proper exception when result not invertible"
+    (cons 'numerical-overflow "")
+    (modulo-expt 10 -1 48))
+)

[-- Attachment #3: Type: text/plain, Size: 139 bytes --]


-- 
Always code as if the guy who ends up maintaining your code will
    be a violent psychopath who knows where you live.  John F. Woods

[-- Attachment #4: Type: text/plain, Size: 142 bytes --]

_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://mail.gnu.org/mailman/listinfo/guile-devel