[libmath-prime-util-perl] 03/55: improve twin prime approximations

[Partha P. Mukherjee]

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

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

commit 19ce02756b3d4f3277659a91861a0b2be2f194e3
Author: Dana Jacobsen <dana at acm.org>
Date:   Thu Apr 24 13:17:27 2014 -0700

    improve twin prime approximations
---
 Changes         |  5 +++++
 t/14-nthprime.t |  5 +++--
 util.c          | 20 ++++++++++----------
 3 files changed, 18 insertions(+), 12 deletions(-)

diff --git a/Changes b/Changes
index 09ee7da..5f8c641 100644
--- a/Changes
+++ b/Changes
@@ -1,5 +1,10 @@
 Revision history for Perl module Math::Prime::Util
 
+0.41  2014-05
+
+    - Small improvement to twin_prime_count_approx and nth_twin_prime_approx.
+
+
 0.40  2014-04-21
 
     [ADDED]
diff --git a/t/14-nthprime.t b/t/14-nthprime.t
index 03fdcab..aedaa4d 100644
--- a/t/14-nthprime.t
+++ b/t/14-nthprime.t
@@ -142,7 +142,8 @@ 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 $approx = nth_twin_prime_approx($n);
+  my $errorp = 100 * abs($nthtpc - $approx) / $nthtpc;
   my $estr = sprintf "%8.6f%%", $errorp;
-  cmp_ok( $errorp, '<=', 2, "nth_twin_prime_approx($n) is $estr");
+  cmp_ok( $errorp, '<=', 2, "nth_twin_prime_approx($n) is $estr (got $approx, expected ~$nthtpc)");
 }
diff --git a/util.c b/util.c
index 6ffba87..47225b8 100644
--- a/util.c
+++ b/util.c
@@ -590,13 +590,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 <    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)));
+    if      (n <    4000) li2 *= 1.2035 * logl(logl(logl(ln)));
+    else if (n <    8000) li2 *= 0.9439 * logl(logl(logl(ln*1000)));
+    else 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)));
+    else if (n < 1000000) li2 *= 0.8640 * logl(logl(logl(ln*16000)));
+    else if (n < 4000000) li2 *= 0.8627 * logl(logl(logl(ln*16000)));
+    else if (n <10000000) li2 *= 0.8536 * logl(logl(logl(ln*16000)));
     return (UV) (two_C2 * li2 + 0.5L);
   }
 }
@@ -985,18 +985,18 @@ UV nth_twin_prime_approx(UV n)
 {
   long double fn = (long double) n;
   long double flogn = logl(n);
+  long double fnlog2n = fn * flogn * flogn;
   UV lo, hi;
 
   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.
+   * These bounds are good for the actual nth_twin_prime values.
    */
-  lo = (UV) (0.0 * fn * flogn * flogn);
-  hi = (UV) (3.0 * fn * flogn * flogn + 3);
+  lo = (UV) (1.2 * fnlog2n);
+  hi = (UV) ( (n >= 1200) ? (1.7 * fnlog2n) : (2.3 * fnlog2n + 5) );
   if (hi <= lo) hi = UV_MAX;
   while (lo < hi) {
     UV mid = lo + (hi-lo)/2;

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