Forwarded patch for modular exponentiation support (GMP powm)

Forwarded patch for modular exponentiation support (GMP powm)

[Rob Browning]

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


Eric has tried to send this to guile-devel, but for some reason his
messages aren't getting through, so I offered to forward it myself.

Thanks


[-- Attachment #2: Type: message/rfc822, Size: 4399 bytes --]

From: Eric Hanchrow <offby1@blarg.net>
To: guile-devel@gnu.org
Subject: [PATCH]: modular exponentation function
Date: Fri, 23 Jan 2004 15:01:52 -0800
Message-ID: <87vfn2i80v.fsf@offby1.atm01.sea.blarg.net>


I've (almost) never written code for guile before, so I'm probably
doing a bunch of things wrong.  But it works for me, and I think it'd
be useful.

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	10 Jan 2004 21:58:24 -0000
@@ -1519,6 +1519,46 @@
 }
 #undef FUNC_NAME
 
+static SCM
+coerce_to_big(SCM n)
+{
+  if (SCM_BIGP(n))
+    return n;
+  else if (SCM_INUMP(n))
+    {
+      SCM bigger = scm_i_mkbig();
+      mpz_init_set_ui(SCM_I_BIG_MPZ(bigger), SCM_INUM(n));
+      return bigger;
+    }
+  else
+    scm_wrong_type_arg("mexpt", 1, n);
+}
+
+SCM_DEFINE(scm_modular_expt, "mexpt", 3, 0, 0,
+           (SCM n, SCM k, SCM m),
+            "Return @var{n} raised to the non-negative integer exponent\n"
+	    "@var{k}, modulo @var{m}.\n"
+	    "\n"
+	    "@lisp\n"
+	    "(mexpt 2 3 5)\n"
+	    "   @result{} 3\n"
+	    "@end lisp")
+#define FUNC_NAME s_scm_modular_expt
+{
+  SCM result = scm_i_mkbig();
+  n = coerce_to_big(n);
+  k = coerce_to_big(k);
+  m = coerce_to_big(m);
+
+  mpz_powm (SCM_I_BIG_MPZ (result),
+            SCM_I_BIG_MPZ (n),
+            SCM_I_BIG_MPZ (k),
+            SCM_I_BIG_MPZ (m));
+  
+  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	10 Jan 2004 21:58:24 -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
+++ mexp.test	2004-01-10 13:54:58.000000000 -0800
@@ -0,0 +1,27 @@
+;;;; 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
+
+(with-test-prefix "mexpt"
+   (pass-if "right answer with fixnums"
+     (= 1 (mexpt 17 23 47)))
+   
+   (pass-if "right answer with bignums"
+              (= 508153794507026 (mexpt 111122223333444455556666 111122223333444455556666 1234123412341234)))
+)

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

-- 
If you can't change your underwear, can you be sure you have any?

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



-- 
Rob Browning
rlb @defaultvalue.org and @debian.org; previously @cs.utexas.edu
GPG starting 2002-11-03 = 14DD 432F AE39 534D B592  F9A0 25C8 D377 8C7E 73A4

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

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

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

[Kevin Ryde]

Eric Hanchrow <offby1@blarg.net> writes:
>
> +SCM_DEFINE(scm_modular_expt, "mexpt", 3, 0, 0,

Perhaps call it modulo-expt or something, to make it vaguely match the
"modulo" function.

> +  n = coerce_to_big(n);
> +  k = coerce_to_big(k);
> +  m = coerce_to_big(m);

It might be better to make mpz_t temporaries explicitly, since they
can be cleared immediately instead of making the garbage collector
pick them up later.

> +  mpz_powm (SCM_I_BIG_MPZ (result),
> +            SCM_I_BIG_MPZ (n),
> +            SCM_I_BIG_MPZ (k),
> +            SCM_I_BIG_MPZ (m));

Don't forget to check for m==0 and throw an scm error for that.

Also, if m is negative then mpz_powm will object if there's no inverse
of n modulo m, which probably should be done as an scm error rather
than the way gmp provokes a divide by zero.  The case can be detected
by mpz_invert, though it will be wasteful to do the same as what
mpz_powm is about to do.


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

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

[Kevin Ryde]

I wrote:
>
> Also, if m is negative

Oops, I meant the exponent k there.

> The case can be detected
> by mpz_invert, though it will be wasteful to do the same as what
> mpz_powm is about to do.

Actually, I think what you'll want is to call mpz_invert explicitly,
then mpz_powm with the inverse (when it exists) and the absolute value
of the exponent.


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

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

[Kevin Ryde]

Eric Hanchrow <offby1@blarg.net> writes:
>
> See if the attached patch has attained sufficient grooviness.

Please direct all guile devel matters to the guile-devel list, not to
me.  If you're having trouble posting directly you might be able to go
through the gmane.org mail<->news gateway.

> And by
> the way: does the function `coerce_to_big' bother you?  It bothers me,
> because it seems necessary to me, and yet there isn't anything like it
> already.  That makes me wonder why I feel it's necessary, whereas
> everyone else got along fine without it; I figure I'm overlooking
> something.

Usually an inum is handled by separate code and one of the gmp "_ui"
functions.  And an inline mpz_t temporary when there's no _ui
function.

> +  /* if the exponent K is negative, and we simply call mpz_powm, we
> +     might well get a divide-by-zero exception.

... when an inverse 1/n mod m doesn't exist (or is not unique) ...

>       Since those are hard
> +     to handle, we'll do the inversion ourselves -- because that way
> +     we get a simple failure code, which is easy to handle. */
> +  
> +  if (-1 == mpz_sgn(k_tmp))
> +    {
> +      needs_inverting = 1;

You can do the inversion before the powm call, ie. call it with 1/n
and abs(k).  Might be a touch simpler.

> +  mpz_powm (SCM_I_BIG_MPZ (result),
> +            n_tmp,
> +            k_tmp,
> +            m_tmp);

You can use the result destination for one of the temporaries, to save
an mpz_init, if it doesn't get too messy.


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

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

[Eric Hanchrow]

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

This isn't exactly a patch -- rather, it's a new file followed by a
patch.  Anyhow I've incorporated most of Kevin Ryde's suggestions.


[-- Attachment #2: semi-patch --]
[-- Type: application/octet-stream, Size: 4819 bytes --]

;;;; 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))
)
\f
Index: guile-core/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
--- guile-core/libguile/numbers.c	6 Jan 2004 21:55:29 -0000	1.221
+++ guile-core/libguile/numbers.c	11 Feb 2004 06:56:11 -0000
@@ -1519,6 +1519,84 @@
 }
 #undef FUNC_NAME
 
+static void
+coerce_to_big(SCM in, mpz_t out)
+{
+  if (SCM_BIGP(in))
+    mpz_set(out, SCM_I_BIG_MPZ(in));
+  else if (SCM_INUMP(in))
+    {
+      mpz_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 non-negative 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
+{
+  SCM result = scm_i_mkbig();
+  mpz_t n_tmp; 
+  mpz_t k_tmp; 
+  mpz_t m_tmp; 
+  int needs_inverting = 0;
+
+  if (SCM_NFALSEP(scm_zero_p(m))) {
+    scm_num_overflow(FUNC_NAME);
+  }
+  
+  mpz_init(n_tmp);
+  mpz_init(k_tmp);
+  mpz_init(m_tmp);
+
+  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
+     might well get a divide-by-zero exception.  Since those are hard
+     to handle, we'll do the inversion ourselves -- because that way
+     we get a simple failure code, which is easy to handle. */
+  
+  if (-1 == mpz_sgn(k_tmp))
+    {
+      needs_inverting = 1;
+      mpz_neg(k_tmp, k_tmp);
+    }
+
+  mpz_powm (SCM_I_BIG_MPZ (result),
+            n_tmp,
+            k_tmp,
+            m_tmp);
+
+  if (needs_inverting)
+    {
+      mpz_t inverted;
+      mpz_init(inverted);
+      if (!mpz_invert (inverted, SCM_I_BIG_MPZ(result), m_tmp))
+        {
+          mpz_clear(inverted);
+          scm_num_overflow(FUNC_NAME);
+        }
+      mpz_set(SCM_I_BIG_MPZ(result), inverted);
+      mpz_clear(inverted);
+    }
+  mpz_clear(m_tmp);
+  mpz_clear(k_tmp);
+  mpz_clear(n_tmp);
+
+  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: guile-core/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
--- guile-core/libguile/numbers.h	18 Nov 2003 19:59:51 -0000	1.77
+++ guile-core/libguile/numbers.h	11 Feb 2004 06:56:11 -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);

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


-- 
Hamburgers!  The cornerstone of any nutritious breakfast.
        -- Jules {From "Pulp Fiction"}

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

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

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

[Kevin Ryde]

Eric Hanchrow <offby1@blarg.net> writes:
>
> +static void
> +coerce_to_big(SCM in, mpz_t out)
> +{
> +  if (SCM_BIGP(in))
> +    mpz_set(out, SCM_I_BIG_MPZ(in));
> +  else if (SCM_INUMP(in))
> +    {
> +      mpz_set_si(out, SCM_INUM(in));

You could use mpz_init_set and mpz_init_set_si here rather than
separate mpz_init's.

> +      mpz_t inverted;
> +      mpz_init(inverted);
> +      if (!mpz_invert (inverted, SCM_I_BIG_MPZ(result), m_tmp))
> +        {
> +          mpz_clear(inverted);
> +          scm_num_overflow(FUNC_NAME);
> +        }
> +      mpz_set(SCM_I_BIG_MPZ(result), inverted);
> +      mpz_clear(inverted);

I'm pretty sure you don't need the extra mpz_t temporary here.  Just
tell mpz_invert to put the result straight into SCM_I_BIG_MPZ(result).

> +  mpz_clear(m_tmp);
> +  mpz_clear(k_tmp);
> +  mpz_clear(n_tmp);

I think you need this under the inversion failure too, to avoid a
memory leak.


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

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

[Eric Hanchrow]

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

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

[Kevin Ryde]

Eric Hanchrow <offby1@blarg.net> writes:
>
> +coerce_to_big (SCM in, mpz_t out)
> +{
>    ...
> +    scm_wrong_type_arg ("modulo-expt", 1, in);

I think there's another memory leak here.  Eg. if k is an invalid type
then n_tmp previously inited is not cleared.

> +  if (SCM_NFALSEP (scm_zero_p (m)))

Probably don't need to call scm_zero_p for that, just
"m==SCM_MAKINUM(0)" I think, since inexacts etc will be rejected by
the coerce.


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

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

[Eric Hanchrow]

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

Once again.  Kevin, you'll notice I have reverted to my earlier
method of initializing the mpz_t temporaries in the main function,
and setting them in to `coerce_to_big'.  That's because I cannot think
of any other way to have them available to be cleaned up if there's
an error.



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

? libguile/guile.info
? libguile/little-test.scm
? libguile/rydes-criteria.txt
? 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	14 Feb 2004 23:59:05 -0000
@@ -1519,6 +1519,120 @@
 }
 #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_INUMP (in))
+    mpz_set_si (out, SCM_INUM (in));
+  else
+    return 0;
+
+  return 1;
+}
+
+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; 
+    
+  /* 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.  No more than one of these
+     variables may be non-zero.
+  */
+  int report_overflow = 0;
+  int position_of_wrong_type = 0;
+
+  SCM result = SCM_UNDEFINED;
+
+  mpz_init (n_tmp);
+  mpz_init (k_tmp);
+  mpz_init (m_tmp);
+    
+  if (SCM_EQ_P (m, SCM_INUM0))
+    {
+      report_overflow = 1;
+      goto cleanup;
+    }
+  
+  if (!coerce_to_big (n, n_tmp))
+    {
+      position_of_wrong_type = 1;
+      goto cleanup;
+    }
+
+  if (!coerce_to_big (k, k_tmp))
+    {
+      position_of_wrong_type = 2;
+      goto cleanup;
+    }
+
+  if (!coerce_to_big (m, m_tmp))
+    {
+      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);
+ 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,
+                        (1 == position_of_wrong_type
+                         ? n
+                         : (2 == position_of_wrong_type
+                            ? k
+                            : m)));
+  
+  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	14 Feb 2004 23:59:05 -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-14 10:21:55.000000000 -0800
@@ -0,0 +1,58 @@
+;;;; 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))
+
+   (pass-if-exception
+    "Proper exception with wrong type argument"
+    (cons 'wrong-type-arg "")
+    (modulo-expt "Sam" 23 10))
+
+   (pass-if-exception
+    "Proper exception with wrong type argument"
+    (cons 'wrong-type-arg "")
+    (modulo-expt 17 9.9 10))
+   
+    (pass-if-exception
+    "Proper exception with wrong type argument"
+    (cons 'wrong-type-arg "")
+    (modulo-expt 17 23 'Ethel))
+   
+)

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


-- 
But users will not now with glad cries glom on to a language that
gives them no more than what Scheme or Pascal gave them.

        -- Guy Steele, http://www.sun.com/research/jtech/pubs/98-oopsla-growing.ps

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

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

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

[Kevin Ryde]

Eric Hanchrow <offby1@blarg.net> writes:
>
> +  if (!coerce_to_big (k, k_tmp))
> +    {
> +      position_of_wrong_type = 2;
> +      goto cleanup;

I'd be inclined to put the necessary clears inline there,

	if (! coerce_and_init (k, k_tmp))
	  {
	    mpz_clear (n_tmp);
	    SCM_WRONG_TYPE_ARG (SCM_ARG2, k);
	  }

	if (! coerce_and_init (m, m_tmp))
	  {
	    mpz_clear (n_tmp);
	    mpz_clear (k_tmp);
	    SCM_WRONG_TYPE_ARG (SCM_ARG3, m);
	  }

Which seems a bit cleaner to me than setting flags.

Or maybe it's just my prejudices against the sort of idiocies with
flags that I remember advocated in Pascal texts. :-)


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

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

[Marius Vollmer]

Eric Hanchrow <offby1@blarg.net> writes:

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

Sorry for the later answer, guys.  Kevin, what do you say, is this
patch OK?

-- 
GPG: D5D4E405 - 2F9B BCCC 8527 692A 04E3  331E FAF8 226A D5D4 E405


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

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

[Kevin Ryde]

Marius Vollmer <mvo@zagadka.de> writes:
>
> Kevin, what do you say, is this patch OK?

There's a newer version fixing memory leaks on errors.


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

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

[Marius Vollmer]

Kevin Ryde <user42@zip.com.au> writes:

> Marius Vollmer <mvo@zagadka.de> writes:
>>
>> Kevin, what do you say, is this patch OK?
>
> There's a newer version fixing memory leaks on errors.

Yes, is that one OK?  If so, could you apply it?

-- 
GPG: D5D4E405 - 2F9B BCCC 8527 692A 04E3  331E FAF8 226A D5D4 E405


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

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

[Kevin Ryde]

Marius Vollmer <mvo@zagadka.de> writes:
>
> Yes, is that one OK?  If so, could you apply it?

Yep, if Eric posts the latest and greatest and his paperwork is in order.


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

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

[Eric Hanchrow]

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

I'm never sure if my posts are getting to this list, so my apologies
if this appears twice.

Here's the latest & greatest (which I think is unchanged from my last
post).  Alas, I have not filled in the paperwork, and can't remember
what I'm supposed to do.  Is there a standard form, or what?

[-- Attachment #2: modular exponentiation patch --]
[-- Type: application/octet-stream, Size: 6322 bytes --]

Index: libguile/numbers.c
===================================================================
RCS file: /cvsroot/guile/guile/guile-core/libguile/numbers.c,v
retrieving revision 1.224
diff -w -u -r1.224 numbers.c
--- libguile/numbers.c	28 Feb 2004 20:01:13 -0000	1.224
+++ libguile/numbers.c	21 Mar 2004 22:50:50 -0000
@@ -1528,6 +1528,120 @@
 }
 #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_INUMP (in))
+    mpz_set_si (out, SCM_INUM (in));
+  else
+    return 0;
+
+  return 1;
+}
+
+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; 
+    
+  /* 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_EQ_P (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);
+ 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);
+}
+#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	21 Mar 2004 22:50:50 -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-14 10:21:55.000000000 -0800
@@ -0,0 +1,58 @@
+;;;; 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))
+
+   (pass-if-exception
+    "Proper exception with wrong type argument"
+    (cons 'wrong-type-arg "")
+    (modulo-expt "Sam" 23 10))
+
+   (pass-if-exception
+    "Proper exception with wrong type argument"
+    (cons 'wrong-type-arg "")
+    (modulo-expt 17 9.9 10))
+   
+    (pass-if-exception
+    "Proper exception with wrong type argument"
+    (cons 'wrong-type-arg "")
+    (modulo-expt 17 23 'Ethel))
+   
+)

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


-- 
Governing the U.S. is like playing 200 simultaneous chess
matches (while whiny columnists second-guess every move on every
board).
        -- Nicholas Kristof, New York Times columnist

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

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

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

[Kevin Ryde]

Eric Hanchrow <offby1@blarg.net> writes:
>
> Alas, I have not filled in the paperwork, and can't remember
> what I'm supposed to do.  Is there a standard form, or what?

A copyright assignment is preferred, but you can read about the
choices at

	http://www.gnu.org/prep/maintain_5.html#SEC5

Most people go with "past and future" if they're making ongoing
contributions, because it all takes a while going back and forth.  The
form for that is below.  You send it to fsf-records, not this list or
anything.  fsf-records will help you with employer disclaimers etc.




---------------------------------------------------------------------------

request-assign.future:

Please email the following information to fsf-records@gnu.org, and we
will send you the assignment form for your past and future changes.
Please use your full name as the subject line of the message.


[What is the name of the program or package you're contributing to?]


[Did you copy any files or text written by someone else in these changes?
Even if that material is free software, we need to know about it.]


[Do you have an employer who might have a basis to claim to own
your changes?  Do you attend a school which might make such a claim?]


[For the copyright registration, what country are you a citizen of?]


[What year were you born?]


[Please write your email address here.]


[Please write your snail address here.]





[Which files have you changed so far, and which new files have you written
so far?]





---------------------------------------------------------------------------


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

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

[Kevin Ryde]

Eric Hanchrow <offby1@blarg.net> writes:
>
> Here's the latest & greatest

It just occurred to me to wonder if negative divisors should give a
negative remainder, the same way the plain `modulo' does.

mpz_powm ignores the sign of the divisor and always gives back a
positive remainder.


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