[libmath-prime-util-perl] 25/55: invmod edge cases, add tests
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 174b33419add2a4e9162f4db5969741d2777950f
Author: Dana Jacobsen <dana at acm.org>
Date: Mon May 5 12:36:43 2014 -0700
invmod edge cases, add tests
---
TODO | 2 --
XS.xs | 67 ++++++++++++++++++++++++++-------------------
lib/Math/Prime/Util.pm | 13 ++++++---
lib/Math/Prime/Util/PP.pm | 5 +++-
lib/Math/Prime/Util/PPFE.pm | 6 ++++
t/19-moebius.t | 28 +++++++++++++++++++
t/92-release-pod-coverage.t | 2 +-
7 files changed, 87 insertions(+), 36 deletions(-)
diff --git a/TODO b/TODO
index af89a74..f8132c7 100644
--- a/TODO
+++ b/TODO
@@ -72,5 +72,3 @@
- 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 5dba66b..0e21827 100644
--- a/XS.xs
+++ b/XS.xs
@@ -737,9 +737,8 @@ divisor_sum(IN SV* svn, ...)
void
znorder(IN SV* sva, IN SV* svn)
ALIAS:
- invmod = 1
- jordan_totient = 2
- legendre_phi = 3
+ jordan_totient = 1
+ legendre_phi = 2
PREINIT:
int astatus, nstatus;
PPCODE:
@@ -752,25 +751,22 @@ znorder(IN SV* sva, IN SV* svn)
switch (ix) {
case 0: ret = znorder(a, n);
break;
- case 1: ret = modinverse(a, n);
- break;
- case 2: ret = jordan_totient(a, n);
+ case 1: ret = jordan_totient(a, n);
if (ret == 0 && n > 1)
goto overflow;
break;
- case 3:
+ case 2:
default: ret = legendre_phi(a, n);
break;
}
- if (ret == 0 && (ix == 0 || ix == 1)) XSRETURN_UNDEF; /* not defined */
+ if (ret == 0 && ix == 0) XSRETURN_UNDEF; /* not defined */
XSRETURN_UV(ret);
}
overflow:
switch (ix) {
case 0: _vcallsub_with_gmp("znorder"); break;
- case 1: _vcallsub_with_gmp("invmod"); break;
- case 2: _vcallsub_with_pp("jordan_totient"); break;
- case 3:
+ case 1: _vcallsub_with_pp("jordan_totient"); break;
+ case 2:
default: _vcallsub_with_pp("legendre_phi"); break;
}
return; /* skip implicit PUTBACK */
@@ -796,35 +792,50 @@ void
kronecker(IN SV* sva, IN SV* svb)
ALIAS:
valuation = 1
+ invmod = 2
PREINIT:
int astatus, bstatus, abpositive, abnegative;
PPCODE:
astatus = _validate_int(aTHX_ sva, 2);
bstatus = _validate_int(aTHX_ svb, 2);
- if (ix == 0) {
- /* Are both a and b positive? */
- abpositive = astatus == 1 && bstatus == 1;
- /* Will both fit in IVs? We should use a bitmask return. */
- abnegative = !abpositive
- && (astatus != 0 && SvIOK(sva) && !SvIsUV(sva))
- && (bstatus != 0 && SvIOK(svb) && !SvIsUV(svb));
- if (abpositive || abnegative) {
- UV a = my_svuv(sva);
- UV b = my_svuv(svb);
- int k = (abpositive) ? kronecker_uu(a,b) : kronecker_ss(a,b);
- RETURN_NPARITY(k);
- }
- } else {
- if (astatus != 0 && bstatus != 0) {
+ if (astatus != 0 && bstatus != 0) {
+ if (ix == 0) {
+ /* Are both a and b positive? */
+ abpositive = astatus == 1 && bstatus == 1;
+ /* Will both fit in IVs? We should use a bitmask return. */
+ abnegative = !abpositive
+ && (SvIOK(sva) && !SvIsUV(sva))
+ && (SvIOK(svb) && !SvIsUV(svb));
+ if (abpositive || abnegative) {
+ UV a = my_svuv(sva);
+ UV b = my_svuv(svb);
+ int k = (abpositive) ? kronecker_uu(a,b) : kronecker_ss(a,b);
+ RETURN_NPARITY(k);
+ }
+ } else if (ix == 1) {
UV n = (astatus == -1) ? (UV)(-(my_sviv(sva))) : my_svuv(sva);
UV k = (bstatus == -1) ? (UV)(-(my_sviv(svb))) : my_svuv(svb);
XSRETURN_UV( valuation(n, k) );
+ } else {
+ UV a, n, ret = 0;
+ n = (bstatus != -1) ? my_svuv(svb) : (UV)(-(my_sviv(svb)));
+ if (n > 0) {
+ a = (astatus != -1) ? my_svuv(sva)
+ : n * ((UV)(-my_sviv(sva))/n + 1) + my_sviv(sva);
+ if (a > 0) {
+ if (n == 1) XSRETURN_UV(0);
+ ret = modinverse(a, n);
+ }
+ }
+ if (ret == 0) XSRETURN_UNDEF;
+ XSRETURN_UV(ret);
}
}
switch (ix) {
case 0: _vcallsub_with_gmp("kronecker"); break;
- case 1:
- default: _vcallsub_with_gmp("valuation"); break;
+ case 1: _vcallsub_with_gmp("valuation"); break;
+ case 2:
+ default: _vcallsub_with_gmp("invmod"); break;
}
return; /* skip implicit PUTBACK */
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 9306bb1..33f08fb 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -813,7 +813,7 @@ __END__
=encoding utf8
-=for stopwords forprimes forcomposites fordivisors Möbius Deléglise totient moebius mertens liouville znorder irand primesieve uniqued k-tuples von SoE pari yafu fonction qui compte le nombre nombres voor PhD superset sqrt(N) gcd(A^M k-th (10001st primegen libtommath kronecker znprimroot znlog gcd lcm
+=for stopwords forprimes forcomposites fordivisors Möbius Deléglise totient moebius mertens liouville znorder irand primesieve uniqued k-tuples von SoE pari yafu fonction qui compte le nombre nombres voor PhD superset sqrt(N) gcd(A^M k-th (10001st primegen libtommath kronecker znprimroot znlog gcd lcm invmod
=head1 NAME
@@ -1878,13 +1878,18 @@ 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
+=head2 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.
+If not defined, undef is returned. If defined, then the return value
+multiplied by C<a> equals C<1> modulo C<n>.
+
+This results correspond to the Pari result of C<lift(Mod(1/a,n))>. The
+semantics with respect to negative arguments match Pari. Notably, a
+negative C<n> is negated, which is different from Math::BigInt, but in both
+cases the return value is still congruent to C<1> modulo C<n> as expected.
=head2 valuation
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index 138ca09..023d023 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -1486,8 +1486,11 @@ sub lcm {
}
sub invmod {
my($a,$n) = @_;
+ return if $n == 0 || $a == 0;
+ return 0 if $n == 1;
+ $n = -$n if $n < 0; # Pari semantics
if ($n > (~0>>1)) {
- my $invmod = Math::BigInt->new("$a")->bmodinv($n);
+ 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;
diff --git a/lib/Math/Prime/Util/PPFE.pm b/lib/Math/Prime/Util/PPFE.pm
index e631953..6927f80 100644
--- a/lib/Math/Prime/Util/PPFE.pm
+++ b/lib/Math/Prime/Util/PPFE.pm
@@ -322,6 +322,12 @@ sub lcm {
return Math::Prime::Util::PP::lcm(@_);
}
sub invmod {
+ my ($a, $n) = @_;
+ my ($va, $vn) = ($a, $n);
+ $va = -$va if defined $va && $va < 0;
+ $vn = -$vn if defined $vn && $vn < 0;
+ _validate_positive_integer($va);
+ _validate_positive_integer($vn);
return Math::Prime::Util::PP::invmod(@_);
}
diff --git a/t/19-moebius.t b/t/19-moebius.t
index f9011bb..85a9432 100644
--- a/t/19-moebius.t
+++ b/t/19-moebius.t
@@ -7,6 +7,7 @@ use Math::Prime::Util
qw/moebius mertens euler_phi jordan_totient divisor_sum exp_mangoldt
chebyshev_theta chebyshev_psi carmichael_lambda znorder liouville
znprimroot znlog kronecker legendre_phi gcd lcm is_power valuation
+ invmod
/;
my $extra = defined $ENV{EXTENDED_TESTING} && $ENV{EXTENDED_TESTING};
@@ -366,6 +367,23 @@ if ($use64) {
push @{$powers{9}}, 118587876497;
}
+my @invmods = (
+ [ 0, 0, undef],
+ [ 1, 0, undef],
+ [ 0, 1, undef],
+ [ 1, 1, 0],
+ [ 45, 59, 21],
+ [ 42, 2017, 1969],
+ [ 42, -2017, 1969],
+ [ -42, 2017, 48],
+ [ -42, -2017, 48],
+ [ 14, 28474, undef],
+);
+if ($use64) {
+ push @invmods, [ 13, 9223372036854775808, 5675921253449092805 ];
+ push @invmods, [ 14, 18446744073709551615, 17129119497016012214 ];
+}
+
# These are slow with XS, and *really* slow with PP.
if (!$usexs) {
%big_mertens = map { $_ => $big_mertens{$_} }
@@ -399,6 +417,7 @@ plan tests => 0 + 1
+ scalar(@znlogs)
+ scalar(@legendre_sums)
+ scalar(@valuations)
+ + 3 + scalar(@invmods)
+ scalar(keys %powers)
+ scalar(keys %primroots) + 2
+ scalar(keys %jordan_totients)
@@ -609,6 +628,15 @@ foreach my $r (@valuations) {
my($n, $k, $exp) = @$r;
is( valuation($n, $k), $exp, "valuation($n,$k) = $exp" );
}
+###### invmod
+ok(!eval { invmod(undef,11); }, "invmod(undef,11)");
+ok(!eval { invmod(11,undef); }, "invmod(11,undef)");
+ok(!eval { invmod('nan',11); }, "invmod('nan',11)");
+
+foreach my $r (@invmods) {
+ my($a, $n, $exp) = @$r;
+ is( invmod($a,$n), $exp, "invmod($a,$n) = ".((defined $exp)?$exp:"<undef>") );
+}
sub cmp_closeto {
my $got = shift;
diff --git a/t/92-release-pod-coverage.t b/t/92-release-pod-coverage.t
index eecf464..74eb36b 100644
--- a/t/92-release-pod-coverage.t
+++ b/t/92-release-pod-coverage.t
@@ -65,7 +65,7 @@ sub mpu_public_regex {
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
--
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