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
