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

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

commit d3078911763322b75bbf7d254367b8c678f905b6
Author: Dana Jacobsen <d...@acm.org>
Date:   Fri Apr 11 15:43:21 2014 -0700

    Finish nth_twin_prime_approx
---
 TODO                        |  2 --
 lib/Math/Prime/Util.pm      | 10 +++++++++-
 lib/Math/Prime/Util/PP.pm   | 45 +++++++++++++++++++++++++++++++++++++++------
 lib/Math/Prime/Util/PPFE.pm |  5 +++++
 t/13-primecount.t           | 21 +++++++++++++++++++--
 t/14-nthprime.t             | 23 ++++++++++++++++++++++-
 t/92-release-pod-coverage.t |  3 ++-
 util.c                      | 15 ++++++++++-----
 8 files changed, 106 insertions(+), 18 deletions(-)

diff --git a/TODO b/TODO
index 15e4cde..021d8b0 100644
--- a/TODO
+++ b/TODO
@@ -68,5 +68,3 @@
 - Ensure a fast path for Math::GMP from MPU -> MPU:GMP -> GMP, and back.
 
 - znlog better implementation
-
-- Documentation, PP implementation, tests for nth_twin_prime_approx
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index e6026e9..480c333 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -1314,7 +1314,9 @@ for the answer, using some small table acceleration.
 
 Returns an approximation to the twin prime count of C<n>.  This returns
 quickly and has a very small error for large values.  The method used is
-conjecture B of Hardy and Littlewood 1922, as stated in Sebah and Gourdon 2002.
+conjecture B of Hardy and Littlewood 1922, as stated in
+Sebah and Gourdon 2002.  For inputs under 10M, a correction factor is 
+additionally applied to reduce the mean squared error.
 
 
 =head2 nth_prime
@@ -1377,6 +1379,12 @@ correction.
 Returns the Nth twin prime.  This is done via sieving and counting, so
 is not very fast for large values.
 
+=head2 nth_twin_prime_approx
+
+Returns an approximation to the Nth twin prime.  A curve fit is used for
+small inputs (under 1200), while for larger inputs a binary search is done
+on the approximate twin prime count.
+
 
 =head2 is_pseudoprime
 
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index bc3f0bb..be6adee 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -149,12 +149,15 @@ sub _validate_positive_integer {
 }
 
 
-my @_primes_small = (
-   0,2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,
-   101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,
-   193,197,199,211,223,227,229,233,239,241,251,257,263,269,271,277,281,283,
-   293,307,311,313,317,331,337,347,349,353,359,367,373,379,383,389,397,401,
-   409,419,421,431,433,439,443,449,457,461,463,467,479,487,491,499,503,509);
+my @_primes_small = (0,2,3,5);
+{
+  my $sieveref = _sieve_erat_string(5000);
+  my $n = 5;
+  foreach my $s (split("0", substr($$sieveref, 3), -1)) {
+    $n += 2 + 2 * length($s);
+    push @_primes_small, $n if $n <= 5000;
+  }
+}
 my @_prime_next_small = (
    2,2,3,5,5,7,7,11,11,11,11,13,13,17,17,17,17,19,19,23,23,23,23,
    29,29,29,29,29,29,31,31,37,37,37,37,37,37,41,41,41,41,43,43,47,
@@ -459,6 +462,11 @@ sub next_prime {
   _validate_positive_integer($n);
 
   return $_prime_next_small[$n] if $n <= $#_prime_next_small;
+  if ($n < $_primes_small[-1]) {
+    my $i = _tiny_prime_count($n);  # Binary search
+    return $_primes_small[$i+1];
+  }
+
 
   $n = Math::BigInt->new(''.$_[0])
      if ref($n) ne 'Math::BigInt' && $n >= MPU_MAXPRIME;
@@ -1353,9 +1361,17 @@ sub twin_prime_count {
 
 sub twin_prime_count_approx {
   my($n) = @_;
+  return twin_prime_count(3,$n) if $n < 2000;
   $n = _upgrade_to_float($n) if ref($n);
   my $logn = log($n);
   my $li2 = ExponentialIntegral($logn) + 2.8853900817779268147198494 - 
($n/$logn);
+  if    ($n <     4000) { $li2 *= 1.0005 * log(log(log($n*4000))); }
+  elsif ($n <     8000) { $li2 *= 0.9734 * log(log(log($n*4000))); }
+  elsif ($n <    32000) { $li2 *= 0.8967 * log(log(log($n*4000))); }
+  elsif ($n <   200000) { $li2 *= 0.8937 * log(log(log($n*4000))); }
+  elsif ($n <  1000000) { $li2 *= 0.8793 * log(log(log($n*4000))); }
+  elsif ($n <  4000000) { $li2 *= 0.8766 * log(log(log($n*4000))); }
+  elsif ($n < 10000000) { $li2 *= 0.8664 * log(log(log($n*4000))); }
   return int(1.32032363169373914785562422 * $li2 + 0.5);
 }
 
@@ -1369,6 +1385,23 @@ sub nth_twin_prime {
   $nth;
 }
 
+sub nth_twin_prime_approx {
+  my($n) = @_;
+  return nth_twin_prime($n) if $n < 6;
+  $n = _upgrade_to_float($n) if ref($n);
+  my $nlogn2 = $n * log($n) * log($n);
+  return int(5.023 * $nlogn2/log(log(1600*$n*$n))) if $n > 59 && $n < 1200;
+
+  my $lo = int(1.0 * $nlogn2);
+  my $hi = int(3.0 * $nlogn2 + 3);
+  while ($lo < $hi) {
+    my $mid = $lo + (($hi-$lo) >> 1);
+    if (twin_prime_count_approx($mid) < $n) { $lo = $mid+1; }
+    else                                    { $hi = $mid;   }
+  }
+  return $lo;
+}
+
 
 #############################################################################
 
diff --git a/lib/Math/Prime/Util/PPFE.pm b/lib/Math/Prime/Util/PPFE.pm
index 77d917d..ac5c509 100644
--- a/lib/Math/Prime/Util/PPFE.pm
+++ b/lib/Math/Prime/Util/PPFE.pm
@@ -126,6 +126,11 @@ sub nth_twin_prime {
   _validate_positive_integer($n);
   return Math::Prime::Util::PP::nth_twin_prime($n);
 }
+sub nth_twin_prime_approx {
+  my($n) = @_;
+  _validate_positive_integer($n);
+  return Math::Prime::Util::PP::nth_twin_prime_approx($n);
+}
 
 
 sub is_prime {
diff --git a/t/13-primecount.t b/t/13-primecount.t
index 545e8b6..4df4e04 100644
--- a/t/13-primecount.t
+++ b/t/13-primecount.t
@@ -4,7 +4,8 @@ use warnings;
 
 use Test::More;
 use Math::Prime::Util qw/prime_count twin_prime_count
-                        prime_count_lower prime_count_upper 
prime_count_approx/;
+                         prime_count_lower prime_count_upper
+                         prime_count_approx twin_prime_count_approx/;
 
 my $isxs  = Math::Prime::Util::prime_get_config->{'xs'};
 my $use64 = Math::Prime::Util::prime_get_config->{'maxbits'} > 32;
@@ -82,6 +83,16 @@ my %intervals = (
 );
 delete @intervals{ grep { (parse_range($_))[1] > ~0 } keys %intervals };
 
+my %tpcs = (
+              5000 =>           126,
+            500000 =>          4565,
+          50000000 =>        239101,
+        5000000000 =>      14618166,
+      500000000000 =>     986222314,
+    50000000000000 =>   71018282471,
+  5000000000000000 => 5357875276068,
+);
+
 plan tests => 0 + 1
                 + 3*scalar(keys %pivals32)
                 + scalar(keys %pivals_small)
@@ -89,7 +100,7 @@ plan tests => 0 + 1
                 + scalar(keys %intervals)
                 + 1
                 + 5 + 2*$extra # prime count specific methods
-                + 3 + (($isxs && $use64) ? 1 : 0); # twin prime counts
+                + 3 + (($isxs && $use64) ? 1+2*scalar(keys %tpcs) : 0);# twin 
pc
 
 ok( eval { prime_count(13); 1; }, "prime_count in void context");
 
@@ -182,4 +193,10 @@ is(twin_prime_count(10**8,10**8+34587), 137, "twin prime 
count 10^8 to +34587");
 is(twin_prime_count(654321), 5744, "twin prime count 654321");
 if ($isxs && $use64) {
   is(twin_prime_count(1000000000123456), 1177209242446, "twin prime count 
1000000000123456");
+  while (my($n, $tpc) = each (%tpcs)) {
+    is(twin_prime_count($n), $tpc, "twin prime count $n");
+    my $errorp = 100 * abs($tpc - twin_prime_count_approx($n)) / $tpc;
+    my $estr = sprintf "%8.6f%%", $errorp;
+    cmp_ok( $errorp, '<=', 2, "twin_prime_count_approx($n) is $estr");
+  }
 }
diff --git a/t/14-nthprime.t b/t/14-nthprime.t
index c026157..03fdcab 100644
--- a/t/14-nthprime.t
+++ b/t/14-nthprime.t
@@ -4,7 +4,8 @@ use warnings;
 
 use Test::More;
 use Math::Prime::Util qw/primes nth_prime nth_twin_prime
-                         nth_prime_lower nth_prime_upper nth_prime_approx/;
+                         nth_prime_lower nth_prime_upper
+                         nth_prime_approx nth_twin_prime_approx/;
 
 my $use64 = Math::Prime::Util::prime_get_config->{'maxbits'} > 32;
 my $usexs = Math::Prime::Util::prime_get_config->{'xs'};
@@ -62,6 +63,19 @@ my %nthprimes_small = map { $_ => $nthprimes32{$_} }
 
 my @small_primes = (0, @{primes($nth_small_prime)});
 
+my %ntpcs = (
+             5 =>                   29,
+            50 =>                 1487,
+           500 =>                32411,
+          5000 =>               557519,
+         50000 =>              8264957,
+        500000 =>            115438667,
+       5000000 =>           1523975909,
+      50000000 =>          19358093939,
+     500000000 =>         239211160649,
+);
+
+
 plan tests => 0 + 2*scalar(keys %pivals32)
                 + 1
                 + 3*scalar(keys %nthprimes32)
@@ -69,6 +83,7 @@ plan tests => 0 + 2*scalar(keys %pivals32)
                 + $use64 * 3 * scalar(keys %nthprimes64)
                 + 3   # nth_prime_lower with max index
                 + 3   # nth_twin_prime
+                + scalar(keys %ntpcs)   # nth_twin_prime_approx
                 + (($extra && $use64 && $usexs) ? 1 : 0);
 
 
@@ -125,3 +140,9 @@ if ($extra && $use64 && $usexs) {
 is( nth_twin_prime(0), 0, "nth_twin_prime(0) = 0" );
 is( nth_twin_prime(17), 239, "239 = 17th twin prime" );
 is( nth_twin_prime(1234), 101207, "101207 = 1234'th twin prime" );
+
+while (my($n, $nthtpc) = each (%ntpcs)) {
+  my $errorp = 100 * abs($nthtpc - nth_twin_prime_approx($n)) / $nthtpc;
+  my $estr = sprintf "%8.6f%%", $errorp;
+  cmp_ok( $errorp, '<=', 2, "nth_twin_prime_approx($n) is $estr");
+}
diff --git a/t/92-release-pod-coverage.t b/t/92-release-pod-coverage.t
index 9b83b6e..b63ddd6 100644
--- a/t/92-release-pod-coverage.t
+++ b/t/92-release-pod-coverage.t
@@ -58,7 +58,8 @@ sub mpu_public_regex {
       prime_count
       prime_count_lower prime_count_upper prime_count_approx
       nth_prime nth_prime_lower nth_prime_upper nth_prime_approx
-      twin_prime_count twin_prime_count_approx nth_twin_prime
+      twin_prime_count twin_prime_count_approx
+      nth_twin_prime nth_twin_prime_approx
       random_prime random_ndigit_prime random_nbit_prime random_strong_prime
       random_proven_prime random_proven_prime_with_cert
       random_maurer_prime random_maurer_prime_with_cert
diff --git a/util.c b/util.c
index 34bdf29..a33bcbf 100644
--- a/util.c
+++ b/util.c
@@ -589,10 +589,13 @@ UV twin_prime_count_approx(UV n)
     long double logn = logl(ln);
     long double li2 = _XS_ExponentialIntegral(logn) + two_over_log_two-ln/logn;
     /* try to minimize MSE */
-    if      (n <    7000) li2 *= 0.8991 * logl(logl(logl(ln*4000)));
-    else if (n <  200000) li2 *= 0.8938 * logl(logl(logl(ln*4000)));
-    else if (n < 1000000) li2 *= 0.8794 * logl(logl(logl(ln*4000)));
-    else if (n < 2000000) li2 *= 0.8808 * logl(logl(logl(ln*4000)));
+    if      (n <    4000) li2 *= 1.0005 * logl(logl(logl(ln*4000)));
+    if      (n <    8000) li2 *= 0.9734 * logl(logl(logl(ln*4000)));
+    if      (n <   32000) li2 *= 0.8967 * logl(logl(logl(ln*4000)));
+    else if (n <  200000) li2 *= 0.8937 * logl(logl(logl(ln*4000)));
+    else if (n < 1000000) li2 *= 0.8793 * logl(logl(logl(ln*4000)));
+    else if (n < 4000000) li2 *= 0.8766 * logl(logl(logl(ln*4000)));
+    else if (n <10000000) li2 *= 0.8664 * logl(logl(logl(ln*4000)));
     return (UV) (two_C2 * li2 + 0.5L);
   }
 }
@@ -982,11 +985,13 @@ UV nth_twin_prime_approx(UV n)
 
   if (n < 6)
     return nth_twin_prime(n);
+  if (n > 59 && n < 1200)  /* Curve fit */
+    return (UV) (5.023 * fn*flogn*flogn / logl(logl(1600*fn*fn)));
 
   /* Binary search on the TPC estimate.
    * Good results require that the TPC estimate is both fast and accurate.
    */
-  lo = (UV) (1.0 * fn * flogn * flogn);
+  lo = (UV) (0.0 * fn * flogn * flogn);
   hi = (UV) (3.0 * fn * flogn * flogn + 3);
   if (hi <= lo) hi = UV_MAX;
   while (lo < hi) {

-- 
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

Reply via email to