[libmath-prime-util-perl] 14/43: Add twin_prime_count_approx

[Partha P. Mukherjee]

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 edc0a9edb932351ae71b360f57d970889493c100
Author: Dana Jacobsen <dana at acm.org>
Date:   Thu Mar 20 18:42:07 2014 -0700

    Add twin_prime_count_approx
---
 Changes                |  1 +
 XS.xs                  | 11 +++++++----
 lib/Math/Prime/Util.pm |  9 ++++++++-
 util.c                 | 11 +++++++++++
 util.h                 |  1 +
 5 files changed, 28 insertions(+), 5 deletions(-)

diff --git a/Changes b/Changes
index 3cba617..af0b19a 100644
--- a/Changes
+++ b/Changes
@@ -7,6 +7,7 @@ Revision history for Perl module Math::Prime::Util
     - random_shawe_taylor_prime              FIPS 186-4 random proven prime
     - random_shawe_taylor_prime_with_cert    as above with certificate
     - twin_prime_count                       counts twin primes in range
+    - twin_prime_count_approx                fast approximation to Pi_2(n)
 
     [FUNCTIONALITY AND PERFORMANCE]
 
diff --git a/XS.xs b/XS.xs
index 0e6da43..cd9bfb8 100644
--- a/XS.xs
+++ b/XS.xs
@@ -606,6 +606,7 @@ next_prime(IN SV* svn)
     prime_count_upper = 6
     prime_count_lower = 7
     prime_count_approx = 8
+    twin_prime_count_approx = 9
   PPCODE:
     if (_validate_int(aTHX_ svn, 0)) {
       UV n = my_svuv(svn);
@@ -623,8 +624,9 @@ next_prime(IN SV* svn)
           case 5: ret = nth_prime_approx(n); break;
           case 6: ret = prime_count_upper(n); break;
           case 7: ret = prime_count_lower(n); break;
-          case 8:
-          default:ret = prime_count_approx(n); break;
+          case 8 :ret = prime_count_approx(n); break;
+          case 9:
+          default:ret = twin_prime_count_approx(n); break;
         }
         XSRETURN_UV(ret);
       }
@@ -642,8 +644,9 @@ next_prime(IN SV* svn)
       case 5:  _vcallsub_with_pp("nth_prime_approx");   break;
       case 6:  _vcallsub_with_pp("prime_count_upper");  break;
       case 7:  _vcallsub_with_pp("prime_count_lower");  break;
-      case 8:
-      default: _vcallsub_with_pp("prime_count_approx"); break;
+      case 8:  _vcallsub_with_pp("prime_count_approx"); break;
+      case 9:
+      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 bee3b52..f3785b6 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -30,9 +30,10 @@ our @EXPORT_OK =
       forprimes forcomposites fordivisors
       prime_iterator prime_iterator_object
       next_prime  prev_prime
-      prime_count twin_prime_count
+      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
       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
@@ -1307,6 +1308,12 @@ range).
 There is no useful formula known for this, unlike prime counts.  We sieve
 for the answer, using some small table acceleration.
 
+=head2 twin_prime_count_approx
+
+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.
+
 
 =head2 nth_prime
 
diff --git a/util.c b/util.c
index cda7e5a..8249387 100644
--- a/util.c
+++ b/util.c
@@ -577,6 +577,17 @@ UV prime_count_approx(UV n)
   return (UV) (_XS_RiemannR( (long double) n ) + 0.5 );
 }
 
+/* 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);
+}
+
 UV prime_count_lower(UV n)
 {
   long double fn, flogn, lower, a;
diff --git a/util.h b/util.h
index c88b8f3..e730a4b 100644
--- a/util.h
+++ b/util.h
@@ -21,6 +21,7 @@ extern UV  prime_count_lower(UV x);
 extern UV  prime_count_upper(UV x);
 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 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