[libmath-prime-util-perl] 07/11: Redo factor loop
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:44:17 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 7d6012cdbe91163a76e45c2b7074ca51802ee14c
Author: Dana Jacobsen <dana at acm.org>
Date: Wed Jun 6 14:21:09 2012 -0600
Redo factor loop
---
TODO | 4 +
XS.xs | 178 +++++++++++++++++++--------------------------
examples/test-factoring.pl | 9 ++-
factor.c | 80 ++++----------------
factor.h | 4 +-
lib/Math/Prime/Util.pm | 83 +++++++++++++--------
util.c | 32 +++++---
7 files changed, 177 insertions(+), 213 deletions(-)
diff --git a/TODO b/TODO
index e155880..9c34838 100644
--- a/TODO
+++ b/TODO
@@ -8,6 +8,8 @@
- 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
@@ -18,3 +20,5 @@
- input validation (in XS, or do we need to make Perl wrappers for everything?)
- examples and benchmarks
+
+- consider having next prime return 0 when it would overflow
diff --git a/XS.xs b/XS.xs
index 7027934..0115743 100644
--- a/XS.xs
+++ b/XS.xs
@@ -227,12 +227,12 @@ erat_simple_primes(IN UV low, IN UV high)
void
factor(IN UV n)
- PREINIT:
PPCODE:
- /* Exit if n is 0 or 1 */
- if (n < 2) {
- XPUSHs(sv_2mortal(newSVuv( n )));
+ if (n < 4) {
+ XPUSHs(sv_2mortal(newSVuv( n ))); /* If n is 0-3, we're done. */
} else {
+ UV factor_stack[MPU_MAX_FACTORS+1];
+ int nstack = 0;
/* Quick trial divisions. We could do tricky gcd magic here. */
while ( (n% 2) == 0 ) { n /= 2; XPUSHs(sv_2mortal(newSVuv( 2 ))); }
while ( (n% 3) == 0 ) { n /= 3; XPUSHs(sv_2mortal(newSVuv( 3 ))); }
@@ -241,125 +241,99 @@ factor(IN UV n)
while ( (n%11) == 0 ) { n /= 11; XPUSHs(sv_2mortal(newSVuv( 11 ))); }
while ( (n%13) == 0 ) { n /= 13; XPUSHs(sv_2mortal(newSVuv( 13 ))); }
while ( (n%17) == 0 ) { n /= 17; XPUSHs(sv_2mortal(newSVuv( 17 ))); }
- if (n < (19*19)) {
- if (n != 1) XPUSHs(sv_2mortal(newSVuv( n )));
- } else {
- UV startfactor = 19;
- UV factor_list[MPU_MAX_FACTORS+1];
- int n_list_factors = 1;
- int i;
- if (n > (UVCONST(0xFFFFFFFF)) ) { /* tune this */
- /* Use SQUFOF to crack some factors, recurse on large factors */
- n_list_factors = squfof_factor(n, factor_list, 800000, (373*373)+1);
- } else {
- factor_list[0] = n;
- }
- for (i = 0; i < n_list_factors; i++) {
- /* trial factorization for each item in the list */
- UV f, m, limit;
- n = factor_list[i];
- /* We pulled out everything through this point, so must be prime */
- if (n < (19*19)) {
- XPUSHs(sv_2mortal(newSVuv( n )));
- continue;
- /* If sufficiently large, run a prob prime test on it. This saves
- * us a huge amount of work on big primes. We could also look into
- * some possible shortcuts. What this really needs is to move
- * the prime test up above SQUFOF. */
- } else if ( (n > 40000000) && is_prob_prime(n) ) {
- XPUSHs(sv_2mortal(newSVuv( n )));
- continue;
+ do { /* loop over each remaining factor */
+ while ( (n >= (19*19)) && (!is_prime(n)) ) {
+ /* n is composite, so factor it. */
+ int split_success = 0;
+ if (n > UVCONST(10000000) ) { /* tune this */
+ /* For sufficiently large numbers, try more complex methods. */
+ /* SQUFOF (succeeds ~98% of the time) */
+ split_success = squfof_factor(n, factor_stack+nstack, 64*4096)-1;
+ assert( (split_success == 0) || (split_success == 1) );
+ /* a few rounds of Pollard rho (succeeds 99+% 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) );
+ }
}
- f = startfactor;
- m = startfactor % 30;
- limit = sqrt((double) n);
- while (f <= limit) {
- if ( (n%f) == 0 ) {
- do {
- n /= f; XPUSHs(sv_2mortal(newSVuv( f )));
- } while ( (n%f) == 0 );
- limit = sqrt((double) n);
+ if (split_success) {
+ nstack++;
+ n = factor_stack[nstack];
+ } else {
+ /* trial divisions */
+ UV f = 19;
+ UV m = 19;
+ UV limit = sqrt((double) n);
+ while (f <= limit) {
+ if ( (n%f) == 0 ) {
+ do {
+ n /= f; XPUSHs(sv_2mortal(newSVuv( f )));
+ } while ( (n%f) == 0 );
+ limit = sqrt((double) n);
+ }
+ f += wheeladvance30[m];
+ m = nextwheel30[m];
}
- f += wheeladvance30[m];
- m = nextwheel30[m];
}
- if (n != 1)
- XPUSHs(sv_2mortal(newSVuv( n )));
}
- }
+ if (n != 1) XPUSHs(sv_2mortal(newSVuv( n )));
+ if (nstack > 0) n = factor_stack[nstack-1];
+ } while (nstack-- > 0);
+ }
+
+#define SIMPLE_FACTOR(func, n, rounds) \
+ if (n <= 1) { \
+ XPUSHs(sv_2mortal(newSVuv( n ))); \
+ } else { \
+ while ( (n% 2) == 0 ) { n /= 2; XPUSHs(sv_2mortal(newSVuv( 2 ))); } \
+ while ( (n% 3) == 0 ) { n /= 3; XPUSHs(sv_2mortal(newSVuv( 3 ))); } \
+ while ( (n% 5) == 0 ) { n /= 5; XPUSHs(sv_2mortal(newSVuv( 5 ))); } \
+ if (n == 1) { /* done */ } \
+ else if (is_prime(n)) { XPUSHs(sv_2mortal(newSVuv( n ))); } \
+ else { \
+ UV factors[MPU_MAX_FACTORS+1]; \
+ int i, nfactors; \
+ nfactors = func(n, factors, rounds); \
+ for (i = 0; i < nfactors; i++) { \
+ XPUSHs(sv_2mortal(newSVuv( factors[i] ))); \
+ } \
+ } \
}
void
-fermat_factor(IN UV n)
- PREINIT:
- UV factors[MPU_MAX_FACTORS+1];
- int nfactors;
- int i;
+trial_factor(IN UV n)
PPCODE:
- nfactors = fermat_factor(n, factors);
- if (nfactors < 1)
- croak("No factors from fermat_factor");
- for (i = 0; i < nfactors; i++) {
- XPUSHs(sv_2mortal(newSVuv( factors[i] )));
- }
+ SIMPLE_FACTOR(trial_factor, n, UV_MAX);
void
-squfof_factor(IN UV n)
- PREINIT:
- UV factors[MPU_MAX_FACTORS+1];
- int nfactors;
- int i;
+fermat_factor(IN UV n, IN UV maxrounds = 64*1024*1024)
PPCODE:
- /* have squfof recurse on each factor found */
- nfactors = squfof_factor(n, factors, 800000, 2);
- if (nfactors < 1)
- croak("No factors from squfof_factor");
- for (i = 0; i < nfactors; i++) {
- XPUSHs(sv_2mortal(newSVuv( factors[i] )));
- }
+ SIMPLE_FACTOR(fermat_factor, n, maxrounds);
+void
+squfof_factor(IN UV n, IN UV maxrounds = 16*1024*1024)
+ PPCODE:
+ SIMPLE_FACTOR(squfof_factor, n, maxrounds);
void
-pbrent_factor(IN UV n)
- PREINIT:
- UV factors[MPU_MAX_FACTORS+1];
- int nfactors;
- int i;
+pbrent_factor(IN UV n, IN UV maxrounds = 64*1024*1024)
PPCODE:
- nfactors = pbrent_factor(n, factors, 50000000);
- if (nfactors < 1)
- croak("No factors from pbrent_factor");
- for (i = 0; i < nfactors; i++) {
- XPUSHs(sv_2mortal(newSVuv( factors[i] )));
- }
+ SIMPLE_FACTOR(pbrent_factor, n, maxrounds);
void
-prho_factor(IN UV n)
- PREINIT:
- UV factors[MPU_MAX_FACTORS+1];
- int nfactors;
- int i;
+prho_factor(IN UV n, IN UV maxrounds = 64*1024*1024)
PPCODE:
- nfactors = prho_factor(n, factors, 50000000);
- if (nfactors < 1)
- croak("No factors from prho_factor");
- for (i = 0; i < nfactors; i++) {
- XPUSHs(sv_2mortal(newSVuv( factors[i] )));
- }
+ SIMPLE_FACTOR(prho_factor, n, maxrounds);
void
-pminus1_factor(IN UV n)
- PREINIT:
- UV factors[MPU_MAX_FACTORS+1];
- int nfactors;
- int i;
+pminus1_factor(IN UV n, IN UV maxrounds = 64*1024*1024)
PPCODE:
- nfactors = pminus1_factor(n, factors, 50000000);
- if (nfactors < 1)
- croak("No factors from pminus1_factor");
- for (i = 0; i < nfactors; i++) {
- XPUSHs(sv_2mortal(newSVuv( factors[i] )));
- }
+ SIMPLE_FACTOR(pminus1_factor, n, maxrounds);
int
miller_rabin(IN UV n, ...)
diff --git a/examples/test-factoring.pl b/examples/test-factoring.pl
index f8c5d3d..a6b2664 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 = ~0;
-my $randmax = 1_000_000_000_000;
+my $nrandom = 1000000;
+my $randmax = (~0 == 0xFFFFFFFF) ? ~0 : 1_000_000_000_000;
print "OK for first 1";
my $dig = 1;
@@ -25,13 +25,14 @@ foreach my $n (2 .. $nlinear) {
}
}
print " numbers\n";
-print "Testing random numbers from $nlinear to ", $nlinear+$randmax, "\n";
+print "Testing random numbers from $nlinear to ", $randmax, "\n";
while ($nrandom-- > 0) {
- my $n = $nlinear + int(rand($randmax));
+ my $n = $nlinear + 1 + int(rand($randmax - $nlinear));
my @mfxs = sort { $a<=>$b } prime_factors($n);
my @mpu = sort { $a<=>$b } factor($n);
die "failure for $n" unless scalar @mfxs == scalar @mpu;
for (0 .. $#mfxs) { die "failure for $n" unless $mfxs[$_] == $mpu[$_]; }
print "." if ($nrandom % 1024) == 0;
}
+print "\n";
diff --git a/factor.c b/factor.c
index 546123e..afbc7db 100644
--- a/factor.c
+++ b/factor.c
@@ -199,9 +199,14 @@ int is_prob_prime(UV n)
#if 1
/* Better basis from: http://miller-rabin.appspot.com/ */
if (n < UVCONST(9080191)) {
- bases[0] = 31; bases[1] = 73; nbases = 2;
+ bases[0] = 31;
+ bases[1] = 73;
+ nbases = 2;
} else if (n < UVCONST(4759123141)) {
- bases[0] = 2; bases[1] = 7; bases[2] = 61; nbases = 3;
+ bases[0] = 2;
+ bases[1] = 7;
+ bases[2] = 61;
+ nbases = 3;
} else if (n < UVCONST(105936894253)) {
bases[0] = 2;
bases[1] = 1005905886;
@@ -263,23 +268,12 @@ int is_prob_prime(UV n)
* Very fast for small numbers, grows rapidly.
* SQUFOF is better for numbers nearing the 64-bit limit.
*/
-int fermat_factor(UV n, UV *factors)
+int fermat_factor(UV n, UV *factors, UV rounds)
{
int nfactors = 0;
IV sqn, x, y, r;
- if (n < 2) {
- factors[0] = n;
- return 1;
- }
-
- while ((n & 1) == 0) {
- factors[nfactors++] = 2;
- n /= 2;
- }
-
- if (n == 1)
- return nfactors;
+ assert( (n >= 3) && ((n%2) != 0) );
sqn = sqrt((double) n);
x = 2 * sqn + 1;
@@ -313,18 +307,7 @@ int pbrent_factor(UV n, UV *factors, UV rounds)
int nfactors = 0;
UV a, f, Xi, Xm, i;
- if (n < 2) {
- factors[0] = n;
- return 1;
- }
-
- while ((n & 1) == 0) {
- factors[nfactors++] = 2;
- n /= 2;
- }
-
- if (n == 1)
- return nfactors;
+ assert( (n >= 3) && ((n%2) != 0) );
Xi = 2;
Xm = 2;
@@ -361,18 +344,7 @@ int prho_factor(UV n, UV *factors, UV rounds)
int nfactors = 0;
UV a, f, t, U, V, i;
- if (n < 2) {
- factors[0] = n;
- return 1;
- }
-
- while ((n & 1) == 0) {
- factors[nfactors++] = 2;
- n /= 2;
- }
-
- if (n == 1)
- return nfactors;
+ assert( (n >= 3) && ((n%2) != 0) );
switch (n%4) {
case 0: a = 5; break;
@@ -411,18 +383,7 @@ int pminus1_factor(UV n, UV *factors, UV rounds)
int nfactors = 0;
UV f, b, i;
- if (n < 2) {
- factors[0] = n;
- return 1;
- }
-
- while ((n & 1) == 0) {
- factors[nfactors++] = 2;
- n /= 2;
- }
-
- if (n == 1)
- return nfactors;
+ assert( (n >= 3) && ((n%2) != 0) );
b = 13;
for (i = 1; i < rounds; i++) {
@@ -450,7 +411,7 @@ static void enqu(IV q, IV *iter) {
if (++qpoint >= 100) *iter = -1;
}
-int squfof_factor(UV n, UV *factors, UV rounds, UV refactor_above)
+int squfof_factor(UV n, UV *factors, UV rounds)
{
int nfactors = 0;
UV rounds2 = rounds/16;
@@ -459,18 +420,7 @@ int squfof_factor(UV n, UV *factors, UV rounds, UV refactor_above)
IV jter, iter;
int reloop;
- if ( (n < 2) ) {
- factors[0] = n;
- return 1;
- }
-
- while ((n & 1) == 0) {
- factors[nfactors++] = 2;
- n /= 2;
- }
-
- if (n == 1)
- return nfactors;
+ assert( (n >= 3) && ((n%2) != 0) );
/* TODO: What value of n leads to overflow? */
@@ -583,7 +533,7 @@ int squfof_factor(UV n, UV *factors, UV rounds, UV refactor_above)
/* printf(" squfof found %lu = %lu * %lu in %ld/%ld rounds\n", n, p, q, jter, iter); */
-#if 0
+#if 1
factors[nfactors++] = p;
factors[nfactors++] = q;
#else
diff --git a/factor.h b/factor.h
index 6c995e7..4c1e314 100644
--- a/factor.h
+++ b/factor.h
@@ -7,9 +7,9 @@
extern int trial_factor(UV n, UV *factors, UV maxtrial);
-extern int fermat_factor(UV n, UV *factors);
+extern int fermat_factor(UV n, UV *factors, UV rounds);
-extern int squfof_factor(UV n, UV *factors, UV rounds, UV refactor_above);
+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 998665a..3f7dd18 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -362,10 +362,13 @@ between C<2> and C<n - 2>. Returns 2 is C<n> is definitely prime, 1 if C<n>
is probably prime, and 0 if C<n> is definitely composite. Since this is
just the Miller-Rabin test, a value of 2 is only returned for inputs of
2 and 3, which are shortcut. If 0 is returned, then the number really is a
-composite. If 1 is returned, we aren't sure.
+composite. If 1 is returned, there is a good chance the number is prime
+(depending on the input and the bases), but we cannot be sure.
This is usually used in combination with other tests to make either stronger
-tests or deterministic results for numbers less than some verified limit.
+tests (e.g. the strong BPSW test) or deterministic results for numbers less
+than some verified limit (such as the C<is_prob_prime> function in this module).
+
=head2 is_prob_prime
@@ -376,7 +379,7 @@ Takes a positive number as input and returns back either 0 (composite),
2 (definitely prime), or 1 (probably prime).
This is done with a tuned set of Miller-Rabin tests such that the result
-should be deterministic for 64-bit input. Either 2, 3, 4, 5, or 7 Miller-Rabin
+will be deterministic for 64-bit input. Either 2, 3, 4, 5, or 7 Miller-Rabin
tests are performed (no more than 3 for 32-bit input), and the result will
then always be 0 (composite) or 2 (prime). A later implementation may switch
to a BPSW test, depending on speed.
@@ -423,17 +426,25 @@ in place because you still have an object.
my @factors = factor(3_369_738_766_071_892_021);
-Produces the prime factors of a positive number input. They will typically
-but necessarily be in numerical order. The special cases of C<n = 0> and
-C<n = 1> will return C<n>, which guarantees multiplying the factors together
-will always result in the input value, though those are the only cases where
+Produces the prime factors of a positive number input. They may not be in
+numerical order. The special cases of C<n = 0> and C<n = 1> will
+return C<n>, which guarantees multiplying the factors together will
+always result in the input value, though those are the only cases where
the returned factors are not prime.
-The current algorithm is to use trial division for 32-bit numbers, while for
-larger numbers a small set of trial divisions is performed, followed by a
-single run of SQUFOF, then trial division of the results. This results in
-faster factorization of most large numbers. More sophisticated methods could
-be used.
+The current algorithm is to use trial division for small numbers, while large
+numbers go through a sequence of small trials, SQUFOF, and Pollard's Rho.
+This process is repeated for each non-prime factor.
+
+
+=head2 trial_factor
+
+ my @factors = trial_factor($n);
+
+Produces the prime factors of a positive number input. The factors will be
+in numerical order. The special cases of C<n = 0> and C<n = 1> will return
+C<n>, while with all other inputs the factors are guaranteed to be prime.
+For large inputs this will be very slow.
=head2 fermat_factor
@@ -442,16 +453,17 @@ be used.
Produces factors, not necessarily prime, of the positive number input. The
particular algorithm is Knuth's algorithm C. For small inputs this will be
very fast, but it slows down quite rapidly as the number of digits increases.
-If there is a factor close to the midpoint (e.g. a semiprime p*q where p and
-q are the same number of digits), then this will be very fast.
+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 squfof_factor
my @factors = squfof_factor($n);
-Produces factors, not necessarily prime, of the positive number input. It
-is possible the function will be unable to find a factor, in which case a
-single element, the input, is returned.
+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 function typically runs very fast.
=head2 prho_factor
@@ -459,20 +471,20 @@ single element, the input, is returned.
=head2 pminus1_factor
-Attempts to find a factor using one of the probabilistic algorigthms of
+ my @factors = prho_factor($n);
+
+ # Use a very small number of rounds
+ my @factors = prho_factor($n, 1000);
+
+Produces factors, not necessarily prime, of the positive number input. An
+optional number of rounds can be given as a second parameter. These attempt
+to find a single factor using one of the probabilistic algorigthms of
Pollard Rho, Brent's modification of Pollard Rho, or Pollard's C<p - 1>.
These are more specialized algorithms usually used for pre-factoring very
-large inputs, or checking very large inputs for naive mistakes. If given
-a prime input, or if just unlucky, these will take a long time to return
-back the single input value. There are cases where these can result in
-finding a factor or very large inputs in remarkably short time, similar to
-how Fermat's method works very well for factors near C<sqrt(n)>. They are
-also amenable to massively parallel searching.
-
-For 64-bit input, these are unlikely to be of much use. An optimized SQUFOF
-implementation takes under 20 milliseconds to find a factor for any 62-bit
-input on modern desktop computers. Lightweight quadratic sieves are
-typically much faster for inputs in the 19+ digit range.
+large inputs, or checking very large inputs for naive mistakes. If the
+input is prime or they run out of rounds, they will return the single
+input value. On some inputs they will take a very long time, while on
+others they succeed in a remarkably short time.
=head1 LIMITATIONS
@@ -497,7 +509,7 @@ Pi(10^10) = 455,052,511.
6.6 primegen (optimized Sieve of Atkin)
11.2 Tomás Oliveira e Silva's segmented sieve v1 (May 2003)
- 5.9 My SoE called in segments
+ 5.9 My segmented SoE
15.6 My Sieve of Eratosthenes using a mod-30 wheel
17.2 A slightly modified verion of Terje Mathisen's mod-30 sieve
35.5 Basic Sieve of Eratosthenes on odd numbers
@@ -522,6 +534,17 @@ 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.
+The presentation here:
+ L<http://math.boisestate.edu/~liljanab/BOISECRYPTFall09/Jacobsen.pdf>
+has a lot of data on 64-bit and GMP factoring performance I collected in 2009.
+Assuming you do not know anything about the inputs, trial division and
+optimized Fermat work very well for small numbers (<= 10 digits), while
+native SQUFOF is typically the method of choice for 11-18 digits (I've
+seen claims that a lightweight QS can be faster for 15+ digits). Some form
+of Quadratic Sieve is usually used for inputs in the 19-100 digit range, and
+beyond that is the Generalized Number Field Sieve. For serious factoring,
+I recommend looking info C<yafu>, C<msieve>, C<Pari>, and C<GGNFS>.
+
=head1 AUTHORS
diff --git a/util.c b/util.c
index 95c0a00..1295d94 100644
--- a/util.c
+++ b/util.c
@@ -46,9 +46,30 @@ static UV count_zero_bits(const unsigned char* m, UV nbytes)
}
+static int _is_trial_prime7(UV x)
+{
+ UV q, i;
+ i = 7;
+ while (1) { /* trial division, skipping multiples of 2/3/5 */
+ q = x/i; if (q<i) return 1; if (x==(q*i)) return 0; i += 4;
+ q = x/i; if (q<i) return 1; if (x==(q*i)) return 0; i += 2;
+ q = x/i; if (q<i) return 1; if (x==(q*i)) return 0; i += 4;
+ q = x/i; if (q<i) return 1; if (x==(q*i)) return 0; i += 2;
+ q = x/i; if (q<i) return 1; if (x==(q*i)) return 0; i += 4;
+ q = x/i; if (q<i) return 1; if (x==(q*i)) return 0; i += 6;
+ q = x/i; if (q<i) return 1; if (x==(q*i)) return 0; i += 2;
+ q = x/i; if (q<i) return 1; if (x==(q*i)) return 0; i += 6;
+ }
+ return 2;
+}
+
static int _is_prime7(UV x)
{
UV q, i;
+
+ if (x >= UVCONST(100000000))
+ return is_prob_prime(x); /* We know this works for all 64-bit n */
+
i = 7;
while (1) { /* trial division, skipping multiples of 2/3/5 */
q = x/i; if (q<i) return 1; if (x==(q*i)) return 0; i += 4;
@@ -96,16 +117,7 @@ int is_prime(UV n)
if (n <= get_prime_cache(0, &sieve))
return ((sieve[d] & mtab) == 0);
-#if 0
- /* Trial division, mod 30 */
return _is_prime7(n);
-#else
- if (n < UVCONST(100000000)) {
- return _is_prime7(n);
- } else {
- return is_prob_prime(n); /* We know this works for all 64-bit n */
- }
-#endif
}
@@ -125,7 +137,7 @@ UV next_trial_prime(UV n)
d = n/30;
m = n - d*30;
m = nextwheel30[m]; if (m == 1) d++;
- while (!_is_prime7(d*30+m)) {
+ while (!_is_trial_prime7(d*30+m)) {
m = nextwheel30[m]; if (m == 1) d++;
}
return(d*30+m);
--
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