[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