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