[libmath-prime-util-perl] 11/11: Next prime overflow, HOLF factoring, compare with Pari

Partha P. Mukherjee ppm-guest at moszumanska.debian.org
Thu May 21 18:44:18 UTC 2015


This is an automated email from the git hooks/post-receive script.

ppm-guest pushed a commit to annotated tag v0.03
in repository libmath-prime-util-perl.

commit 62a9793e9a6534e7d04f505d41e18ec04d37e3cd
Author: Dana Jacobsen <dana at acm.org>
Date:   Thu Jun 7 05:25:11 2012 -0600

    Next prime overflow, HOLF factoring, compare with Pari
---
 Changes                            |  8 ++++--
 TODO                               | 15 ++++++-----
 XS.xs                              | 26 +++++++++----------
 examples/bench-factor-extra.t      | 35 +++++++++++++++----------
 examples/bench-factor-semiprime.pl | 11 ++++++++
 examples/bench-is-prime.pl         |  4 ++-
 examples/bench-nthprime.pl         |  3 ++-
 examples/test-factoring.pl         |  4 +--
 factor.c                           | 52 ++++++++++++++++++++++++++++++++++----
 factor.h                           |  2 ++
 lib/Math/Prime/Util.pm             | 48 ++++++++++++++++++++++++++++++-----
 t/12-nextprime.t                   |  5 +++-
 util.c                             |  7 +++++
 util.h                             |  5 ++++
 14 files changed, 172 insertions(+), 53 deletions(-)

diff --git a/Changes b/Changes
index 6a51075..cd77577 100644
--- a/Changes
+++ b/Changes
@@ -1,9 +1,13 @@
 Revision history for Perl extension Math::Prime::Util.
 
 0.03  6 June 2012
-    - Speed up factoring a little by moving some work into the XS routine.
+    - Speed up factoring.
     - fixed powmod routine, speedup for smaller numbers
-    - Add Miller-Rabin and deterministic probable prime functions.
+    - Add Miller-Rabin and deterministic probable prime functions.  These
+      are now used for is_prime and factoring, giving a big speedup for
+      numbers > 32-bit.
+    - Add HOLF factoring (just for demo)
+    - Next prime returns 0 on overflow
 
 0.02  5 June 2012
     - Back off new_ok to new/isa_ok to keep Test::More requirements low.
diff --git a/TODO b/TODO
index 9c34838..9ce9eb8 100644
--- a/TODO
+++ b/TODO
@@ -1,16 +1,12 @@
 
 - Examine behavior near 32-bit limit on 32-bit machines.
+  (done for factoring)
 
 - GMP versions of all routines.
 
 - segment sieve should itself use a segment for its primes.
   Today we'd need sqrt(2^64) max = 140MB.  Segmenting would yield under 1MB.
 
-- factoring (Fermat, Pollard Rho, HOLF, etc.)
-- Speed up basic factoring
-- Since is_prime does trials for small numbers, we should skip it in factor
-  for small numbers.  We're just duplicating work.
-
 - Add test to check maxbits in compiled library vs. Perl
 
 - Li(n)
@@ -19,6 +15,11 @@
 
 - input validation (in XS, or do we need to make Perl wrappers for everything?)
 
-- examples and benchmarks
+- Faster SQUFOF
 
-- consider having next prime return 0 when it would overflow
+- random_prime( {bits => 32,
+                 digits => 4,
+                 between => '4 and 100',
+                 between_nth => '26 and 42', } )
+  for super huge bases, we can use approx functions to get center, generate
+  a range, and then do random selection within it.
diff --git a/XS.xs b/XS.xs
index f034c51..f15eb6e 100644
--- a/XS.xs
+++ b/XS.xs
@@ -245,19 +245,14 @@ factor(IN UV n)
           int split_success = 0;
           if (n > UVCONST(60000000) ) {  /* tune this */
             /* For sufficiently large n, try more complex methods. */
-            /* SQUFOF (succeeds ~98% of the time) */
-            split_success = squfof_factor(n, factor_stack+nstack, 64*4096)-1;
+            /* SQUFOF (succeeds 98-99.9%) */
+            split_success = squfof_factor(n, factor_stack+nstack, 256*1024)-1;
             assert( (split_success == 0) || (split_success == 1) );
-            /* a few rounds of Pollard rho (succeeds 99+% of the rest) */
-            if (1 && !split_success) {
+            /* a few rounds of Pollard rho (succeeds most of the rest) */
+            if (!split_success) {
               split_success = prho_factor(n, factor_stack+nstack, 400)-1;
               assert( (split_success == 0) || (split_success == 1) );
             }
-            /* Fermat (Knuth) -- highly debatable with no round limit */
-            if (0 && !split_success) {
-              split_success = fermat_factor(n, factor_stack+nstack,1000)-1;
-              assert( (split_success == 0) || (split_success == 1) );
-            }
           }
           if (split_success) {
             nstack++;
@@ -315,22 +310,27 @@ fermat_factor(IN UV n, IN UV maxrounds = 64*1024*1024)
     SIMPLE_FACTOR(fermat_factor, n, maxrounds);
 
 void
-squfof_factor(IN UV n, IN UV maxrounds = 16*1024*1024)
+holf_factor(IN UV n, IN UV maxrounds = 64*1024*1024)
+  PPCODE:
+    SIMPLE_FACTOR(holf_factor, n, maxrounds);
+
+void
+squfof_factor(IN UV n, IN UV maxrounds = 4*1024*1024)
   PPCODE:
     SIMPLE_FACTOR(squfof_factor, n, maxrounds);
 
 void
-pbrent_factor(IN UV n, IN UV maxrounds = 64*1024*1024)
+pbrent_factor(IN UV n, IN UV maxrounds = 4*1024*1024)
   PPCODE:
     SIMPLE_FACTOR(pbrent_factor, n, maxrounds);
 
 void
-prho_factor(IN UV n, IN UV maxrounds = 64*1024*1024)
+prho_factor(IN UV n, IN UV maxrounds = 4*1024*1024)
   PPCODE:
     SIMPLE_FACTOR(prho_factor, n, maxrounds);
 
 void
-pminus1_factor(IN UV n, IN UV maxrounds = 64*1024*1024)
+pminus1_factor(IN UV n, IN UV maxrounds = 4*1024*1024)
   PPCODE:
     SIMPLE_FACTOR(pminus1_factor, n, maxrounds);
 
diff --git a/examples/bench-factor-extra.t b/examples/bench-factor-extra.t
index 555ac97..18962c4 100755
--- a/examples/bench-factor-extra.t
+++ b/examples/bench-factor-extra.t
@@ -1,23 +1,25 @@
-#!/usr/bin/perl -w
-
+#!/usr/bin/env perl
 use strict;
+use warnings;
 use Math::Prime::Util;
 use Benchmark qw/:all/;
 use List::Util qw/min max/;
 my $count = shift || -2;
+my $is64bit = (~0 > 4294967295);
+my $maxdigits = ($is64bit) ? 20 : 10;  # Noting the range is limited for max.
 
 srand(29);
-my $rounds = 100;
-my $sqrounds = 64*1024;
-test_at_digits($_) for (3..10);
+my $rounds = 400;
+my $sqrounds = 256*1024;
+my $hrounds = 100000;
+test_at_digits($_) for ( 3 .. $maxdigits );
 
 
 sub test_at_digits {
   my $digits = shift;
 
   die "Digits has to be >= 1" unless $digits >= 1;
-  die "Digits has to be <= 10" if (~0 == 4294967295) && ($digits > 10);
-  die "Digits has to be <= 19" if $digits > 19;
+  die "Digits has to be <= $maxdigits" if $digits > $maxdigits;
 
   my @nums = genrand($digits, 1000);
   #my @nums = gensemi($digits, 1000, 23);
@@ -26,15 +28,18 @@ sub test_at_digits {
 
   # Determine success rates
   my %nfactored;
+  my $tfac = 0;
+  my $calc_nfacs = sub { ((scalar grep { $_ > 5 } @_) > 1) ? 1 : 0 };
   for (@nums) {
-    $nfactored{'prho'} += int((scalar grep { $_ > 5 } Math::Prime::Util::prho_factor($_, $rounds))/2);
-    $nfactored{'pbrent'} += int((scalar grep { $_ > 5 } Math::Prime::Util::pbrent_factor($_, $rounds))/2);
-    $nfactored{'pminus1'} += int((scalar grep { $_ > 5 } Math::Prime::Util::pminus1_factor($_, $rounds))/2);
-    $nfactored{'fermat'} += int((scalar grep { $_ > 5 } Math::Prime::Util::fermat_factor($_, $rounds))/2);
-    $nfactored{'squfof'} += int((scalar grep { $_ > 5 } Math::Prime::Util::squfof_factor($_, $sqrounds))/2);
-    #$nfactored{'trial'} += int((scalar grep { $_ > 5 } Math::Prime::Util::trial_factor($_))/2);
+    $tfac += $calc_nfacs->(Math::Prime::Util::factor($_));
+    $nfactored{'prho'} += $calc_nfacs->(Math::Prime::Util::prho_factor($_, $rounds));
+    $nfactored{'pbrent'} += $calc_nfacs->(Math::Prime::Util::pbrent_factor($_, $rounds));
+    $nfactored{'pminus1'} += $calc_nfacs->(Math::Prime::Util::pminus1_factor($_, $rounds));
+    $nfactored{'squfof'} += $calc_nfacs->(Math::Prime::Util::squfof_factor($_, $sqrounds));
+    #$nfactored{'trial'} += $calc_nfacs->(Math::Prime::Util::trial_factor($_));
+    #$nfactored{'fermat'} += $calc_nfacs->(Math::Prime::Util::fermat_factor($_, $rounds));
+    $nfactored{'holf'} += $calc_nfacs->(Math::Prime::Util::holf_factor($_, $hrounds));
   }
-  my $tfac = $nfactored{'fermat'};
 
   print "factoring 1000 random $digits-digit numbers ($min_num - $max_num)\n";
   print "Factorizations: ",
@@ -48,10 +53,12 @@ sub test_at_digits {
     "pbrent"  => sub { Math::Prime::Util::pbrent_factor($_, $rounds) for @nums },
     "pminus1" => sub { Math::Prime::Util::pminus1_factor($_, $rounds) for @nums },
     "fermat"  => sub { Math::Prime::Util::fermat_factor($_, $rounds) for @nums },
+    "holf"    => sub { Math::Prime::Util::holf_factor($_, $hrounds) for @nums },
     "squfof"  => sub { Math::Prime::Util::squfof_factor($_, $sqrounds) for @nums },
     "trial"   => sub { Math::Prime::Util::trial_factor($_) for @nums },
   };
   delete $lref->{'fermat'} if $digits >= 9;
+  #delete $lref->{'holf'} if $digits >= 9;
   cmpthese($count, $lref);
   print "\n";
 }
diff --git a/examples/bench-factor-semiprime.pl b/examples/bench-factor-semiprime.pl
index e78a20c..090adec 100755
--- a/examples/bench-factor-semiprime.pl
+++ b/examples/bench-factor-semiprime.pl
@@ -6,7 +6,9 @@ srand(377);
 
 use Math::Prime::Util qw/factor/;
 use Math::Factor::XS qw/prime_factors/;
+use Math::Pari qw/factorint/;
 use Benchmark qw/:all/;
+use Data::Dumper;
 my $digits = shift || 15;
 my $count = shift || -2;
 
@@ -59,8 +61,17 @@ foreach my $sp (@semiprimes) {
   die "wrong for $sp\n" unless ($#factors == 1) && ($factors[0] * $factors[1]) == $sp;
 }
 print "OK\n";
+print "Verifying Math::Pari $Math::Pari::VERSION ...";
+foreach my $sp (@semiprimes) {
+  my @factors;
+  my ($pn,$pc) = @{factorint($sp)};
+  push @factors, (int($pn->[$_])) x $pc->[$_] for (0 .. $#{$pn});
+  die "wrong for $sp\n" unless ($#factors == 1) && ($factors[0] * $factors[1]) == $sp;
+}
+print "OK\n";
 
 cmpthese($count,{
     'MPU'   => sub { factor($_) for @semiprimes; },
     'MFXS'  => sub { prime_factors($_) for @semiprimes; },
+    'Pari'  => sub { foreach my $n (@semiprimes) { my @factors; my ($pn,$pc) = @{factorint($n)}; push @factors, (int($pn->[$_])) x $pc->[$_] for (0 .. $#{$pn}); } }
 });
diff --git a/examples/bench-is-prime.pl b/examples/bench-is-prime.pl
index 82382fc..cc9694a 100755
--- a/examples/bench-is-prime.pl
+++ b/examples/bench-is-prime.pl
@@ -4,13 +4,14 @@ use strict;
 #use Math::Primality;
 use Math::Prime::XS;
 use Math::Prime::Util;
+use Math::Pari;
 #use Math::Prime::FastSieve;
 use Benchmark qw/:all/;
 use List::Util qw/min max/;
 my $count = shift || -5;
 
 srand(29);
-test_at_digits($_) for (5..10);
+test_at_digits($_) for (5..15);
 
 
 sub test_at_digits {
@@ -35,6 +36,7 @@ sub test_at_digits {
     #'M::P::FS' => sub { $sieve->isprime($_) for @nums },
     'M::P::U' => sub { Math::Prime::Util::is_prime($_) for @nums },
     'MPU prob' => sub { Math::Prime::Util::is_prob_prime($_) for @nums },
+    'Math::Pari' => sub { Math::Pari::isprime($_) for @nums },
   });
   print "\n";
 }
diff --git a/examples/bench-nthprime.pl b/examples/bench-nthprime.pl
index 9ef89ab..55c239a 100644
--- a/examples/bench-nthprime.pl
+++ b/examples/bench-nthprime.pl
@@ -1,8 +1,9 @@
 #!/usr/bin/env perl
 use strict;
 use warnings;
-use Math::Prime::Util qw/nth_prime/;
+use Math::Prime::Util qw/nth_prime prime_precalc/;
 use Devel::TimeThis;
+#prime_precalc(100000000);
 
 foreach my $e (3 .. 9) {
   my $n = 10 ** $e;
diff --git a/examples/test-factoring.pl b/examples/test-factoring.pl
index a6b2664..64fe487 100755
--- a/examples/test-factoring.pl
+++ b/examples/test-factoring.pl
@@ -7,8 +7,8 @@ use Math::Prime::Util qw/factor/;
 use Math::Factor::XS qw/prime_factors/;
 
 my $nlinear = 1000000;
-my $nrandom = 1000000;
-my $randmax = (~0 == 0xFFFFFFFF) ? ~0 : 1_000_000_000_000;
+my $nrandom = shift || 1000000;
+my $randmax = ~0;
 
 print "OK for first 1";
 my $dig = 1;
diff --git a/factor.c b/factor.c
index 6a4503e..3ef5039 100644
--- a/factor.c
+++ b/factor.c
@@ -311,6 +311,43 @@ int fermat_factor(UV n, UV *factors, UV rounds)
   return 1;
 }
 
+/* Hart's One Line Factorization.
+ * Translation from my GMP, missing perfect square calc and premult.
+ */
+int holf_factor(UV n, UV *factors, UV rounds)
+{
+  UV i, s, m, f;
+
+  assert( (n >= 3) && ((n%2) != 0) );
+
+  for (i = 1; i <= rounds; i++) {
+    s = sqrt(n*i);                      /* TODO: overflow here */
+    if ( (s*s) != (n*i) )  s++;
+    m = powsqr(s, n);
+    /* Cheaper would be:
+     *     if (m is probably not a perfect sequare)  continue;
+     *     f = sqrt(m);
+     *     if (f*f == m) { yay }
+     */
+    if ( ((m&2)!= 0) || ((m&7)==5) || ((m&11) == 8) )
+      continue;
+    f = sqrt(m);                        /* expensive */
+    if ( (f*f) == m ) {
+      f = gcd_ui( (s>f) ? s-f : f-s, n);
+      /* This should always succeed, but with overflow concerns.... */
+      if ((f == 1) || (f == n))
+        break;
+      factors[0] = f;
+      factors[1] = n/f;
+      assert( factors[0] * factors[1] == n );
+      return 2;
+    }
+  }
+  factors[0] = n;
+  return 1;
+}
+
+
 /* Pollard / Brent
  *
  * Probabilistic.  If you give this a prime number, it will loop
@@ -332,7 +369,7 @@ int pbrent_factor(UV n, UV *factors, UV rounds)
     default: a = 11; break;
   }
 
-  for (i = 1; i < rounds; i++) {
+  for (i = 1; i <= rounds; i++) {
     Xi = powsqradd(Xi, a, n);
     f = gcd_ui( (Xi>Xm) ? Xi-Xm : Xm-Xi, n);
     if ( (f != 1) && (f != n) ) {
@@ -355,6 +392,7 @@ int pbrent_factor(UV n, UV *factors, UV rounds)
  */
 int prho_factor(UV n, UV *factors, UV rounds)
 {
+  int in_loop = 0;
   UV a, f, i;
   UV U = 7;
   UV V = 7;
@@ -369,12 +407,15 @@ int prho_factor(UV n, UV *factors, UV rounds)
     default: a =  3; break;
   }
 
-  for (i = 1; i < rounds; i++) {
+  for (i = 1; i <= rounds; i++) {
     U = powsqradd(U, a, n);
     V = powsqradd(V, a, n);
     V = powsqradd(V, a, n);
     f = gcd_ui( (U > V) ? U-V : V-U, n);
-    if ( (f != 1) && (f != n) ) {
+    if (f == n) {
+      if (in_loop++)     /* Mark that we've been here */
+        break;           /* Exit now if we're cycling */
+    } else if (f != 1) {
       factors[0] = f;
       factors[1] = n/f;
       assert( factors[0] * factors[1] == n );
@@ -397,7 +438,7 @@ int pminus1_factor(UV n, UV *factors, UV rounds)
 
   assert( (n >= 3) && ((n%2) != 0) );
 
-  for (i = 1; i < rounds; i++) {
+  for (i = 1; i <= rounds; i++) {
     kf = powmod(kf, i, n);
     if (kf == 0) kf = n;
     f = gcd_ui(kf-1, n);
@@ -412,8 +453,9 @@ int pminus1_factor(UV n, UV *factors, UV rounds)
   return 1;
 }
 
+
 /* My modification of Ben Buhrow's modification of Bob Silverman's SQUFOF code.
- * I like Jason P's code a lot -- I should put it in. */
+ */
 static IV qqueue[100+1];
 static IV qpoint;
 static void enqu(IV q, IV *iter) {
diff --git a/factor.h b/factor.h
index 4c1e314..57cedb0 100644
--- a/factor.h
+++ b/factor.h
@@ -9,6 +9,8 @@ extern int trial_factor(UV n, UV *factors, UV maxtrial);
 
 extern int fermat_factor(UV n, UV *factors, UV rounds);
 
+extern int holf_factor(UV n, UV *factors, UV rounds);
+
 extern int squfof_factor(UV n, UV *factors, UV rounds);
 
 extern int pbrent_factor(UV n, UV *factors, UV maxrounds);
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 3f7dd18..6a155c2 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -206,8 +206,14 @@ methods, is_prime, prime_count, nth_prime, approximations and bounds for
 the prime_count and nth prime, next_prime and prev_prime, factoring utilities,
 and more.
 
-The default sieving and factoring are intended to be the fastest on CPAN,
-including Math::Prime::XS, Math::Prime::FastSieve, and Math::Factor::XS.
+All routines currently work in native integers (32-bit or 64-bit).  Bignum
+support may be added later.  If you need bignum support for these types of
+functions inside Perl now, I recommend L<Math::Pari>.
+
+The default sieving and factoring are intended to be (and currently are)
+the fastest on CPAN, including L<Math::Prime::XS>, L<Math::Prime::FastSieve>,
+and L<Math::Factor::XS>.  L<Math::Pari> is slower in some things, faster in
+others.
 
 
 
@@ -246,7 +252,9 @@ using wheel factorization, or a segmented sieve.
 
   $n = next_prime($n);
 
-Returns the next prime greater than the input number.
+Returns the next prime greater than the input number.  0 is returned if the
+next prime is larger than a native integer type (the last primes being
+C<4,294,967,291> in 32-bit Perl and C<18,446,744,073,709,551,557> in 64-bit).
 
 
 =head2 prev_prime
@@ -456,6 +464,18 @@ very fast, but it slows down quite rapidly as the number of digits increases.
 It is very fast for inputs with a factor close to the midpoint
 (e.g. a semiprime p*q where p and q are the same number of digits).
 
+=head2 holf_factor
+
+  my @factors = holf_factor($n);
+
+Produces factors, not necessarily prime, of the positive number input.  An
+optional number of rounds can be given as a second parameter.  It is possible
+the function will be unable to find a factor, in which case a single element,
+the input, is returned.  This uses Hart's One Line Factorization with no
+premultiplier.  It is an interesting alternative to Fermat's algorithm,
+and there are some inputs it can rapidly factor.  In the long run it has the
+same advantages and disadvantages as Fermat's method.
+
 =head2 squfof_factor
 
   my @factors = squfof_factor($n);
@@ -497,6 +517,7 @@ not being big number aware.
 
 Perl versions earlier than 5.8.0 have issues with 64-bit.  The test suite will
 try to determine if your Perl is broken.  This will show up in factoring tests.
+Perl 5.6.2 32-bit works fine, as do later versions with 32-bit and 64-bit.
 
 
 =head1 PERFORMANCE
@@ -508,6 +529,7 @@ Pi(10^10) = 455,052,511.
        5.6  Tomás Oliveira e Silva's segmented sieve v2 (Sep 2010)
        6.6  primegen (optimized Sieve of Atkin)
       11.2  Tomás Oliveira e Silva's segmented sieve v1 (May 2003)
+      17.0  Pari 2.3.5 (primepi)
 
        5.9  My segmented SoE
       15.6  My Sieve of Eratosthenes using a mod-30 wheel
@@ -529,10 +551,22 @@ Perl modules, counting the primes to C<800_000_000> (800 million), in seconds:
    [hours]  Math::Primality             0.04
 
 
-Factoring performance depends on the input, and I'm still tuning it.  Compared
-to Math::Factor::XS, it is faster for small numbers, a little slower in
-the 7-9 digit range, and then gets much faster.  14+ digit semiprimes are
-factored 10 - 100x faster.
+C<is_prime> is slightly faster than M::P::XS for small inputs, but is much
+faster with larger ones  (5x faster for 8-digit, over 50x for 14-digit).
+Math::Pari is relatively slow for small inputs, but becomes faster in the
+13-digit range.
+
+
+Factoring performance depends on the input, and the algorithm choices used
+are still being tuned.  Compared to Math::Factor::XS, it is a tiny bit faster
+for most input under 10M or so, and rapidly gets faster.  For numbers
+larger than 32 bits it's 10-100x faster (depending on the number -- a power
+of two will be identical, while a semiprime with large factors will be on
+the extreme end).  Pari's underlying algorithms and code are very
+sophisticated, and will always be more so than this module, and of course
+supports bignums which is a huge advantage.  Small numbers factor much, much
+faster with Math::Prime::Util.  Pari passes M::P::U in speed somewhere in the
+16 digit range and rapidly increases its lead.
 
 The presentation here:
  L<http://math.boisestate.edu/~liljanab/BOISECRYPTFall09/Jacobsen.pdf>
diff --git a/t/12-nextprime.t b/t/12-nextprime.t
index 27fb925..20c8518 100644
--- a/t/12-nextprime.t
+++ b/t/12-nextprime.t
@@ -8,7 +8,7 @@ use Math::Prime::Util qw/next_prime prev_prime/;
 my $use64 = Math::Prime::Util::_maxbits > 32;
 my $extra = defined $ENV{RELEASE_TESTING} && $ENV{RELEASE_TESTING};
 
-plan tests => 499*2 + 3*2 + 6;
+plan tests => 499*2 + 3*2 + 6 + 2;
 
 my @small_primes = qw/
 2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71
@@ -62,3 +62,6 @@ is( next_prime(19660), 19661, "next prime of 19660 is 19661" );
 is( prev_prime(19662), 19661, "prev prime of 19662 is 19661" );
 is( prev_prime(19660), 19609, "prev prime of 19660 is 19609" );
 is( prev_prime(19610), 19609, "prev prime of 19610 is 19609" );
+
+is( prev_prime(2), 0, "Previous prime of 2 returns 0" );
+is( next_prime(~0-5), 0, "Next prime of ~0-5 returns 0" );
diff --git a/util.c b/util.c
index 3d0c116..cc920a5 100644
--- a/util.c
+++ b/util.c
@@ -194,6 +194,13 @@ UV next_prime(UV n)
     n = sieve_size;
   }
 
+  /* Overflow */
+#if BITS_PER_WORD == 32
+  if (n > UVCONST(4294967291))  return 0;
+#else
+  if (n > UVCONST(18446744073709551557))  return 0;
+#endif
+
   d = n/30;
   m = n - d*30;
   m = nextwheel30[m];  if (m == 1) d++;
diff --git a/util.h b/util.h
index 9737e8c..b74d8af 100644
--- a/util.h
+++ b/util.h
@@ -24,6 +24,11 @@ extern unsigned char* get_prime_segment(void);
 extern void           free_prime_segment(void);
 
 /* Above this value, is_prime will do deterministic Miller-Rabin */
+/* With 64-bit math, we can do much faster mulmods from 2^16-2^32 */
+#if BITS_PER_WORD == 32
 #define MPU_PROB_PRIME_BEST  UVCONST(100000000)
+#else
+#define MPU_PROB_PRIME_BEST  UVCONST(100000)
+#endif
 
 #endif

-- 
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