This is an automated email from the git hooks/post-receive script.

ppm-guest pushed a commit to annotated tag v0.41
in repository libmath-prime-util-perl.

commit 1dc5d3c1b2b0f1eecf1b0a2dfea67e7b2f160898
Author: Dana Jacobsen <d...@acm.org>
Date:   Fri May 2 18:10:20 2014 -0700

    Replace complicated invmod with simpler version that works properly for all 
n
---
 lib/Math/Prime/Util.pm |  6 +++---
 util.c                 | 44 ++++++++++++--------------------------------
 xt/pari-compare.pl     | 11 +++++++++++
 3 files changed, 26 insertions(+), 35 deletions(-)

diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 533e36a..51dce96 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -1165,9 +1165,9 @@ Objects can be passed to functions, and allow early loop 
exits.
   forcomposites { say } 2000,2020;
 
 Given a block and either an end number or a start and end pair, calls the
-block for each composite in the inclusive range.  The composites are the
-numbers greater than 1 which are not prime:
-C<4, 6, 8, 9, 10, 12, 14, 15, ...>
+block for each composite in the inclusive range.  The composites,
+L<OEIS A002808|http://oeis.org/A002808>, are the numbers greater than 1
+which are not prime:  C<4, 6, 8, 9, 10, 12, 14, 15, ...>
 
 
 =head2 fordivisors
diff --git a/util.c b/util.c
index 6286f44..3937473 100644
--- a/util.c
+++ b/util.c
@@ -1481,38 +1481,18 @@ UV znprimroot(UV n) {
   return 0;
 }
 
-/* Calculate 1/a mod p.  From William Hart. */
-UV modinverse(UV a, UV p) {
-  IV u1 = 1, u3 = a;
-  IV v1 = 0, v3 = p;
-  IV t1 = 0, t3 = 0;
-  IV quot;
-  while (v3) {
-    quot = u3 - v3;
-    if (u3 < (v3<<2)) {
-      if (quot < v3) {
-        if (quot < 0) {
-          t1 = u1; u1 = v1; v1 = t1;
-          t3 = u3; u3 = v3; v3 = t3;
-        } else {
-          t1 = u1 - v1; u1 = v1; v1 = t1;
-          t3 = u3 - v3; u3 = v3; v3 = t3;
-        }
-      } else if (quot < (v3<<1)) {
-        t1 = u1 - (v1<<1); u1 = v1; v1 = t1;
-        t3 = u3 - (v3<<1); u3 = v3; v3 = t3;
-      } else {
-        t1 = u1 - v1*3; u1 = v1; v1 = t1;
-        t3 = u3 - v3*3; u3 = v3; v3 = t3;
-      }
-    } else {
-      quot = u3 / v3;
-      t1 = u1 - v1*quot; u1 = v1; v1 = t1;
-      t3 = u3 - v3*quot; u3 = v3; v3 = t3;
-    }
- }
- if (u1 < 0) u1 += p;
- return u1;
+/* Calculate 1/a mod n. */
+UV modinverse(UV a, UV n) {
+  IV t = 0;  UV nt = 1;
+  UV r = n;  UV nr = a;
+  while (nr != 0) {
+    UV quot = r / nr;
+    { UV tmp = nt;  nt = t - quot*nt;  t = tmp; }
+    { UV tmp = nr;  nr = r - quot*nr;  r = tmp; }
+  }
+  if (r > 1) return 1;  /* No inverse */
+  if (t < 0) t += n;
+  return t;
 }
 
 UV divmod(UV a, UV b, UV n) {   /* a / b  mod n */
diff --git a/xt/pari-compare.pl b/xt/pari-compare.pl
index 666cf99..5234c26 100755
--- a/xt/pari-compare.pl
+++ b/xt/pari-compare.pl
@@ -196,6 +196,17 @@ while (1) {
 
   # TODO: carmichael lambda?  Pari doesn't have it.
 
+  # TODO: if we ever export out modular inverse, here's a test:
+  #    use Math::Pari qw/PARI lift/;
+  #    for (1..1000) {
+  #       my $p = random_nbit_prime(64);
+  #       for $i (2..100) {
+  #         my $mpu = znorder($i,$p);
+  #         my $pari = lift(PARI "Mod(1/$i,$p)");
+  #         die "$i $p" unless $pari == $mpu;
+  #       }
+  #     }
+
   { my $m = int(rand($n-1));
     my $mn = PARI "Mod($m,$n)";
     my $order = znorder($m, $n);

-- 
Alioth's /usr/local/bin/git-commit-notice on 
/srv/git.debian.org/git/pkg-perl/packages/libmath-prime-util-perl.git

_______________________________________________
Pkg-perl-cvs-commits mailing list
Pkg-perl-cvs-commits@lists.alioth.debian.org
http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-perl-cvs-commits

Reply via email to