[libmath-prime-util-perl] 47/50: Have nth_prime use Lehmer prime count on lower bound. 100x speedup for big numbers.

[Partha P. Mukherjee]

This is an automated email from the git hooks/post-receive script.

ppm-guest pushed a commit to annotated tag v0.13
in repository libmath-prime-util-perl.

commit 9e13087df19c0dd0b64ef94c3e9ed45cad6b6a8a
Author: Dana Jacobsen <dana at acm.org>
Date:   Sun Nov 18 02:42:52 2012 -0800

    Have nth_prime use Lehmer prime count on lower bound.  100x speedup for big numbers.
---
 Changes                    |  4 ++++
 TODO                       |  2 ++
 examples/bench-nthprime.pl |  4 ++--
 t/14-nthprime.t            |  2 +-
 util.c                     | 34 +++++++++++++++++++++++-----------
 5 files changed, 32 insertions(+), 14 deletions(-)

diff --git a/Changes b/Changes
index 3a29363..b6debc8 100644
--- a/Changes
+++ b/Changes
@@ -34,6 +34,10 @@ Revision history for Perl extension Math::Prime::Util.
     - prime_count will use Lehmer prime counting algorithm for largish
       sizes (above 4 million).  This is MUCH faster than sieving.
 
+    - nth_prime now uses the fast Lehmer prime count below the lower limit,
+      then sieves up from there.  This makes a big speed difference for inputs
+      over 10^6 or so -- over 100x faster for 10^9 and up.
+
 0.11  23 July 2012
     - Turn off threading tests on Cygwin, as threads on some Cygwin platforms
       give random panics (my Win7 64-bit works fine, XP 32-bit does not).
diff --git a/TODO b/TODO
index 7da582a..c704281 100644
--- a/TODO
+++ b/TODO
@@ -33,3 +33,5 @@
 - Add Lehmer in PP (I have a version, just needs some polishing).
 
 - Move AKS tests into their own test, and add some provable prime tests.
+
+- Add a 'verbose' option to general config.
diff --git a/examples/bench-nthprime.pl b/examples/bench-nthprime.pl
index 9d5742a..7656b95 100755
--- a/examples/bench-nthprime.pl
+++ b/examples/bench-nthprime.pl
@@ -11,10 +11,10 @@ my $count = shift || -5;
 
 srand(29);
 my @darray;
-push @darray, [gendigits($_,int(2700/($_*$_*$_)))]  for (2 .. 9);
+push @darray, [gendigits($_,int(5400/($_*$_*$_)))]  for (2 .. 10);
 
 my $sum;
-foreach my $digits (3 .. 9) {
+foreach my $digits (3 .. 10) {
   my @digarray = @{$darray[$digits-2]};
   my $numitems = scalar @digarray;
   my $timing = cmpthese(
diff --git a/t/14-nthprime.t b/t/14-nthprime.t
index 6979b9d..c4d3a9d 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 { ($_ <= 2000000) || $extra }
+                      #grep { ($_ <= 2000000) || $extra }
                       keys %nthprimes32;
 
 my @small_primes = (0, @{primes($nth_small_prime)});
diff --git a/util.c b/util.c
index d6aa254..979fe51 100644
--- a/util.c
+++ b/util.c
@@ -31,6 +31,7 @@ extern long double fabsl(long double);
 #include "sieve.h"
 #include "factor.h"
 #include "cache.h"
+#include "lehmer.h"
 
 static const unsigned char byte_zeros[256] =
   {8,7,7,6,7,6,6,5,7,6,6,5,6,5,5,4,7,6,6,5,6,5,5,4,6,5,5,4,5,4,4,3,
@@ -555,20 +556,31 @@ UV _XS_nth_prime(UV n)
   upper_limit = _XS_nth_prime_upper(n);
   MPUassert(upper_limit > 0, "nth_prime got an upper limit of 0");
 
-  /* Get the primary cache, and ensure it is at least this large.  If the
-   * input is small enough, get a sieve covering the range.  Otherwise, we'll
-   * walk segments.  Make sure we have enough primes in the cache so the
-   * segmented siever won't have to keep resieving.
+  /* For relatively small values, generate a sieve and count the results.
+   * For larger values, compute a lower bound, use Lehmer's algorithm to get
+   * a fast prime count, then start segment sieving from there.
    */
-  if (upper_limit <= (1*1024*1024*30))
+  if (upper_limit <= 1*1024*1024*30) {
+    /* Generate a sieve and count. */
     segment_size = get_prime_cache(upper_limit, &cache_sieve) / 30;
-  else
-    segment_size = get_prime_cache(sqrt(upper_limit), &cache_sieve) / 30;
+    /* Count up everything in the cached sieve. */
+    if (segment_size > 0)
+      count += count_segment_maxcount(cache_sieve, segment_size, target, &p);
+    release_prime_cache(cache_sieve);
+  } else {
+    double fn = n;
+    double flogn = log(fn);
+    double flog2n = log(flogn);   /* Dusart 2010, page 2, n >= 3 */
+    UV lower_limit = fn * (flogn + flog2n - 1.0 + ((flog2n-2.10)/flogn));
+    segment_size = lower_limit / 30;
+    lower_limit = 30 * segment_size - 1;
+    count = _XS_lehmer_pi(lower_limit) - 3;
+    MPUassert(count <= target, "Pi(nth_prime_lower(n))) > n");
+
+    /* Make sure the segment siever won't have to keep resieving. */
+    prime_precalc(sqrt(upper_limit));
+  }
 
-  /* Count up everything in the cached sieve. */
-  if (segment_size > 0)
-    count += count_segment_maxcount(cache_sieve, segment_size, target, &p);
-  release_prime_cache(cache_sieve);
   if (count == target)
     return p;
 

-- 
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-perl/packages/libmath-prime-util-perl.git