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 <d...@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

_______________________________________________
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