[libmath-prime-util-perl] 08/11: Factoring is isprime updates
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 52935e2c6616e4c7e35e74e4621d9050fb78216b
Author: Dana Jacobsen <dana at acm.org>
Date: Wed Jun 6 17:30:28 2012 -0600
Factoring is isprime updates
---
XS.xs | 10 +--
examples/bench-factor-semiprime.pl | 8 +-
examples/bench-is-prime.pl | 18 +++--
factor.c | 156 ++++++++++++++++---------------------
util.c | 36 ++++++++-
util.h | 4 +
6 files changed, 128 insertions(+), 104 deletions(-)
diff --git a/XS.xs b/XS.xs
index 0115743..53c1a5f 100644
--- a/XS.xs
+++ b/XS.xs
@@ -242,16 +242,15 @@ factor(IN UV n)
while ( (n%13) == 0 ) { n /= 13; XPUSHs(sv_2mortal(newSVuv( 13 ))); }
while ( (n%17) == 0 ) { n /= 17; XPUSHs(sv_2mortal(newSVuv( 17 ))); }
do { /* loop over each remaining factor */
- while ( (n >= (19*19)) && (!is_prime(n)) ) {
- /* n is composite, so factor it. */
+ while ( (n >= (19*19)) && (!is_definitely_prime(n)) ) {
int split_success = 0;
- if (n > UVCONST(10000000) ) { /* tune this */
- /* For sufficiently large numbers, try more complex methods. */
+ 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;
assert( (split_success == 0) || (split_success == 1) );
/* a few rounds of Pollard rho (succeeds 99+% of the rest) */
- if (!split_success) {
+ if (1 && !split_success) {
split_success = prho_factor(n, factor_stack+nstack, 400)-1;
assert( (split_success == 0) || (split_success == 1) );
}
@@ -279,6 +278,7 @@ factor(IN UV n)
f += wheeladvance30[m];
m = nextwheel30[m];
}
+ break; /* We just factored n via trial division. Exit loop. */
}
}
if (n != 1) XPUSHs(sv_2mortal(newSVuv( n )));
diff --git a/examples/bench-factor-semiprime.pl b/examples/bench-factor-semiprime.pl
index 70f28e8..e78a20c 100755
--- a/examples/bench-factor-semiprime.pl
+++ b/examples/bench-factor-semiprime.pl
@@ -2,6 +2,7 @@
use strict;
use warnings;
$| = 1; # fast pipes
+srand(377);
use Math::Prime::Util qw/factor/;
use Math::Factor::XS qw/prime_factors/;
@@ -12,8 +13,11 @@ my $count = shift || -2;
my @min_factors_by_digit = (2,2,3,3,5,11,17,47,97);
my $smallest_factor_allowed = $min_factors_by_digit[$digits];
$smallest_factor_allowed = $min_factors_by_digit[-1] unless defined $smallest_factor_allowed;
-my $numprimes = 50;
-die "Digits has to be >= 2 and <= 19" unless $digits >= 2 && $digits <= 19;
+my $numprimes = 100;
+
+die "Digits has to be >= 2" unless $digits >= 2;
+die "Digits has to be <= 10" if (~0 == 4294967295) && ($digits > 10);
+die "Digits has to be <= 19" if $digits > 19;
# Construct some semiprimes of the appropriate number of digits
# There are much cleverer ways of doing this, using randomly selected
diff --git a/examples/bench-is-prime.pl b/examples/bench-is-prime.pl
index 8dff8f7..82382fc 100755
--- a/examples/bench-is-prime.pl
+++ b/examples/bench-is-prime.pl
@@ -9,23 +9,25 @@ use Benchmark qw/:all/;
use List::Util qw/min max/;
my $count = shift || -5;
-test_at_mag($_) for (4, 8, 12, 16);
+srand(29);
+test_at_digits($_) for (5..10);
-sub test_at_mag {
- my $magnitude = shift;
+sub test_at_digits {
+ my $digits = shift;
+ die "Digits must be > 0" unless $digits > 0;
- my $base = 10 ** ($magnitude-1);
- my $max = 10 ** $magnitude;
- my $digits = $magnitude+1;
- my @nums = map { $base+int(rand($max-$base)) } (1 .. 200);
+ my $base = ($digits == 1) ? 0 : int(10 ** ($digits-1));
+ my $max = int(10 ** $digits);
+ $max = ~0 if $max > ~0;
+ my @nums = map { $base+int(rand($max-$base)) } (1 .. 1000);
my $min_num = min @nums;
my $max_num = max @nums;
#my $sieve = Math::Prime::FastSieve::Sieve->new(10 ** $magnitude + 1);
#Math::Prime::Util::prime_precalc(10 ** $magnitude + 1);
- print "is_prime for random $digits-digit numbers ($min_num - $max_num)\n";
+ print "is_prime for 1000 random $digits-digit numbers ($min_num - $max_num)\n";
cmpthese($count,{
#'Math::Primality' => sub { Math::Primality::is_prime($_) for @nums },
diff --git a/factor.c b/factor.c
index afbc7db..2dfc9c0 100644
--- a/factor.c
+++ b/factor.c
@@ -89,6 +89,9 @@ static UV gcd_ui(UV x, UV y) {
return x;
}
+/* if n is smaller than this, you can multiply without overflow */
+#define HALF_WORD (UVCONST(1) << (BITS_PER_WORD/2))
+
static UV mulmod(UV a, UV b, UV m) {
UV p;
UV r = 0;
@@ -107,31 +110,36 @@ static UV mulmod(UV a, UV b, UV m) {
return r;
}
-/* n^power + a mod m */
-static UV powmodadd(UV n, UV power, UV add, UV m) {
- UV t = 1;
- while (power) {
- if (power & 1)
- t = mulmod(t, n, m);
- n = mulmod(n, n, m);
- power >>= 1;
- }
- t = ((m-t) > add) ? t+add : t+add-m; /* (t+a) % m where t < m */
- return t;
-}
-
/* n^power mod m */
static UV powmod(UV n, UV power, UV m) {
UV t = 1;
- while (power) {
- if (power & 1)
- t = mulmod(t, n, m);
- n = mulmod(n, n, m);
- power >>= 1;
+ if (m < HALF_WORD) {
+ n %= m;
+ while (power) {
+ if (power & 1)
+ t = (t*n)%m;
+ n = (n*n)%m;
+ power >>= 1;
+ }
+ } else {
+ while (power) {
+ if (power & 1)
+ t = mulmod(t, n, m);
+ n = (n < HALF_WORD) ? (n*n)%m : mulmod(n, n, m);
+ power >>= 1;
+ }
}
return t;
}
+/* n + a mod m */
+static UV addmod(UV n, UV a, UV m) {
+ return ((m-n) > a) ? n+a : n+a-m;
+}
+
+/* n^power + a mod m */
+#define powmodadd(n, p, a, m) addmod(powmod(n,p,m),a,m)
+
/* Miller-Rabin probabilistic primality test
* Returns 1 if probably prime relative to the bases, 0 if composite.
@@ -162,7 +170,7 @@ int miller_rabin(UV n, const UV *bases, UV nbases)
if ( (x == 1) || (x == (n-1)) ) continue;
for (r = 0; r < s; r++) {
- x = powmod(x, 2, n);
+ x = (x < HALF_WORD) ? (x*x) % n : mulmod(x, x, n);
if (x == 1) {
return 0;
} else if (x == (n-1)) {
@@ -270,7 +278,6 @@ int is_prob_prime(UV n)
*/
int fermat_factor(UV n, UV *factors, UV rounds)
{
- int nfactors = 0;
IV sqn, x, y, r;
assert( (n >= 3) && ((n%2) != 0) );
@@ -289,12 +296,13 @@ int fermat_factor(UV n, UV *factors, UV rounds)
} while (r > 0);
}
r = (x-y)/2;
- if (r != 1)
- factors[nfactors++] = r;
- n /= r;
- if (n != 1)
- factors[nfactors++] = n;
- return nfactors;
+ if ( (r != 1) && (r != n) ) {
+ factors[0] = r;
+ factors[1] = n/r;
+ return 2;
+ }
+ factors[0] = n;
+ return 1;
}
/* Pollard / Brent
@@ -304,13 +312,12 @@ int fermat_factor(UV n, UV *factors, UV rounds)
*/
int pbrent_factor(UV n, UV *factors, UV rounds)
{
- int nfactors = 0;
- UV a, f, Xi, Xm, i;
+ UV a, f, i;
+ UV Xi = 2;
+ UV Xm = 2;
assert( (n >= 3) && ((n%2) != 0) );
- Xi = 2;
- Xm = 2;
switch (n%4) {
case 0: a = 1; break;
case 1: a = 3; break;
@@ -323,15 +330,15 @@ int pbrent_factor(UV n, UV *factors, UV rounds)
Xi = powmodadd(Xi, 2, a, n);
f = gcd_ui(Xi - Xm, n);
if ( (f != 1) && (f != n) ) {
- factors[nfactors++] = f;
- factors[nfactors++] = n/f;
- return nfactors;
+ factors[0] = f;
+ factors[1] = n/f;
+ return 2;
}
if ( (i & (i-1)) == 0) /* i is a power of 2 */
Xm = Xi;
}
- factors[nfactors++] = n;
- return nfactors;
+ factors[0] = n;
+ return 1;
}
/* Pollard's Rho
@@ -341,8 +348,9 @@ int pbrent_factor(UV n, UV *factors, UV rounds)
*/
int prho_factor(UV n, UV *factors, UV rounds)
{
- int nfactors = 0;
- UV a, f, t, U, V, i;
+ UV a, f, i;
+ UV U = 7;
+ UV V = 7;
assert( (n >= 3) && ((n%2) != 0) );
@@ -354,9 +362,6 @@ int prho_factor(UV n, UV *factors, UV rounds)
default: a = 3; break;
}
- U = 7;
- V = 7;
-
for (i = 1; i < rounds; i++) {
U = powmodadd(U, 2, a, n);
V = powmodadd(V, 2, a, n);
@@ -364,13 +369,13 @@ int prho_factor(UV n, UV *factors, UV rounds)
f = gcd_ui( (U > V) ? U-V : V-U, n);
if ( (f != 1) && (f != n) ) {
- factors[nfactors++] = f;
- factors[nfactors++] = n/f;
- return nfactors;
+ factors[0] = f;
+ factors[1] = n/f;
+ return 2;
}
}
- factors[nfactors++] = n;
- return nfactors;
+ factors[0] = n;
+ return 1;
}
/* Pollard's P-1
@@ -380,26 +385,22 @@ int prho_factor(UV n, UV *factors, UV rounds)
*/
int pminus1_factor(UV n, UV *factors, UV rounds)
{
- int nfactors = 0;
- UV f, b, i;
+ UV f, i;
+ UV b = 13;
assert( (n >= 3) && ((n%2) != 0) );
- b = 13;
for (i = 1; i < rounds; i++) {
b = powmod(b, i, n);
f = gcd_ui(b-1, n);
- if (f == n) {
- factors[nfactors++] = n;
- return nfactors;
- } else if (f != 1) {
- factors[nfactors++] = f;
- factors[nfactors++] = n/f;
- return nfactors;
- }
+ if (n == 1) continue;
+ if (f == n) break; /* or we could continue */
+ factors[0] = f;
+ factors[1] = n/f;
+ return 2;
}
- factors[nfactors++] = n;
- return nfactors;
+ factors[0] = f;
+ return 1;
}
/* My modification of Ben Buhrow's modification of Bob Silverman's SQUFOF code.
@@ -413,7 +414,6 @@ static void enqu(IV q, IV *iter) {
int squfof_factor(UV n, UV *factors, UV rounds)
{
- int nfactors = 0;
UV rounds2 = rounds/16;
UV temp;
IV iq,ll,l2,p,pnext,q,qlast,r,s,t,i;
@@ -430,9 +430,9 @@ int squfof_factor(UV n, UV *factors, UV rounds)
p = s;
temp = n - (s*s); /* temp = n - floor(sqrt(n))^2 */
if (temp == 0) {
- factors[nfactors++] = s;
- factors[nfactors++] = s;
- return nfactors;
+ factors[0] = s;
+ factors[1] = s;
+ return 2;
}
q = temp; /* q = excess of n over next smaller square */
@@ -453,8 +453,7 @@ int squfof_factor(UV n, UV *factors, UV rounds)
else if (q <= l2)
enqu(q,&jter);
if (jter < 0) {
- factors[nfactors++] = n;
- return nfactors;
+ factors[0] = n; return 1;
}
}
@@ -477,8 +476,7 @@ int squfof_factor(UV n, UV *factors, UV rounds)
} /* end of main loop */
if (jter >= rounds) {
- factors[nfactors++] = n;
- return nfactors;
+ factors[0] = n; return 1;
}
qlast = r;
@@ -516,38 +514,22 @@ int squfof_factor(UV n, UV *factors, UV rounds)
}
if (iter >= rounds2) {
- factors[nfactors++] = n;
- return nfactors;
+ factors[0] = n; return 1;
}
if ((q & 1) == 0) q/=2; /* q was factor or 2*factor */
if ( (q == 1) || (q == n) ) {
- factors[nfactors++] = n;
- return nfactors;
+ factors[0] = n; return 1;
}
p = n/q;
- /* smallest factor first */
- if (q < p) { t = p; p = q; q = t; }
/* printf(" squfof found %lu = %lu * %lu in %ld/%ld rounds\n", n, p, q, jter, iter); */
-#if 1
- factors[nfactors++] = p;
- factors[nfactors++] = q;
-#else
- if ( (p >= refactor_above) && (is_prob_prime(p) == 0) )
- nfactors += squfof_factor(p, factors+nfactors, rounds, refactor_above);
- else
- factors[nfactors++] = p;
-
- if ( (q >= refactor_above) && (is_prob_prime(q) == 0) )
- nfactors += squfof_factor(q, factors+nfactors, rounds, refactor_above);
- else
- factors[nfactors++] = q;
-#endif
- return nfactors;
+ factors[0] = p;
+ factors[1] = q;
+ return 2;
}
diff --git a/util.c b/util.c
index 1295d94..44ebf7f 100644
--- a/util.c
+++ b/util.c
@@ -46,6 +46,7 @@ static UV count_zero_bits(const unsigned char* m, UV nbytes)
}
+/* Does trial division, assuming x not divisible by 2, 3, or 5 */
static int _is_trial_prime7(UV x)
{
UV q, i;
@@ -63,11 +64,12 @@ static int _is_trial_prime7(UV x)
return 2;
}
+/* Does trial division or prob tests, assuming x not divisible by 2, 3, or 5 */
static int _is_prime7(UV x)
{
UV q, i;
- if (x >= UVCONST(100000000))
+ if (x > MPU_PROB_PRIME_BEST)
return is_prob_prime(x); /* We know this works for all 64-bit n */
i = 7;
@@ -97,6 +99,7 @@ static const unsigned char prime_is_small[] =
0x00,0x20,0x8a,0x00,0x20,0x8a,0x00,0x00,0x88,0x80,0x00,0x02,0x22,0x08,0x02};
#define NPRIME_IS_SMALL (sizeof(prime_is_small)/sizeof(prime_is_small[0]))
+/* Return of 2 if n is prime, 0 if not. Do it fast. */
int is_prime(UV n)
{
UV d, m;
@@ -104,6 +107,32 @@ int is_prime(UV n)
const unsigned char* sieve;
if ( n < (NPRIME_IS_SMALL*8))
+ return ((prime_is_small[n/8] >> (n%8)) & 1) ? 2 : 0;
+
+ d = n/30;
+ m = n - d*30;
+ mtab = masktab30[m]; /* Bitmask in mod30 wheel */
+
+ /* Return 0 if a multiple of 2, 3, or 5 */
+ if (mtab == 0)
+ return 0;
+
+ if (n <= get_prime_cache(0, &sieve))
+ return ((sieve[d] & mtab) == 0) ? 2 : 0;
+
+ return _is_prime7(n);
+}
+
+/* Shortcut, asking for a very quick response of 1 = prime, 0 = dunno.
+ * No trial divisions will be done, making this useful for factoring.
+ */
+int is_definitely_prime(UV n)
+{
+ UV d, m;
+ unsigned char mtab;
+ const unsigned char* sieve;
+
+ if ( n < (NPRIME_IS_SMALL*8))
return ((prime_is_small[n/8] >> (n%8)) & 1);
d = n/30;
@@ -117,7 +146,10 @@ int is_prime(UV n)
if (n <= get_prime_cache(0, &sieve))
return ((sieve[d] & mtab) == 0);
- return _is_prime7(n);
+ if (n > MPU_PROB_PRIME_BEST)
+ return (is_prob_prime(n) == 2);
+
+ return 0;
}
diff --git a/util.h b/util.h
index 1a89605..9737e8c 100644
--- a/util.h
+++ b/util.h
@@ -4,6 +4,7 @@
#include "ptypes.h"
extern int is_prime(UV x);
+extern int is_definitely_prime(UV x);
extern UV next_trial_prime(UV x);
extern UV next_prime(UV x);
extern UV prev_prime(UV x);
@@ -22,4 +23,7 @@ extern UV nth_prime(UV x);
extern unsigned char* get_prime_segment(void);
extern void free_prime_segment(void);
+/* Above this value, is_prime will do deterministic Miller-Rabin */
+#define MPU_PROB_PRIME_BEST UVCONST(100000000)
+
#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