[libmath-prime-util-perl] 16/23: Add PP Lehmer prime count, including nth_prime speedup. Fix some no-XS-with-GMP issues.
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:45:56 UTC 2015
This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.14
in repository libmath-prime-util-perl.
commit a49ae24c7b0ad32dbe9b52b8eea5ad4459b84412
Author: Dana Jacobsen <dana at acm.org>
Date: Thu Nov 29 02:30:59 2012 -0800
Add PP Lehmer prime count, including nth_prime speedup. Fix some no-XS-with-GMP issues.
---
Changes | 3 +
TODO | 6 +-
lib/Math/Prime/Util.pm | 19 ++++--
lib/Math/Prime/Util/PP.pm | 162 +++++++++++++++++++++++++++++++++++++++-------
mulmod.h | 7 ++
t/14-nthprime.t | 2 +-
6 files changed, 164 insertions(+), 35 deletions(-)
diff --git a/Changes b/Changes
index 86c6544..f309135 100644
--- a/Changes
+++ b/Changes
@@ -19,6 +19,9 @@ Revision history for Perl extension Math::Prime::Util.
- Allow environment variables MPU_NO_XS and MPU_NO_GMP to turn off XS and
GMP support respectively if they are defined and equal to 1.
+ - Lehmer prime count for Pure Perl code, including use in nth_prime.
+ Almost 10x speedup on test suite if using only Pure Perl.
+
0.13 19 November 2012
- Fix an issue with prime count, and make prime count available as a
diff --git a/TODO b/TODO
index daade6c..f01fe4b 100644
--- a/TODO
+++ b/TODO
@@ -46,7 +46,5 @@
RiemannZeta
RiemannR
-- Add Lehmer in PP (I have a version, just needs some polishing). This is a
- high priority, as the test suite now assumes prime_count is fast for large
- inputs.
-
+- Dynamically use a mulmodadd in PP aks, just like the new C code does.
+ This will mean it'll work for full-size native ints.
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 0df6425..30041a2 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -254,10 +254,12 @@ sub primes {
if ($high > $_XS_MAXVAL) {
if ($_HAVE_GMP) {
$sref = Math::Prime::Util::GMP::primes($low,$high);
- # Convert the returned strings into BigInts
- croak "Internal error: large value without bigint loaded."
- unless defined $Math::BigInt::VERSION;
- @$sref = map { Math::BigInt->new("$_") } @$sref;
+ if ($high > ~0) {
+ # Convert the returned strings into BigInts
+ croak "Internal error: large value without bigint loaded."
+ unless defined $Math::BigInt::VERSION;
+ @$sref = map { Math::BigInt->new("$_") } @$sref;
+ }
return $sref;
}
return Math::Prime::Util::PP::primes($low,$high);
@@ -905,6 +907,9 @@ sub next_prime {
return _XS_next_prime($n) if $n <= $_XS_MAXVAL
&& (!defined $bigint::VERSION || $n < $_Config{'maxprime'} );
+ # Try to stick to the plan with respect to maximum return values.
+ return 0 if ref($_[0]) ne 'Math::BigInt' && $n >= $_Config{'maxprime'};
+
if ($_HAVE_GMP) {
# If $n is a bigint object, try to make the return value the same
return (ref($_[0]) eq 'Math::BigInt')
@@ -954,8 +959,12 @@ sub prime_count {
}
return _XS_prime_count($low,$high);
}
+ # We can relax these constraints if MPU::GMP gets a Lehmer implementation.
return Math::Prime::Util::GMP::prime_count($low,$high) if $_HAVE_GMP
- && defined &Math::Prime::Util::GMP::prime_count;
+ && defined &Math::Prime::Util::GMP::prime_count
+ && ( (ref($high) eq 'Math::BigInt')
+ || (($high-$low) < int($low/1_000_000))
+ );
return Math::Prime::Util::PP::prime_count($low,$high);
}
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index 5689490..73569c8 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -170,21 +170,17 @@ sub is_prime {
sub _sieve_erat_string {
my($end) = @_;
+ $end-- if ($end & 1) == 0;
+ my $s_end = $end >> 1;
- # Prefill with 3 and 5 already marked.
- my $whole = int( ($end>>1) / 15);
+ my $whole = int( $s_end / 15); # Prefill with 3 and 5 already marked.
croak "Sieve too large" if $whole > 1_145_324_612; # ~32 GB string
my $sieve = "100010010010110" . "011010010010110" x $whole;
- # Make exactly the number of entries requested, never more.
- substr($sieve, ($end>>1)+1) = '';
- my $n = 7;
- while ( ($n*$n) <= $end ) {
- my $s = $n*$n;
- my $filter_s = $s >> 1;
- my $filter_end = $end >> 1;
- while ($filter_s <= $filter_end) {
- substr($sieve, $filter_s, 1) = "1";
- $filter_s += $n;
+ substr($sieve, $s_end+1) = ''; # Ensure we don't make too many entries
+ my ($n, $limit) = ( 7, int(sqrt($end)) );
+ while ( $n <= $limit ) {
+ for (my $s = ($n*$n) >> 1; $s <= $s_end; $s += $n) {
+ substr($sieve, $s, 1) = '1';
}
do { $n += 2 } while substr($sieve, $n>>1, 1);
}
@@ -386,6 +382,96 @@ sub prev_prime {
#$d*30+$m;
}
+#############################################################################
+# Lehmer prime count
+#
+sub _mapes {
+ my($v, $ma) = @_;
+ return $v if $ma == 0;
+ my $val = $v-int($v/2);
+ $val += -int($v/3)+int($v/6) if $ma >= 2;
+ $val += -int($v/5)+int($v/10)+int($v/15)-int($v/30) if $ma >= 3;
+ $val += -int($v/7)+int($v/14)+int($v/21)-int($v/42)+int($v/35)-int($v/70)-int($v/105)+int($v/210) if $ma >= 4;
+ $val += -int($v/11)+int($v/22)+int($v/33)-int($v/66)+int($v/55)-int($v/110)-int($v/165)+int($v/330)+int($v/77)-int($v/154)-int($v/231)+int($v/462)-int($v/385)+int($v/770)+int($v/1155)-int($v/2310) if $ma >= 5;
+ $val += -int($v/13)+int($v/26)+int($v/39)-int($v/78)+int($v/65)-int($v/130)-int($v/195)+int($v/390)+int($v/91)-int($v/182)-int($v/273)+int($v/546)-int($v/455)+int($v/910)+int($v/1365)-int($v/2730)+int($v/143)-int($v/286)-int($v/429)+int($v/858)-int($v/715)+int($v/1430)+int($v/2145)-int($v/4290)-int($v/1001)+int($v/2002)+int($v/3003)-int($v/6006)+int($v/5005)-int($v/10010)-int($v/15015)+int($v/30030) if $ma >= 6;
+ return $val;
+}
+
+sub _legendre_phi {
+ my ($x, $a, $primes) = @_;
+ return _mapes($x,$a) if $a <= 6;
+ return ($x > 0 ? 1 : 0) if $x < $primes->[$a];
+
+ my $sum = 0;
+ my %vals = ( $x => 1 );
+ while ($a > 6) {
+ my $primea = $primes->[$a-1];
+ my %newvals;
+ while (my($v,$c) = each %vals) {
+ next if $c == 0;
+ # next if $v < $primea;
+ $newvals{$v} += $c;
+ my $sval = int($v / $primea);
+ if ($sval >= $primea) {
+ $newvals{$sval} -= $c;
+ } else {
+ $sum -= $c;
+ }
+ }
+ %vals = %newvals;
+ $a--;
+ }
+ while (my($v,$c) = each %vals) {
+ next if $c == 0;
+ $sum += $c * _mapes($v, $a);
+ }
+ return $sum;
+}
+
+sub _sieve_prime_count {
+ my $high = shift;
+ return (0,0,1,2,2,3,3)[$high] if $high < 7;
+ $high-- if ($high % 2) == 0; # Make high go to odd number.
+
+ my $sieveref = _sieve_erat($high);
+ my $count = 1 + $$sieveref =~ tr/0//;
+ return $count;
+}
+
+sub _lehmer_pi {
+ my $x = shift;
+ return _sieve_prime_count($x) if $x < 1_000;
+ my $z = int(sqrt($x+0.5));
+ my $a = _lehmer_pi(int(sqrt($z)+0.5));
+ my $b = _lehmer_pi($z);
+ my $c = _lehmer_pi(int($x**(1/3)+0.5));
+
+ # Generate at least b primes.
+ my $nth_prime_upper = ($b <= 10) ? 29 : int($b*(log($b) + log(log($b)))) + 1;
+ my $primes = primes( $nth_prime_upper );
+
+ my $sum = int(($b + $a - 2) * ($b - $a + 1) / 2);
+ $sum += _legendre_phi($x, $a, $primes);
+
+ # Make sure these calls are fast.
+ # 1M for 10**8, 32M for 10**10, 5600M for 10**12
+ prime_precalc( int($x / $primes->[$a]) );
+
+ foreach my $i ($a+1 .. $b) {
+ my $w = int($x / $primes->[$i-1]);
+ $sum = $sum - _sieve_prime_count($w);
+ if ($i <= $c) {
+ my $bi = _sieve_prime_count(int(sqrt($w)+0.5));
+ foreach my $j ($i .. $bi) {
+ $sum = $sum - _sieve_prime_count(int($w / $primes->[$j-1])) + $j - 1;
+ }
+ }
+ }
+ $sum;
+}
+#############################################################################
+
+
sub prime_count {
my($low,$high) = @_;
if (!defined $high) {
@@ -405,10 +491,9 @@ sub prime_count {
return $count if $low > $high;
if ( ref($low) eq 'Math::BigInt' || ref($high) eq 'Math::BigInt'
- || $high > 16_000_000_000
+ || ($high-$low) < 10
|| ($high-$low) < int($low/1_000_000) ) {
- # Too big to sieve.
- my $count = 0;
+ # Trial primes seems best.
my $curprime = next_prime($low-1);
while ($curprime <= $high) {
$count++;
@@ -417,16 +502,20 @@ sub prime_count {
return $count;
}
- my $sieveref;
- if ($low == 3) {
- $sieveref = _sieve_erat($high);
- } else {
- $sieveref = _sieve_segment($low,$high);
+ # TODO: Needs tuning
+ if ($high > 50_000) {
+ if ( ($high / ($high-$low+1)) < 100 ) {
+ $count += _lehmer_pi($high);
+ $count -= ($low == 3) ? 1 : _lehmer_pi($low-1);
+ return $count;
+ }
}
- $count += $$sieveref =~ tr/0//;
+ return (_sieve_prime_count($high) - 1 + $count) if $low == 3;
- $count;
+ my $sieveref = _sieve_segment($low,$high);
+ $count += $$sieveref =~ tr/0//;
+ return $count;
}
@@ -446,9 +535,21 @@ sub nth_prime {
my $prime = 0;
+ my $count = 1;
+ my $start = 3;
+
+ my $logn = log($n);
+ my $loglogn = log($logn);
+ my $nth_prime_upper = ($n <= 10) ? 29 : int($n*($logn + $loglogn)) + 1;
+ if ($nth_prime_upper > 100000) {
+ # Use fast Lehmer prime count combined with lower bound to get close.
+ my $nth_prime_lower = int($n * ($logn + $loglogn - 1.0 + (($loglogn-2.10)/$logn)));
+ $nth_prime_lower-- unless $nth_prime_lower % 2;
+ $count = _lehmer_pi($nth_prime_lower);
+ $start = $nth_prime_lower + 2;
+ }
+
{
- my $count = 1;
- my $start = 3;
# Make sure incr is an even number.
my $incr = ($n < 1000) ? 1000 : ($n < 10000) ? 10000 : 100000;
my $sieveref;
@@ -889,7 +990,18 @@ sub _test_anr {
sub is_aks_prime {
my $n = shift;
- $n = Math::BigInt->new("$n") unless ref($n) eq 'Math::BigInt';
+
+ if (ref($n) ne 'Math::BigInt') {
+ if (!defined $MATH::BigInt::VERSION) {
+ eval { require Math::BigInt; Math::BigInt->import(try=>'GMP,Pari'); 1; }
+ or do { croak "Cannot load Math::BigInt "; }
+ }
+ if (!defined $MATH::BigFloat::VERSION) {
+ eval { require Math::BigFloat; Math::BigFloat->import(); 1; }
+ or do { croak "Cannot load Math::BigFloat "; }
+ }
+ $n = Math::BigInt->new("$n");
+ }
return 0 if $n < 2;
return 0 if _is_perfect_power($n);
diff --git a/mulmod.h b/mulmod.h
index bce8eb7..92e8cb3 100644
--- a/mulmod.h
+++ b/mulmod.h
@@ -35,6 +35,13 @@
:"r"(a), "r"(b), "r"(c) /* input */
:"%rax", "%rdx" /* clobbered registers */
);
+ /* A version for _MSC_VER:
+ *
+ * __asm { mov rax, qword ptr a
+ * mul qword ptr b
+ * div qword ptr c
+ * mov qword ptr d, rdx }
+ */
return d;
}
#define mulmod(a,b,m) _mulmod(a,b,m)
diff --git a/t/14-nthprime.t b/t/14-nthprime.t
index 28abc61..55bca37 100644
--- a/t/14-nthprime.t
+++ b/t/14-nthprime.t
@@ -47,7 +47,7 @@ my %nthprimes64 = (
100000000000000000 => 4185296581467695669,
);
my %nthprimes_small = map { $_ => $nthprimes32{$_} }
- #grep { ($_ <= 2_000_000) || $extra }
+ grep { ($_ <= 10_000_000) || $extra }
keys %nthprimes32;
my @small_primes = (0, @{primes($nth_small_prime)});
--
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