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 25db44d65d5ff718c7f1ea1abf0bd349f17ee2f9
Author: Dana Jacobsen <d...@acm.org>
Date:   Sun May 4 02:11:01 2014 -0700

    Add invmod
---
 Changes                     |  1 +
 TODO                        |  2 ++
 XS.xs                       | 18 +++++++++++-------
 lib/Math/Prime/Util.pm      | 10 +++++++++-
 lib/Math/Prime/Util/PP.pm   | 19 +++++++++++++++++++
 lib/Math/Prime/Util/PPFE.pm |  3 +++
 util.c                      |  2 +-
 xt/pari-compare.pl          | 37 +++++++++++++++++++++++++------------
 8 files changed, 71 insertions(+), 21 deletions(-)

diff --git a/Changes b/Changes
index 92fa57a..8d49f17 100644
--- a/Changes
+++ b/Changes
@@ -5,6 +5,7 @@ Revision history for Perl module Math::Prime::Util
     [ADDED]
 
     - valuation(n,k)                         how many times does k divide n?
+    - invmod(a,n)                            inverse of a modulo n
 
     [FUNCTIONALITY AND PERFORMANCE]
 
diff --git a/TODO b/TODO
index f8132c7..af89a74 100644
--- a/TODO
+++ b/TODO
@@ -72,3 +72,5 @@
 - We don't use legendre_phi for other functions any more, but it'd be nice
   to speed it up using some ideas from the Ohana 2011 SAGE branch.  For example
   (10**13,10**5) takes 2.5x longer, albeit with 6x less memory.
+
+- tests for invmod, inc negatives, undef, non-numeric, high-bit, and no inv.
diff --git a/XS.xs b/XS.xs
index 9960fde..5dba66b 100644
--- a/XS.xs
+++ b/XS.xs
@@ -737,8 +737,9 @@ divisor_sum(IN SV* svn, ...)
 void
 znorder(IN SV* sva, IN SV* svn)
   ALIAS:
-    jordan_totient = 1
-    legendre_phi = 2
+    invmod = 1
+    jordan_totient = 2
+    legendre_phi = 3
   PREINIT:
     int astatus, nstatus;
   PPCODE:
@@ -750,23 +751,26 @@ znorder(IN SV* sva, IN SV* svn)
       UV ret;
       switch (ix) {
         case 0:  ret = znorder(a, n);
-                 if (ret == 0) XSRETURN_UNDEF;  /* not defined */
                  break;
-        case 1:  ret = jordan_totient(a, n);
+        case 1:  ret = modinverse(a, n);
+                 break;
+        case 2:  ret = jordan_totient(a, n);
                  if (ret == 0 && n > 1)
                    goto overflow;
                  break;
-        case 2:
+        case 3:
         default: ret = legendre_phi(a, n);
                  break;
       }
+      if (ret == 0 && (ix == 0 || ix == 1))  XSRETURN_UNDEF;  /* not defined */
       XSRETURN_UV(ret);
     }
     overflow:
     switch (ix) {
       case 0:  _vcallsub_with_gmp("znorder");  break;
-      case 1:  _vcallsub_with_pp("jordan_totient");  break;
-      case 2:
+      case 1:  _vcallsub_with_gmp("invmod");  break;
+      case 2:  _vcallsub_with_pp("jordan_totient");  break;
+      case 3:
       default: _vcallsub_with_pp("legendre_phi"); break;
     }
     return; /* skip implicit PUTBACK */
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index d58a230..9306bb1 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -40,7 +40,7 @@ our @EXPORT_OK =
       random_maurer_prime random_maurer_prime_with_cert
       random_shawe_taylor_prime random_shawe_taylor_prime_with_cert
       primorial pn_primorial consecutive_integer_lcm
-      gcd lcm factor factor_exp all_factors divisors valuation
+      gcd lcm factor factor_exp all_factors divisors valuation invmod
       moebius mertens euler_phi jordan_totient exp_mangoldt liouville
       partitions
       chebyshev_theta chebyshev_psi
@@ -1878,6 +1878,14 @@ follow the semantics of Mathematica, Pari, and Perl 6, 
re:
   lcm(0, n) = 0              Any zero in list results in zero return
   lcm(n,-m) = lcm(n, m)      We use the absolute values
 
+=head invmod
+
+  say "The inverse of 42 mod 2017 = ", invmod(42,2017);
+
+Given two integers C<a> and C<n>, return the inverse of C<a> modulo C<n>.
+If not defined, undef is returned.  If defined, then multiplying the return
+value by C<a>, modulo C<n>, will equal 1.
+
 =head2 valuation
 
   say "$n is divisible by 2 ", valuation($n,2), " times.";
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index 95cdf9f..08966c1 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -1484,6 +1484,25 @@ sub lcm {
   $lcm = _bigint_to_int($lcm) if $lcm->bacmp(''.~0) <= 0;
   return $lcm;
 }
+sub invmod {
+  my($a,$n) = @_;
+  if ($n > (~0>>1)) {
+    my $invmod = Math::BigInt->new("$a")->bmodinv($n);
+    return if !defined $invmod || $invmod->is_nan;
+    $invmod = _bigint_to_int($invmod) if $invmod->bacmp(''.~0) <= 0;
+    return $invmod;
+  }
+  my($t,$nt,$r,$nr) = (0, 1, $n, $a);
+  use integer;
+  while ($nr != 0) {
+    my $quot = int($r/$nr);
+    ($nt,$t) = ($t-$quot*$nt,$nt);
+    ($nr,$r) = ($r-$quot*$nr,$nr);
+  }
+  return if $r > 1;
+  $t += $n if $t < 0;
+  $t;
+}
 
 # no validation, x is allowed to be negative, y must be >= 0
 sub _gcd_ui {
diff --git a/lib/Math/Prime/Util/PPFE.pm b/lib/Math/Prime/Util/PPFE.pm
index eab46ff..e631953 100644
--- a/lib/Math/Prime/Util/PPFE.pm
+++ b/lib/Math/Prime/Util/PPFE.pm
@@ -321,6 +321,9 @@ sub gcd {
 sub lcm {
   return Math::Prime::Util::PP::lcm(@_);
 }
+sub invmod {
+  return Math::Prime::Util::PP::invmod(@_);
+}
 
 sub legendre_phi {
   my($x, $a) = @_;
diff --git a/util.c b/util.c
index 3937473..69580fe 100644
--- a/util.c
+++ b/util.c
@@ -1490,7 +1490,7 @@ UV modinverse(UV a, UV n) {
     { 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 (r > 1) return 0;  /* No inverse */
   if (t < 0) t += n;
   return t;
 }
diff --git a/xt/pari-compare.pl b/xt/pari-compare.pl
index 5234c26..243319d 100755
--- a/xt/pari-compare.pl
+++ b/xt/pari-compare.pl
@@ -77,9 +77,21 @@ foreach my $n (0 .. $small) {
 
   { my $m = int(rand($n-1));
     my $mn = PARI "Mod($m,$n)";
+    my $invmod = invmod($m, $n);
+    if (defined $invmod) {
+      die "invmod($m, $n)" unless Math::Pari::lift(PARI "Mod(1/$m,$n)") == 
$invmod;
+    } else {
+      eval { PARI "Mod(1/$m,$n)" };
+      die "invmod($m, $n) defined in Pari" unless $@ =~ /impossible inverse/
+        || ($m == 0 && $@ =~ /division by zero/);
+    }
+  }
+
+  { my $m = int(rand($n-1));
+    my $mn = PARI "Mod($m,$n)";
     my $order = znorder($m, $n);
     if (defined $order) {
-      die "znorder($m, $n)" unless Math::Pari::znorder($mn) == znorder($m,$n);
+      die "znorder($m, $n)" unless Math::Pari::znorder($mn) == $order
     } else {
       eval { Math::Pari::znorder($mn); };
       die "znorder($m, $n) defined in Pari" unless $@ =~ /not an element/;
@@ -196,22 +208,23 @@ 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 $invmod = invmod($m, $n);
+    if (defined $invmod) {
+      die "invmod($m, $n)" unless Math::Pari::lift(PARI "Mod(1/$m,$n)") == 
$invmod;
+    } else {
+      eval { PARI "Mod(1/$m,$n)" };
+      die "invmod($m, $n) defined in Pari" unless $@ =~ /impossible inverse/
+        || ($m == 0 && $@ =~ /division by zero/);
+    }
+  }
 
   { my $m = int(rand($n-1));
     my $mn = PARI "Mod($m,$n)";
     my $order = znorder($m, $n);
     if (defined $order) {
-      die "znorder($m, $n)" unless Math::Pari::znorder($mn) == znorder($m,$n);
+      die "znorder($m, $n)" unless Math::Pari::znorder($mn) == $order;
     } else {
       eval { Math::Pari::znorder($mn); };
       die "znorder($m, $n) defined in Pari" unless $@ =~ /not an element/;

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