This is an automated email from the git hooks/post-receive script. ppm-guest pushed a commit to annotated tag v0.29 in repository libmath-prime-util-perl.
commit 70ba979e592fbe95a3e2f7f79003599c22dd1ae0 Author: Dana Jacobsen <d...@acm.org> Date: Tue May 28 20:57:57 2013 -0700 Don't use MULTICALL yet -- memory leak --- Changes | 3 ++ XS.xs | 70 +++++++++++++++++++------------- factor.c | 120 ++++++++++++++++++++++++++++++++++++++++++++++++------- factor.h | 16 ++++---- t/32-iterators.t | 13 ++++++ 5 files changed, 172 insertions(+), 50 deletions(-) diff --git a/Changes b/Changes index ba643c7..38bbf9f 100644 --- a/Changes +++ b/Changes @@ -4,6 +4,9 @@ Revision history for Perl extension Math::Prime::Util. - Fix a signed vs. unsigned char issue in ranged moebius. + - forprimes now uses a segmented sieve. This (1) allows arbitrary 64-bit + ranges with good memory use, and (2) allows nesting on threaded perls. + - Added XS strong Lucas test, but not using Selfridge parameters. Hence we only use it internally. diff --git a/XS.xs b/XS.xs index 0bbf285..2d23cd5 100644 --- a/XS.xs +++ b/XS.xs @@ -399,6 +399,9 @@ pminus1_factor(IN UV n, IN UV B1 = 1*1024*1024, IN UV B2 = 0) SIMPLE_FACTOR(pminus1_factor, n, B1, B2); int +_XS_is_pseudoprime(IN UV n, IN UV a) + +int _XS_miller_rabin(IN UV n, ...) PREINIT: UV bases[64]; @@ -429,6 +432,9 @@ int _XS_is_strong_lucas_pseudoprime(IN UV n) int +_XS_is_frobenius_underwood_pseudoprime(IN UV n) + +int _XS_is_prob_prime(IN UV n) int @@ -681,39 +687,49 @@ forprimes (SV* block, IN SV* svbeg, IN SV* svend = 0) croak("Not a subroutine reference"); SAVESPTR(GvSV(PL_defgv)); svarg = newSVuv(0); - ctx = start_segment_primes(beg, end, &segment); - if (!CvISXSUB(cv)) { - dMULTICALL; - I32 gimme = G_VOID; - PUSH_MULTICALL(cv); - while (beg < 7) { - beg = (beg <= 2) ? 2 : (beg <= 3) ? 3 : 5; - { sv_setuv(svarg, beg); GvSV(PL_defgv) = svarg; MULTICALL; } - beg += 1 + (beg > 2); + /* Handle early part */ + while (beg < 7) { + dSP; + beg = (beg <= 2) ? 2 : (beg <= 3) ? 3 : 5; + if (beg <= end) { + sv_setuv(svarg, beg); + GvSV(PL_defgv) = svarg; + PUSHMARK(SP); + call_sv((SV*)cv, G_VOID|G_DISCARD); } + beg += 1 + (beg > 2); + } + if (beg <= end) { + ctx = start_segment_primes(beg, end, &segment); while (next_segment_primes(ctx, &seg_base, &seg_low, &seg_high)) { - START_DO_FOR_EACH_SIEVE_PRIME( segment, seg_low - seg_base, seg_high - seg_base ) - { sv_setuv(svarg, seg_base + p); GvSV(PL_defgv) = svarg; MULTICALL; } - END_DO_FOR_EACH_SIEVE_PRIME - } + /* I'm getting a memory leak in the MULTICALL and I'm having no luck + * finding out why it is happening. Forget MULTICALL for now. */ + if (0 && !CvISXSUB(cv)) { + dMULTICALL; + I32 gimme = G_VOID; + PUSH_MULTICALL(cv); + START_DO_FOR_EACH_SIEVE_PRIME( segment, seg_low - seg_base, seg_high - seg_base ) { + sv_setuv(svarg, seg_base + p); + GvSV(PL_defgv) = svarg; + MULTICALL; + } END_DO_FOR_EACH_SIEVE_PRIME #ifdef PERL_HAS_BAD_MULTICALL_REFCOUNT - if (CvDEPTH(multicall_cv) > 1) - SvREFCNT_inc(multicall_cv); + if (CvDEPTH(multicall_cv) > 1) + SvREFCNT_inc(multicall_cv); #endif - POP_MULTICALL; - } else { - while (beg < 7) { - beg = (beg <= 2) ? 2 : (beg <= 3) ? 3 : 5; - { dSP; sv_setuv(svarg, beg); GvSV(PL_defgv) = svarg; PUSHMARK(SP); call_sv((SV*)cv, G_VOID|G_DISCARD); } - beg += 1 + (beg > 2); - } - while (next_segment_primes(ctx, &seg_base, &seg_low, &seg_high)) { - START_DO_FOR_EACH_SIEVE_PRIME( segment, seg_low - seg_base, seg_high - seg_base ) - { dSP; sv_setuv(svarg, seg_base + p); GvSV(PL_defgv) = svarg; PUSHMARK(SP); call_sv((SV*)cv, G_VOID|G_DISCARD); } - END_DO_FOR_EACH_SIEVE_PRIME + POP_MULTICALL; + } else { + START_DO_FOR_EACH_SIEVE_PRIME( segment, seg_low - seg_base, seg_high - seg_base ) { + dSP; + sv_setuv(svarg, seg_base + p); + GvSV(PL_defgv) = svarg; + PUSHMARK(SP); + call_sv((SV*)cv, G_VOID|G_DISCARD); + } END_DO_FOR_EACH_SIEVE_PRIME + } } + end_segment_primes(ctx); } - end_segment_primes(ctx); SvREFCNT_dec(svarg); #endif XSRETURN_UNDEF; diff --git a/factor.c b/factor.c index 6533ef7..db8bfba 100644 --- a/factor.c +++ b/factor.c @@ -308,6 +308,25 @@ static int jacobi_iu(IV in, UV m) { } +/* Fermat pseudoprime */ +int _XS_is_pseudoprime(UV n, UV a) +{ + UV x; + UV const nm1 = n-1; + + if (n <= 3) return 0; + if (a < 2) croak("Base %"UVuf" is invalid", a); + if (a >= n) { + a %= n; + if ( a <= 1 || a == nm1 ) + return 1; + } + + x = powmod(a, nm1, n); /* x = a^(n-1) mod n */ + return (x == 1); +} + + /* Miller-Rabin probabilistic primality test * Returns 1 if probably prime relative to the bases, 0 if composite. * Bases must be between 2 and n-2 @@ -405,28 +424,31 @@ int _XS_is_prob_prime(UV n) #else -#if 0 - /* TODO: Must verify this Lucas with Feitsma database */ - if (1 && n >= UVCONST(4294967295)) { + /* Verified with Feitsma database. No counterexamples below 2^64. + * This is faster than multiple M-R routines once we're over 32-bit */ + if (n >= UVCONST(4294967295)) { prob_prime = _SPRP2(n) && _XS_is_strong_lucas_pseudoprime(n); return 2*prob_prime; } -#endif - /* Better bases from http://miller-rabin.appspot.com/, 23 May 2013 */ + /* Better bases from http://miller-rabin.appspot.com/, 27 May 2013 */ + /* Verify with: + * cat /local/share/spsps-below-2-to-64.txt | perl -MMath::Prime::Util=:all + * -nE 'chomp; next unless is_strong_pseudoprime($_, @bases); say;' + */ if (n < UVCONST(341531)) { - bases[0] = UVCONST(9345883071009581737); + bases[0] = UVCONST( 9345883071009581737 ); nbases = 1; - } else if (n < UVCONST(716169301)) { - bases[0] = 15; - bases[1] = UVCONST( 13393019396194701 ); + } else if (n < UVCONST(885594169)) { + bases[0] = UVCONST( 725270293939359937 ); + bases[1] = UVCONST( 3569819667048198375 ); nbases = 2; - } else if (n < UVCONST(273919523041)) { /* 29+ bits */ - bases[0] = 15; - bases[1] = UVCONST( 7363882082 ); - bases[2] = UVCONST( 992620450144556 ); + } else if (n < UVCONST(350269456337)) { + bases[0] = UVCONST( 4230279247111683200 ); + bases[1] = UVCONST( 14694767155120705706 ); + bases[2] = UVCONST( 16641139526367750375 ); nbases = 3; - } else if (n < UVCONST(55245642489451)) { /* 37+ bits */ + } else if (n < UVCONST(55245642489451)) { /* 38+ bits */ bases[0] = 2; bases[1] = UVCONST( 141889084524735 ); bases[2] = UVCONST( 1199124725622454117 ); @@ -1158,3 +1180,73 @@ int racing_squfof_factor(UV n, UV *factors, UV rounds) factors[0] = n; return 1; } + + +/****************************************************************************/ + +int _XS_is_frobenius_underwood_pseudoprime(UV n) +{ + int bit; + UV x, result, multiplier, a, b, np1, len, t1, t2, na; + IV t; + + if (n < 2) return 0; + if (n < 4) return 1; + if ((n % 2) == 0) return 0; + if (is_perfect_square(n,0)) return 0; + if (n == UV_MAX) return 0; + + x = 0; + t = -1; + while ( jacobi_iu( t, n ) != -1 ) { + x++; + t = (IV)(x*x) - 4; + } + result = addmod( addmod(x, x, n), 5, n); + multiplier = addmod(x, 2, n); + + a = 1; + b = 2; + np1 = n+1; + { UV v = np1; len = 1; while (v >>= 1) len++; } + + if (x == 0) { + for (bit = len-2; bit >= 0; bit--) { + t2 = addmod(b, b, n); + na = mulmod(a, t2, n); + t1 = addmod(b, a, n); + t2 = addmod(b, n-a, n); /* subtract */ + b = mulmod(t1, t2, n); + a = na; + if ( (np1 >> bit) & UVCONST(1) ) { + t1 = mulmod(a, 2, n); + na = addmod(t1, b, n); + t1 = addmod(b, b, n); + b = addmod(t1, n-a, n); /* subtract */ + a = na; + } + } + } else { + for (bit = len-2; bit >= 0; bit--) { + t1 = mulmod(a, x, n); + t2 = addmod(b, b, n); + t1 = addmod(t1, t2, n); + na = mulmod(a, t1, n); + t1 = addmod(b, a, n); + t2 = addmod(b, n-a, n); /* subtract */ + b = mulmod(t1, t2, n); + a = na; + if ( (np1 >> bit) & UVCONST(1) ) { + t1 = mulmod(a, multiplier, n); + na = addmod(t1, b, n); + t1 = addmod(b, b, n); + b = addmod(t1, n-a, n); /* subtract */ + a = na; + } + } + } + if (_XS_get_verbose()>1) printf("%"UVuf" is %s with x = %"UVuf"\n", n, (a == 0 && b == result) ? "probably prime" : "composite", x); + if (a == 0 && b == result) + return 1; + return 0; +} diff --git a/factor.h b/factor.h index 77608c5..29722bf 100644 --- a/factor.h +++ b/factor.h @@ -10,22 +10,20 @@ extern int factor(UV n, UV *factors); 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 racing_squfof_factor(UV n, UV *factors, UV rounds); - - extern int pbrent_factor(UV n, UV *factors, UV maxrounds, UV a); - extern int prho_factor(UV n, UV *factors, UV maxrounds); - extern int pminus1_factor(UV n, UV *factors, UV B1, UV B2); +extern int squfof_factor(UV n, UV *factors, UV rounds); +extern int racing_squfof_factor(UV n, UV *factors, UV rounds); + +extern int _XS_is_pseudoprime(UV n, UV a); extern int _XS_miller_rabin(UV n, const UV *bases, int nbases); +extern int _SPRP2(UV n); extern int _XS_is_prob_prime(UV n); -extern int _XS_is_prime_tom(UV n, int t); +extern int _XS_is_strong_lucas_pseudoprime(UV n); +extern int _XS_is_frobenius_underwood_pseudoprime(UV n); extern UV _XS_divisor_sum(UV n); diff --git a/t/32-iterators.t b/t/32-iterators.t index bfc285b..7e4ae1e 100644 --- a/t/32-iterators.t +++ b/t/32-iterators.t @@ -8,6 +8,7 @@ use Math::Prime::Util qw/primes forprimes prime_iterator/; my $use64 = Math::Prime::Util::prime_get_config->{'maxbits'} > 32; plan tests => 8 + 5 + + 11 + 3 + 7 + 2; @@ -22,6 +23,18 @@ ok(!eval { forprimes { 1 } "abc"; }, "forprimes abc"); ok(!eval { forprimes { 1 } 2, "abc"; }, "forprimes 2, abc"); ok(!eval { forprimes { 1 } 5.6; }, "forprimes abc"); +{my @t; forprimes {push @t,$_} 1; is_deeply( [@t], [], "forprimes 1" ); } +{my @t; forprimes {push @t,$_} 2; is_deeply( [@t], [2], "forprimes 3" ); } +{my @t; forprimes {push @t,$_} 3; is_deeply( [@t], [2,3], "forprimes 3" ); } +{my @t; forprimes {push @t,$_} 4; is_deeply( [@t], [2,3], "forprimes 4" ); } +{my @t; forprimes {push @t,$_} 5; is_deeply( [@t], [2,3,5], "forprimes 5" ); } +{my @t; forprimes {push @t,$_} 3,5; is_deeply( [@t], [3,5], "forprimes 3,5" ); } +{my @t; forprimes {push @t,$_} 3,6; is_deeply( [@t], [3,5], "forprimes 3,6" ); } +{my @t; forprimes {push @t,$_} 3,7; is_deeply( [@t], [3,5,7], "forprimes 3,7" ); } +{my @t; forprimes {push @t,$_} 5,7; is_deeply( [@t], [5,7], "forprimes 5,7" ); } +{my @t; forprimes {push @t,$_} 5,11; is_deeply( [@t], [5,7,11], "forprimes 5,11" ); } +{my @t; forprimes {push @t,$_} 7,11; is_deeply( [@t], [7,11], "forprimes 7,11" ); } + { my @t; forprimes { push @t, $_ } 50; is_deeply( [@t], [2,3,5,7,11,13,17,19,23,29,31,37,41,43,47], "forprimes 50" ); } -- 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