[libmath-prime-util-perl] 23/55: Add invmod
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:53:40 UTC 2015
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 <dana at 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
More information about the Pkg-perl-cvs-commits
mailing list