[libmath-prime-util-perl] 21/55: Replace complicated invmod with simpler version that works properly for all n
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 1dc5d3c1b2b0f1eecf1b0a2dfea67e7b2f160898
Author: Dana Jacobsen <dana at 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
More information about the Pkg-perl-cvs-commits
mailing list