This is an automated email from the git hooks/post-receive script. ppm-guest pushed a commit to annotated tag v0.04 in repository libmath-prime-util-perl.
commit f9b919393f4a49dce665a833b0f49f4faf754cd1 Author: Dana Jacobsen <[email protected]> Date: Thu Jun 7 16:53:16 2012 -0600 Factoring updates --- XS.xs | 37 +++++++++---- examples/bench-factor.pl | 112 +++++++++++++++++++++++++++++++++++++++ examples/bench-is-prime.pl | 9 ++-- factor.c | 129 +++++++++++++++++++++++++-------------------- factor.h | 35 +++--------- util.h | 2 +- 6 files changed, 223 insertions(+), 101 deletions(-) diff --git a/XS.xs b/XS.xs index f15eb6e..9c02513 100644 --- a/XS.xs +++ b/XS.xs @@ -230,20 +230,35 @@ factor(IN UV n) if (n < 4) { XPUSHs(sv_2mortal(newSVuv( n ))); /* If n is 0-3, we're done. */ } else { + UV tlim = 19; /* Below this we've checked */ UV factor_stack[MPU_MAX_FACTORS+1]; int nstack = 0; - /* Quick trial divisions. We could do tricky gcd magic here. */ + /* Quick trial divisions. Crude use of GCD to hopefully go faster. */ 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 ))); } - while ( (n% 7) == 0 ) { n /= 7; XPUSHs(sv_2mortal(newSVuv( 7 ))); } - 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 >= UVCONST(3*3)) && (gcd_ui(n, UVCONST(3234846615) != 1)) ) { + while ( (n% 3) == 0 ) { n /= 3; XPUSHs(sv_2mortal(newSVuv( 3 ))); } + while ( (n% 5) == 0 ) { n /= 5; XPUSHs(sv_2mortal(newSVuv( 5 ))); } + while ( (n% 7) == 0 ) { n /= 7; XPUSHs(sv_2mortal(newSVuv( 7 ))); } + 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 ))); } + while ( (n%19) == 0 ) { n /= 19; XPUSHs(sv_2mortal(newSVuv( 19 ))); } + while ( (n%23) == 0 ) { n /= 23; XPUSHs(sv_2mortal(newSVuv( 23 ))); } + while ( (n%29) == 0 ) { n /= 29; XPUSHs(sv_2mortal(newSVuv( 29 ))); } + tlim = 31; + } + if ( (n >= UVCONST(31*31)) && (gcd_ui(n, UVCONST(95041567) != 1)) ) { + while ( (n%31) == 0 ) { n /= 31; XPUSHs(sv_2mortal(newSVuv( 31 ))); } + while ( (n%37) == 0 ) { n /= 37; XPUSHs(sv_2mortal(newSVuv( 37 ))); } + while ( (n%41) == 0 ) { n /= 41; XPUSHs(sv_2mortal(newSVuv( 41 ))); } + while ( (n%43) == 0 ) { n /= 43; XPUSHs(sv_2mortal(newSVuv( 43 ))); } + while ( (n%47) == 0 ) { n /= 47; XPUSHs(sv_2mortal(newSVuv( 47 ))); } + tlim = 53; + } do { /* loop over each remaining factor */ - while ( (n >= (19*19)) && (!is_definitely_prime(n)) ) { + while ( (n >= (tlim*tlim)) && (!is_definitely_prime(n)) ) { int split_success = 0; - if (n > UVCONST(60000000) ) { /* tune this */ + if (n > UVCONST(10000000) ) { /* tune this */ /* For sufficiently large n, try more complex methods. */ /* SQUFOF (succeeds 98-99.9%) */ split_success = squfof_factor(n, factor_stack+nstack, 256*1024)-1; @@ -259,8 +274,8 @@ factor(IN UV n) n = factor_stack[nstack]; } else { /* trial divisions */ - UV f = 19; - UV m = 19; + UV f = tlim; + UV m = tlim % 30; UV limit = sqrt((double) n); while (f <= limit) { if ( (n%f) == 0 ) { diff --git a/examples/bench-factor.pl b/examples/bench-factor.pl new file mode 100755 index 0000000..f88404a --- /dev/null +++ b/examples/bench-factor.pl @@ -0,0 +1,112 @@ +#!/usr/bin/env perl +use strict; +use warnings; +use Math::Prime::Util qw/factor/; +use Math::Factor::XS qw/prime_factors/; +use Benchmark qw/:all/; +use List::Util qw/min max reduce/; +my $count = shift || -2; +my $is64bit = (~0 > 4294967295); +my $maxdigits = ($is64bit) ? 20 : 10; # Noting the range is limited for max. +my $semiprimes = 0; +my $howmany = 10000; + +srand(29); + +for my $d ( 3 .. $maxdigits ) { + print "Factor $howmany $d-digit numbers\n"; + test_at_digits($d, $howmany); +} + +sub test_at_digits { + my $digits = shift; + die "Digits has to be >= 1" unless $digits >= 1; + die "Digits has to be <= $maxdigits" if $digits > $maxdigits; + my $quantity = shift; + + my @rnd = genrand($digits, $quantity); + my @smp = gensemi($digits, $quantity); + + # verify + foreach my $p (@rnd, @smp) { + my @mpxs = prime_factors($p); push @mpxs, $p if $p < 2; + + verify_factor($p, \@mpxs, [factor($p)], "Math::Prime::Util $Math::Prime::Util::VERSION"); +} + + + #my $min_num = min @nums; + #my $max_num = max @nums; + #my $whatstr = "$digits-digit ", $semiprimes ? "semiprime" : "random"; + #print "factoring 1000 $digits-digit ", + # $semiprimes ? "semiprimes" : "random numbers", + # " ($min_num - $max_num)\n"; + + my $lref = { + "MPU rand" => sub { factor($_) for @rnd }, + "MPU semi" => sub { factor($_) for @smp }, + "MFXS rand" => sub { prime_factors($_) for @rnd }, + "MFXS semi" => sub { prime_factors($_) for @smp }, + }; + cmpthese($count, $lref); +} + +sub verify_factor { + my $n = shift; + my $aref_master = shift; + my $aref_check = shift; + my $name = shift; + + my @master = sort {$a<=>$b} @{$aref_master}; + my @check = sort {$a<=>$b} @{$aref_check}; + die "Factor $n master fail!" unless $n == reduce { $a * $b } @master; + die "Factor $n fail: $name" unless $#check == $#master; + die "Factor $n fail: $name" unless $n == reduce { $a * $b } @check; + for (0 .. $#master) { + die "Factor $n fail: $name" unless $master[$_] == $check[$_]; + } + 1; +} + +sub genrand { + my $digits = shift; + my $num = shift; + + 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 .. $num); + return @nums; +} + +sub gensemi { + my $digits = shift; + my $num = shift; + + my @min_factors_by_digit = (2,2,3,5,7,13,23,47,97); + my $smallest_factor = $min_factors_by_digit[$digits]; + $smallest_factor = $min_factors_by_digit[-1] unless defined $smallest_factor; + + my $base = ($digits == 1) ? 0 : int(10 ** ($digits-1)); + my $max = int(10 ** $digits); + $max = (~0-4) if $max > (~0-4); + my @semiprimes; + + foreach my $i (1 .. $num) { + my @factors; + my $n; + while (1) { + $n = $base + int(rand($max-$base)); + $n += 1 if ($n%2) == 0; + $n += 3 if ($n%3) == 0; + @factors = Math::Prime::Util::factor($n); + next if scalar @factors != 2; + next if $factors[0] < $smallest_factor; + next if $factors[1] < $smallest_factor; + last if scalar @factors == 2; + } + die "ummm... $n != $factors[0] * $factors[1]\n" unless $n == $factors[0] * $factors[1]; + push @semiprimes, $n; + } + return @semiprimes; +} diff --git a/examples/bench-is-prime.pl b/examples/bench-is-prime.pl index cc9694a..69feee5 100755 --- a/examples/bench-is-prime.pl +++ b/examples/bench-is-prime.pl @@ -4,14 +4,17 @@ use strict; #use Math::Primality; use Math::Prime::XS; use Math::Prime::Util; -use Math::Pari; +#use Math::Pari; #use Math::Prime::FastSieve; use Benchmark qw/:all/; use List::Util qw/min max/; my $count = shift || -5; +my $is64bit = (~0 > 4294967295); +my $maxdigits = ($is64bit) ? 20 : 10; # Noting the range is limited for max. + srand(29); -test_at_digits($_) for (5..15); +test_at_digits($_) for (5 .. $maxdigits); sub test_at_digits { @@ -36,7 +39,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 }, + #'Math::Pari' => sub { Math::Pari::isprime($_) for @nums }, }); print "\n"; } diff --git a/factor.c b/factor.c index 3ef5039..59e4314 100644 --- a/factor.c +++ b/factor.c @@ -3,6 +3,8 @@ #include <string.h> #include <assert.h> #include <limits.h> +#define __STDC_LIMIT_MACROS +#include <stdint.h> #include <math.h> #include "factor.h" @@ -77,73 +79,84 @@ int trial_factor(UV n, UV *factors, UV maxtrial) return nfactors; } -static UV gcd_ui(UV x, UV y) { - UV t; - if (y < x) { t = x; x = y; y = t; } +#if UINT64_MAX > UV_MAX - while (y > 0) { - x = x % y; - t = x; x = y; y = t; - } - 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 r = 0; - while (b > 0) { - if (b & 1) { - if (r == 0) { - r = a; - } else { - r = m - r; - r = (a >= r) ? a-r : m-r+a; - } - } - a = (a > (m-a)) ? (a-m)+a : a+a; - b >>= 1; - } - return r; -} + /* We have 64-bit available, but UV is 32-bit. Do the math in 64-bit. + * Even if it is emulated, it should be as fast or faster than us doing it. + */ + #define addmod(n,a,m) (UV)(((uint64_t)(n)+(uint64_t)(a)) % ((uint64_t)(m))) + #define mulmod(a,b,m) (UV)(((uint64_t)(a)*(uint64_t)(b)) % ((uint64_t)(m))) + #define sqrmod(n,m) (UV)(((uint64_t)(n)*(uint64_t)(n)) % ((uint64_t)(m))) -/* n^power mod m */ -static UV powmod(UV n, UV power, UV m) { - UV t = 1; - if (m < HALF_WORD) { - n %= m; - while (power) { - if (power & 1) - t = (t*n)%m; - n = (n*n)%m; - power >>= 1; - } - } else { + static UV powmod(UV n, UV power, UV m) { + UV t = 1; while (power) { if (power & 1) t = mulmod(t, n, m); - n = (n < HALF_WORD) ? (n*n)%m : mulmod(n, n, m); + n = sqrmod(n, m); power >>= 1; } + return t; } - return t; -} -/* n + a mod m */ -static UV addmod(UV n, UV a, UV m) { - return ((m-n) > a) ? n+a : n+a-m; -} +#else + + /* UV is the largest integral type available (that we know of). */ + + /* 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 r = 0; + while (b > 0) { + if (b & 1) { + if (r == 0) { + r = a; + } else { + r = m - r; + r = (a >= r) ? a-r : m-r+a; + } + } + a = (a > (m-a)) ? (a-m)+a : a+a; + b >>= 1; + } + return r; + } -/* n^2 mod m */ -#define powsqr(n, m) (n < HALF_WORD) ? (n*n)%m : mulmod(n,n,m) + #define addmod(n,a,m) ((((m)-(n)) > (a)) ? ((n)+(a)) : ((n)+(a)-(m))) + #define mulmod(a,b,m) (((a)|(b)) < HALF_WORD) ? ((a)*(b))%(m) : _mulmod(a,b,m) + #define sqrmod(n,m) ((n) < HALF_WORD) ? ((n)*(n))%(m) : _mulmod(n,n,m) + + /* n^power mod m */ + static UV powmod(UV n, UV power, UV m) { + UV t = 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 = sqrmod(n,m); + power >>= 1; + } + } + return t; + } + +#endif /* n^power + a mod m */ -#define powmodadd(n, p, a, m) addmod(powmod(n,p,m),a,m) +#define powaddmod(n, p, a, m) addmod(powmod(n,p,m),a,m) /* n^2 + a mod m */ -#define powsqradd(n, a, m) addmod(powsqr(n,m), a, m) +#define sqraddmod(n, a, m) addmod(sqrmod(n,m), a,m) /* Miller-Rabin probabilistic primality test @@ -175,7 +188,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 = powsqr(x, n); + x = sqrmod(x, n); if (x == 1) { return 0; } else if (x == (n-1)) { @@ -323,7 +336,7 @@ int holf_factor(UV n, UV *factors, UV rounds) for (i = 1; i <= rounds; i++) { s = sqrt(n*i); /* TODO: overflow here */ if ( (s*s) != (n*i) ) s++; - m = powsqr(s, n); + m = sqrmod(s, n); /* Cheaper would be: * if (m is probably not a perfect sequare) continue; * f = sqrt(m); @@ -370,7 +383,7 @@ int pbrent_factor(UV n, UV *factors, UV rounds) } for (i = 1; i <= rounds; i++) { - Xi = powsqradd(Xi, a, n); + Xi = sqraddmod(Xi, a, n); f = gcd_ui( (Xi>Xm) ? Xi-Xm : Xm-Xi, n); if ( (f != 1) && (f != n) ) { factors[0] = f; @@ -408,9 +421,9 @@ int prho_factor(UV n, UV *factors, UV rounds) } for (i = 1; i <= rounds; i++) { - U = powsqradd(U, a, n); - V = powsqradd(V, a, n); - V = powsqradd(V, a, n); + U = sqraddmod(U, a, n); + V = sqraddmod(V, a, n); + V = sqraddmod(V, a, n); f = gcd_ui( (U > V) ? U-V : V-U, n); if (f == n) { if (in_loop++) /* Mark that we've been here */ diff --git a/factor.h b/factor.h index 57cedb0..ad0a757 100644 --- a/factor.h +++ b/factor.h @@ -22,34 +22,13 @@ extern int pminus1_factor(UV n, UV *factors, UV maxrounds); extern int miller_rabin(UV n, const UV *bases, UV nbases); extern int is_prob_prime(UV n); -#if 0 -/* Try to make a quick probable prime test */ -static int is_perhaps_prime(UV n) -{ - static const UV bases2[2] = {31, 73}; - static const UV bases3[3] = {2,7,61}; - if ( (n == 2) || (n == 3) || (n == 5) || (n == 7) ) - return 2; - if ( (n<2) || ((n% 2)==0) || ((n% 3)==0) || ((n% 5)==0) || ((n% 7)==0) ) - return 0; - if (n < 121) - return 2; - if (n < UVCONST(9080191)) - return 2*miller_rabin(n, bases2, 2); - else - return miller_rabin(n, bases3, 3) * ((n < UVCONST(4759123141)) ? 2 : 1); +static UV gcd_ui(UV x, UV y) { + UV t; + if (y < x) { t = x; x = y; y = t; } + while (y > 0) { + t = y; y = x % y; x = t; /* y1 <- x0 % y0 ; x1 <- y0 */ + } + return x; } -static int is_perhaps_prime7(UV n) -{ - static const UV bases2[2] = {31, 73}; - static const UV bases3[3] = {2,7,61}; - /* n must be >= 73 */ - if (n < UVCONST(9080191)) - return 2*miller_rabin(n, bases2, 2); - else - return miller_rabin(n, bases3, 3) * ((n < UVCONST(4759123141)) ? 2 : 1); -} -#endif - #endif diff --git a/util.h b/util.h index b74d8af..0ecd27c 100644 --- a/util.h +++ b/util.h @@ -25,7 +25,7 @@ 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 +#if (BITS_PER_WORD == 32) && (UINT64_MAX <= UV_MAX) #define MPU_PROB_PRIME_BEST UVCONST(100000000) #else #define MPU_PROB_PRIME_BEST UVCONST(100000) -- 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 [email protected] http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-perl-cvs-commits
