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 <d...@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

_______________________________________________
Pkg-perl-cvs-commits mailing list
Pkg-perl-cvs-commits@lists.alioth.debian.org
http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-perl-cvs-commits

Reply via email to