[libmath-prime-util-perl] 21/40: Update p+1 factoring, add to tests
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:49:04 UTC 2015
This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.30
in repository libmath-prime-util-perl.
commit d584e6a3f937086e280442bfbb528314d972e638
Author: Dana Jacobsen <dana at acm.org>
Date: Fri Jul 12 12:18:20 2013 -0700
Update p+1 factoring, add to tests
---
TODO | 4 ---
factor.c | 95 ++++++++++++++++++++++++++++----------------------
lib/Math/Prime/Util.pm | 9 +++++
t/50-factoring.t | 3 +-
4 files changed, 64 insertions(+), 47 deletions(-)
diff --git a/TODO b/TODO
index 4b14d75..598ea96 100644
--- a/TODO
+++ b/TODO
@@ -50,10 +50,6 @@
so we can use MULTICALL. Seems to be fixed on 5.18.0. See if this was
broken on 5.16 and previous, or if this was just a 5.17 issue.
-- Tests, documentation, optimizations, and comparisons for p+1 factoring
-
-- Fix p+1 factoring (see latest GMP code).
-
- Switch to new text proofs.
- Sync all the Lucas pseudoprimes with GMP.
diff --git a/factor.c b/factor.c
index bf46e85..b83da1a 100644
--- a/factor.c
+++ b/factor.c
@@ -993,51 +993,62 @@ int pminus1_factor(UV n, UV *factors, UV B1, UV B2)
}
/* Simple Williams p+1 */
-int pplus1_factor(UV n, UV *factors, UV B)
+static void pp1_pow(UV *cX, unsigned long exp, UV n)
{
- UV i, f, b;
- UV P1, P2, P3, U1, U2, U3, V1, V2, V3;
-
- /* Calculate 3 sequences at once */
- P1 = 5 % n;
- P2 = 9 % n;
- P3 = 13 % n;
- U1 = P1; U2 = P2; U3 = P3;
- for (i = 1; i < B; i++) {
- { UV v = i; b = 1; while (v >>= 1) b++; }
- V1 = mulsubmod(P1, P1, 2, n);
- V2 = mulsubmod(P2, P2, 2, n);
- V3 = mulsubmod(P3, P3, 2, n);
- while (b > 1) {
- b--;
- if ( (i >> (b-1)) & UVCONST(1) ) {
- U1 = mulsubmod(U1, V1, P1, n);
- V1 = mulsubmod(V1, V1, 2, n);
- U2 = mulsubmod(U2, V2, P2, n);
- V2 = mulsubmod(V2, V2, 2, n);
- U3 = mulsubmod(U3, V3, P3, n);
- V3 = mulsubmod(V3, V3, 2, n);
- } else {
- V1 = mulsubmod(U1, V1, P1, n);
- U1 = mulsubmod(U1, U1, 2, n);
- V2 = mulsubmod(U2, V2, P2, n);
- U2 = mulsubmod(U2, U2, 2, n);
- V3 = mulsubmod(U3, V3, P3, n);
- U3 = mulsubmod(U3, U3, 2, n);
- }
+ UV X0 = *cX;
+ UV X = *cX;
+ UV Y = mulsubmod(X, X, 2, n);
+ unsigned long bit;
+ {
+ unsigned long v = exp;
+ unsigned long b = 1;
+ while (v >>= 1) b++;
+ bit = 1UL << (b-2);
+ }
+ while (bit) {
+ if ( exp & bit ) {
+ X = mulsubmod(X, Y, X0, n);
+ Y = mulsubmod(Y, Y, 2, n);
+ } else {
+ Y = mulsubmod(X, Y, X0, n);
+ X = mulsubmod(X, X, 2, n);
}
- P1 = U1; P2 = U2; P3 = U3;
- f = gcd_ui(n, submod(U1, 2, n));
- if (f == 1 || f == n)
- f = gcd_ui(n, submod(U2, 2, n));
- if (f == 1 || f == n)
- f = gcd_ui(n, submod(U3, 2, n));
- if (f > 1 && f != n) {
- factors[0] = f;
- factors[1] = n/f;
- MPUassert( factors[0] * factors[1] == n , "incorrect factoring");
- return 2;
+ bit >>= 1;
+ }
+ *cX = X;
+}
+int pplus1_factor(UV n, UV *factors, UV B1)
+{
+ UV X1, X2, f;
+ UV sqrtB1 = isqrt(B1);
+ MPUassert( (n >= 3) && ((n%2) != 0) , "bad n in pminus1_factor");
+
+ X1 = 7 % n;
+ X2 = 11 % n;
+ START_DO_FOR_EACH_PRIME(2, B1) {
+ UV k = p;
+ if (p < sqrtB1) {
+ UV kmin = B1/p;
+ while (k <= kmin)
+ k *= p;
+ }
+ pp1_pow(&X1, k, n);
+ if (X1 != 2) {
+ f = gcd_ui( submod(X1, 2, n) , n);
+ if (f != 1 && f != n) break;
+ }
+ pp1_pow(&X2, k, n);
+ if (X2 != 2) {
+ f = gcd_ui( submod(X2, 2, n) , n);
+ if (f != 1 && f != n) break;
}
+ } END_DO_FOR_EACH_PRIME
+
+ if ( (f != 1) && (f != n) ) {
+ factors[0] = f;
+ factors[1] = n/f;
+ MPUassert( factors[0] * factors[1] == n , "incorrect factoring");
+ return 2;
}
factors[0] = n;
return 1;
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 3191dc0..c8304be 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -114,6 +114,7 @@ BEGIN {
*pbrent_factor = \&Math::Prime::Util::PP::pbrent_factor;
*prho_factor = \&Math::Prime::Util::PP::prho_factor;
*pminus1_factor = \&Math::Prime::Util::PP::pminus1_factor;
+ *pplus1_factor = \&Math::Prime::Util::PP::pminus1_factor; # TODO: implement PP p+1.
};
$_Config{'nobigint'} = 0;
@@ -3813,6 +3814,14 @@ is Pollard's C<p-1> method, using two stages with default smoothness
settings of 1_000_000 for B1, and C<10 * B1> for B2. This method can rapidly
find a factor C<p> of C<n> where C<p-1> is smooth (it has no large factors).
+=head2 pplus1_factor
+
+ my @factors = pplus1_factor($n);
+ my @factors = pplus1_factor($n, 1_000); # set B1 smoothness
+
+Produces factors, not necessarily prime, of the positive number input. This
+is Williams' C<p+1> method, using one stage and two predefined initial points.
+
=head1 MATHEMATICAL FUNCTIONS
diff --git a/t/50-factoring.t b/t/50-factoring.t
index 3073555..8365d7a 100644
--- a/t/50-factoring.t
+++ b/t/50-factoring.t
@@ -74,7 +74,7 @@ my %all_factors = (
0 => [],
);
-plan tests => (2 * scalar @testn) + scalar(keys %all_factors) + 10*8 + 1;
+plan tests => (2 * scalar @testn) + scalar(keys %all_factors) + 10*9 + 1;
foreach my $n (@testn) {
my @f = factor($n);
@@ -106,6 +106,7 @@ extra_factor_test("rsqufof_factor", sub {Math::Prime::Util::rsqufof_factor(shift
extra_factor_test("pbrent_factor", sub {Math::Prime::Util::pbrent_factor(shift)});
extra_factor_test("prho_factor", sub {Math::Prime::Util::prho_factor(shift)});
extra_factor_test("pminus1_factor",sub {Math::Prime::Util::pminus1_factor(shift)});
+extra_factor_test("pplus1_factor", sub {Math::Prime::Util::pplus1_factor(shift)});
# To hit some extra coverage
is_deeply( [Math::Prime::Util::trial_factor(5514109)], [2203,2503], "trial factor 2203*2503" );
--
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