[libmath-prime-util-perl] 31/43: Add nth_twin_prime_approx, tighten twin_prime_count_approx

Partha P. Mukherjee ppm-guest at moszumanska.debian.org
Thu May 21 18:53:08 UTC 2015


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 b7b1849fa58c85db0216d6c1ca3462a230316191
Author: Dana Jacobsen <dana at acm.org>
Date:   Thu Apr 10 18:26:12 2014 -0700

    Add nth_twin_prime_approx, tighten twin_prime_count_approx
---
 Changes                |  1 +
 TODO                   |  2 ++
 XS.xs                  | 27 +++++++++++++++------------
 lib/Math/Prime/Util.pm |  3 ++-
 util.c                 | 44 ++++++++++++++++++++++++++++++++++++++------
 util.h                 |  1 +
 6 files changed, 59 insertions(+), 19 deletions(-)

diff --git a/Changes b/Changes
index 00611f2..c4caa69 100644
--- a/Changes
+++ b/Changes
@@ -9,6 +9,7 @@ Revision history for Perl module Math::Prime::Util
     - twin_prime_count                       counts twin primes in range
     - twin_prime_count_approx                fast approximation to Pi_2(n)
     - nth_twin_prime                         returns the nth twin prime
+    - nth_twin_prime_approx                  estimates the nth twin prime
 
     [FUNCTIONALITY AND PERFORMANCE]
 
diff --git a/TODO b/TODO
index 021d8b0..15e4cde 100644
--- a/TODO
+++ b/TODO
@@ -68,3 +68,5 @@
 - 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/XS.xs b/XS.xs
index 1cbabd2..6795342 100644
--- a/XS.xs
+++ b/XS.xs
@@ -604,10 +604,11 @@ next_prime(IN SV* svn)
     nth_prime_lower = 4
     nth_prime_approx = 5
     nth_twin_prime = 6
-    prime_count_upper = 7
-    prime_count_lower = 8
-    prime_count_approx = 9
-    twin_prime_count_approx = 10
+    nth_twin_prime_approx = 7
+    prime_count_upper = 8
+    prime_count_lower = 9
+    prime_count_approx = 10
+    twin_prime_count_approx = 11
   PPCODE:
     if (_validate_int(aTHX_ svn, 0)) {
       UV n = my_svuv(svn);
@@ -624,10 +625,11 @@ next_prime(IN SV* svn)
           case 4: ret = nth_prime_lower(n); break;
           case 5: ret = nth_prime_approx(n); break;
           case 6: ret = nth_twin_prime(n); break;
-          case 7: ret = prime_count_upper(n); break;
-          case 8: ret = prime_count_lower(n); break;
-          case 9 :ret = prime_count_approx(n); break;
-          case 10:
+          case 7: ret = nth_twin_prime_approx(n); break;
+          case 8: ret = prime_count_upper(n); break;
+          case 9: ret = prime_count_lower(n); break;
+          case 10:ret = prime_count_approx(n); break;
+          case 11:
           default:ret = twin_prime_count_approx(n); break;
         }
         XSRETURN_UV(ret);
@@ -645,10 +647,11 @@ next_prime(IN SV* svn)
       case 4:  _vcallsub_with_pp("nth_prime_lower");    break;
       case 5:  _vcallsub_with_pp("nth_prime_approx");   break;
       case 6:  _vcallsub_with_pp("nth_twin_prime");     break;
-      case 7:  _vcallsub_with_pp("prime_count_upper");  break;
-      case 8:  _vcallsub_with_pp("prime_count_lower");  break;
-      case 9:  _vcallsub_with_pp("prime_count_approx"); break;
-      case 10:
+      case 7:  _vcallsub_with_pp("nth_twin_prime_approx"); break;
+      case 8:  _vcallsub_with_pp("prime_count_upper");  break;
+      case 9:  _vcallsub_with_pp("prime_count_lower");  break;
+      case 10: _vcallsub_with_pp("prime_count_approx"); break;
+      case 11:
       default: _vcallsub_with_pp("twin_prime_count_approx"); break;
     }
     return; /* skip implicit PUTBACK */
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index d7a01af..e6026e9 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -33,7 +33,8 @@ our @EXPORT_OK =
       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 88a5efa..34bdf29 100644
--- a/util.c
+++ b/util.c
@@ -580,12 +580,21 @@ UV prime_count_approx(UV n)
 /* See http://numbers.computation.free.fr/Constants/Primes/twin.pdf, page 5 */
 UV twin_prime_count_approx(UV n)
 {
-  const long double two_C2 = 1.32032363169373914785562422L;
-  const long double two_over_log_two = 2.8853900817779268147198494L;
-  long double ln = (long double) n;
-  long double logn = logl(ln);
-  long double li2 = _XS_ExponentialIntegral(logn) + two_over_log_two - ln/logn;
-  return (UV) (two_C2 * li2 + 0.5L);
+  /* Best would be another estimate for n < ~ 5000 */
+  if (n < 2000) return twin_prime_count(3,n);
+  {
+    const long double two_C2 = 1.32032363169373914785562422L;
+    const long double two_over_log_two = 2.8853900817779268147198494L;
+    long double ln = (long double) 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)));
+    return (UV) (two_C2 * li2 + 0.5L);
+  }
 }
 
 UV prime_count_lower(UV n)
@@ -965,6 +974,29 @@ UV nth_twin_prime(UV n)
   return nth;
 }
 
+UV nth_twin_prime_approx(UV n)
+{
+  long double fn = (long double) n;
+  long double flogn = logl(n);
+  UV lo, hi;
+
+  if (n < 6)
+    return nth_twin_prime(n);
+
+  /* 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);
+  hi = (UV) (3.0 * fn * flogn * flogn + 3);
+  if (hi <= lo) hi = UV_MAX;
+  while (lo < hi) {
+    UV mid = lo + (hi-lo)/2;
+    if (twin_prime_count_approx(mid) < fn) lo = mid+1;
+    else                                   hi = mid;
+  }
+  return lo;
+}
+
 
 /* Return a char array with lo-hi+1 elements. mu[k-lo] = µ(k) for k = lo .. hi.
  * It is the callers responsibility to call Safefree on the result. */
diff --git a/util.h b/util.h
index 2130671..af8205c 100644
--- a/util.h
+++ b/util.h
@@ -23,6 +23,7 @@ extern UV  prime_count_approx(UV x);
 extern UV  twin_prime_count(UV low, UV high);
 extern UV  twin_prime_count_approx(UV n);
 extern UV  nth_twin_prime(UV n);
+extern UV  nth_twin_prime_approx(UV n);
 
 extern int powerof(UV n);
 extern int is_power(UV n, UV a);

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