[libmath-prime-util-perl] 13/43: Add twin_prime_count
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:53:06 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 f12cbbbb60472b86695a64b1d6d77a6960c61177
Author: Dana Jacobsen <dana at acm.org>
Date: Thu Mar 20 18:08:20 2014 -0700
Add twin_prime_count
---
Changes | 1 +
XS.xs | 12 ++++++++--
lib/Math/Prime/Util.pm | 20 +++++++++++++++-
lib/Math/Prime/Util/PP.pm | 15 ++++++++++++
lib/Math/Prime/Util/PPFE.pm | 11 +++++++++
t/13-primecount.t | 11 +++++++--
util.c | 56 +++++++++++++++++++++++++++++++++++++++++++++
util.h | 1 +
8 files changed, 122 insertions(+), 5 deletions(-)
diff --git a/Changes b/Changes
index 7c79f09..3cba617 100644
--- a/Changes
+++ b/Changes
@@ -6,6 +6,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
[FUNCTIONALITY AND PERFORMANCE]
diff --git a/XS.xs b/XS.xs
index 1523068..0e6da43 100644
--- a/XS.xs
+++ b/XS.xs
@@ -305,6 +305,7 @@ void
prime_count(IN SV* svlo, ...)
ALIAS:
_XS_segment_pi = 1
+ twin_prime_count = 2
PREINIT:
int lostatus, histatus;
UV lo, hi;
@@ -321,7 +322,9 @@ prime_count(IN SV* svlo, ...)
hi = my_svuv(ST(1));
}
if (lo <= hi) {
- if (ix == 1 || (hi / (hi-lo+1)) > 100) {
+ if (ix == 2) {
+ count = twin_prime_count(lo, hi);
+ } else if (ix == 1 || (hi / (hi-lo+1)) > 100) {
count = _XS_prime_count(lo, hi);
} else {
count = _XS_LMO_pi(hi);
@@ -331,7 +334,12 @@ prime_count(IN SV* svlo, ...)
}
XSRETURN_UV(count);
}
- _vcallsubn(aTHX_ GIMME_V, VCALL_ROOT, "_generic_prime_count", items);
+ switch (ix) {
+ case 0:
+ case 1: _vcallsubn(aTHX_ GIMME_V, VCALL_ROOT, "_generic_prime_count", items); break;
+ case 2:
+ default:_vcallsub_with_pp("twin_prime_count"); break;
+ }
return; /* skip implicit PUTBACK */
UV
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 651eaaf..bee3b52 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -30,7 +30,7 @@ our @EXPORT_OK =
forprimes forcomposites fordivisors
prime_iterator prime_iterator_object
next_prime prev_prime
- prime_count
+ prime_count twin_prime_count
prime_count_lower prime_count_upper prime_count_approx
nth_prime nth_prime_lower nth_prime_upper nth_prime_approx
random_prime random_ndigit_prime random_nbit_prime random_strong_prime
@@ -1290,6 +1290,24 @@ A slightly faster but much less accurate answer can be obtained by averaging
the upper and lower bounds.
+=head2 twin_prime_count
+
+Similar to prime count, but returns the count of twin primes (primes C<p>
+where C<p+2> is also prime). Takes either a single number indicating a count
+from 2 to the argument, or two numbers indicating a range.
+
+The primes being counted are the first value, so a range of C<(3,5)> will
+return a count of two, because both C<3> and C<5> are counted as twin primes.
+A range of C<(12,13)> will return a count of zero, because neither C<12+2>
+nor C<13+2> are prime. In contrast, C<primesieve> requires all elements of
+a constellation to be within the range to be counted, so would return one for
+the first example (C<5> is not counted because its pair C<7> is not in the
+range).
+
+There is no useful formula known for this, unlike prime counts. We sieve
+for the answer, using some small table acceleration.
+
+
=head2 nth_prime
say "The ten thousandth prime is ", nth_prime(10_000);
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index 9855eac..d2bc825 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -1334,6 +1334,21 @@ sub prime_count_upper {
return int($result);
}
+sub twin_prime_count {
+ my($low,$high) = @_;
+ my $sum = 0;
+ # TODO: I suspect calling primes() on segments would be faster in most cases.
+ if ($high >= $low) {
+ my $p = next_prime($low-1);
+ my $p2 = next_prime($p);
+ while ($p <= $high) {
+ $sum++ if ($p2-$p) == 2;
+ ($p, $p2) = ($p2, next_prime($p2));
+ }
+ }
+ $sum;
+}
+
#############################################################################
diff --git a/lib/Math/Prime/Util/PPFE.pm b/lib/Math/Prime/Util/PPFE.pm
index 70a1464..bda8583 100644
--- a/lib/Math/Prime/Util/PPFE.pm
+++ b/lib/Math/Prime/Util/PPFE.pm
@@ -105,6 +105,17 @@ sub prime_count_approx {
_validate_positive_integer($n);
return Math::Prime::Util::PP::prime_count_approx($n);
}
+sub twin_prime_count {
+ my($low,$high) = @_;
+ if (scalar @_ > 1) {
+ _validate_positive_integer($low);
+ _validate_positive_integer($high);
+ } else {
+ ($low,$high) = (2, $low);
+ _validate_positive_integer($high);
+ }
+ return Math::Prime::Util::PP::twin_prime_count($low,$high);
+}
sub is_prime {
diff --git a/t/13-primecount.t b/t/13-primecount.t
index 2ba8afb..2d87671 100644
--- a/t/13-primecount.t
+++ b/t/13-primecount.t
@@ -3,7 +3,8 @@ use strict;
use warnings;
use Test::More;
-use Math::Prime::Util qw/prime_count prime_count_lower prime_count_upper prime_count_approx/;
+use Math::Prime::Util qw/prime_count twin_prime_count
+ prime_count_lower prime_count_upper prime_count_approx/;
my $isxs = Math::Prime::Util::prime_get_config->{'xs'};
my $use64 = Math::Prime::Util::prime_get_config->{'maxbits'} > 32;
@@ -87,7 +88,8 @@ plan tests => 0 + 1
+ $use64 * 3 * scalar(keys %pivals64)
+ scalar(keys %intervals)
+ 1
- + 5 + 2*$extra; # prime count specific methods
+ + 5 + 2*$extra # prime count specific methods
+ + 3; # twin prime counts
ok( eval { prime_count(13); 1; }, "prime_count in void context");
@@ -173,3 +175,8 @@ if ($extra) {
is(Math::Prime::Util::PP::_lehmer_pi (3456789), 247352, "PP Lehmer count");
is(Math::Prime::Util::PP::_sieve_prime_count(3456789), 247352, "PP sieve count");
}
+
+####### Twin prime counts
+is(twin_prime_count(13,31), 2, "twin prime count 13 to 31");
+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");
diff --git a/util.c b/util.c
index be03129..cda7e5a 100644
--- a/util.c
+++ b/util.c
@@ -848,6 +848,62 @@ UV nth_prime(UV n)
return ( (segbase*30) + p );
}
+UV twin_prime_count(UV beg, UV end)
+{
+ unsigned char* segment;
+ UV sum = 0;
+
+ if (beg <= 2) {
+#if BITS_PER_WORD == 64
+ if(end>=UVCONST(10000000000000)){sum=15834664872;beg=10000000000000;}
+ else if(end>=UVCONST( 1000000000000)){sum= 1870585220;beg= 1000000000000;}
+ else if(end>=UVCONST( 100000000000)){sum= 224376048;beg= 100000000000;}
+ else if(end>=UVCONST( 10000000000)){sum= 27412679;beg= 10000000000;}
+ else if(end>=UVCONST( 1000000000)){sum= 3424506;beg= 1000000000;}
+#else
+ if(end>=UVCONST( 1000000000)){sum= 3424506;beg= 1000000000;}
+#endif
+ }
+
+ if (beg <= 3 && end >= 3) sum++;
+ if (beg <= 5 && end >= 5) sum++;
+ if (beg < 11) beg = 7;
+ if (beg <= end) {
+ /* Make end points odd */
+ beg |= 1;
+ end = (end-1) | 1;
+ /* Cheesy way of counting the partial-byte edges */
+ while ((beg % 30) != 1) {
+ if (_XS_is_prime(beg) && _XS_is_prime(beg+2) && beg <= end) sum++;
+ beg += 2;
+ }
+ while ((end % 30) != 29) {
+ if (_XS_is_prime(end) && _XS_is_prime(end+2) && beg <= end) sum++;
+ end -= 2; if (beg > end) break;
+ }
+ }
+ if (beg <= end) {
+ UV seg_base, seg_low, seg_high;
+ void* ctx = start_segment_primes(beg, end, &segment);
+ while (next_segment_primes(ctx, &seg_base, &seg_low, &seg_high)) {
+ UV p, bytes = (seg_high-seg_low+29)/30;
+ for (p = 0; p < bytes; p++) {
+ UV s = segment[p];
+ if (!(s & 0x0C)) sum++;
+ if (!(s & 0x30)) sum++;
+ if (!(s & 0x80)) {
+ if (p+1 < bytes) { if (!(segment[p+1] & 0x01)) sum++; }
+ else { if (_XS_is_prime(seg_high+2)) sum++; }
+ }
+ }
+ }
+ end_segment_primes(ctx);
+ }
+ return sum;
+}
+
+
+
/* 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. */
#define PGTLO(p,lo) ((p) >= lo) ? (p) : ((p)*(lo/(p)) + ((lo%(p))?(p):0))
diff --git a/util.h b/util.h
index 978e192..c88b8f3 100644
--- a/util.h
+++ b/util.h
@@ -20,6 +20,7 @@ extern UV nth_prime_approx(UV x);
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 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