[libmath-prime-util-perl] 02/05: Implementation
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:44:12 UTC 2015
This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.01
in repository libmath-prime-util-perl.
commit db6e2788b50cb7091802a65eb95875f411aef377
Author: Dana Jacobsen <dana at acm.org>
Date: Sun Jun 3 17:59:45 2012 -0600
Implementation
---
Changes | 4 +
MANIFEST | 26 +++
Makefile.PL | 48 +++++
TODO | 6 +
XS.xs | 216 +++++++++++++++++++
bitarray.h | 54 +++++
lib/Math/Prime/Util.pm | 353 ++++++++++++++++++++++++++++++
sieve.c | 250 ++++++++++++++++++++++
sieve.h | 85 ++++++++
t/01-load.t | 7 +
t/02-can.t | 16 ++
t/03-init.t | 22 ++
t/10-isprime.t | 112 ++++++++++
t/11-primes.t | 96 +++++++++
t/12-nextprime.t | 64 ++++++
t/13-primecount.t | 52 +++++
t/14-nthprime.t | 71 +++++++
t/30-relations.t | 36 ++++
t/90-release-perlcritic.t | 24 +++
t/91-release-pod-syntax.t | 19 ++
t/92-release-pod-coverage.t | 25 +++
util.c | 507 ++++++++++++++++++++++++++++++++++++++++++++
util.h | 22 ++
23 files changed, 2115 insertions(+)
diff --git a/Changes b/Changes
new file mode 100644
index 0000000..0660a2b
--- /dev/null
+++ b/Changes
@@ -0,0 +1,4 @@
+Revision history for Perl extension Math::Prime::Util.
+
+0.01 3 June 2012
+ - Original release
diff --git a/MANIFEST b/MANIFEST
new file mode 100644
index 0000000..a0a6eb4
--- /dev/null
+++ b/MANIFEST
@@ -0,0 +1,26 @@
+Changes
+lib/Math/Prime/Util.pm
+LICENSE
+Makefile.PL
+MANIFEST
+README
+TODO
+XS.xs
+bitarray.h
+sieve.h
+sieve.c
+util.h
+util.c
+t/01-load.t
+t/02-can.t
+t/03-init.t
+t/10-isprime.t
+t/11-primes.t
+t/12-nextprime.t
+t/13-primecount.t
+t/14-nthprime.t
+t/30-relations.t
+t/90-release-perlcritic.t
+t/91-release-pod-syntax.t
+t/92-release-pod-coverage.t
+
diff --git a/Makefile.PL b/Makefile.PL
new file mode 100644
index 0000000..afb943b
--- /dev/null
+++ b/Makefile.PL
@@ -0,0 +1,48 @@
+use ExtUtils::MakeMaker;
+
+WriteMakefile1(
+ NAME => 'Math::Prime::Util',
+ ABSTRACT => 'Utilities related to prime numbers, including fast generators / sievers',
+ VERSION_FROM => 'lib/Math/Prime/Util.pm',
+ LICENSE => 'perl',
+ AUTHOR => 'Dana A Jacobsen <dana at acm.org>',
+
+ LIBS => [''], # e.g., '-lm'
+ DEFINE => '', # e.g., '-DHAVE_SOMETHING'
+ INC => '', # e.g., '-I/usr/include/other'
+ OBJECT => 'sieve.o util.o XS.o',
+
+ PREREQ_PM => {
+ 'Test::More' => '0.45',
+ 'Exporter' => '5.562',
+ 'XSLoader' => '0.01',
+ },
+ # Example:
+ #META_MERGE => {
+ # recommends => { 'Data::BitStream' => 0, },
+ # },
+
+ MIN_PERL_VERSION => 5.006002,
+);
+
+sub WriteMakefile1 {
+ my %params = @_;
+ my $eumm_version = $ExtUtils::MakeMaker::VERSION;
+ $eumm_version = eval $eumm_version;
+
+ if ($params{BUILD_REQUIRES} and $eumm_version < 6.5503) {
+ #EUMM 6.5502 has problems with BUILD_REQUIRES
+ $params{PREREQ_PM}={ %{$params{PREREQ_PM} || {}} , %{$params{BUILD_REQUIRES}} };
+ delete $params{BUILD_REQUIRES};
+ }
+ delete $params{CONFIGURE_REQUIRES} if $eumm_version < 6.52;
+ delete $params{MIN_PERL_VERSION} if $eumm_version < 6.48;
+ delete $params{META_MERGE} if $eumm_version < 6.46;
+ delete $params{META_ADD} if $eumm_version < 6.46;
+ delete $params{LICENSE} if $eumm_version < 6.31;
+ delete $params{AUTHOR} if $] < 5.005;
+ delete $params{ABSTRACT_FROM} if $] < 5.005;
+ delete $params{BINARY_LOCATION} if $] < 5.005;
+
+ WriteMakefile(%params);
+}
diff --git a/TODO b/TODO
new file mode 100644
index 0000000..4ddd74e
--- /dev/null
+++ b/TODO
@@ -0,0 +1,6 @@
+
+- Examine behavior near 32-bit limit on 32-bit machines.
+
+- GMP versions of all routines.
+
+- factoring (Fermat, Pollard Rho, HOLF, etc.)
diff --git a/XS.xs b/XS.xs
new file mode 100644
index 0000000..a2ecd9b
--- /dev/null
+++ b/XS.xs
@@ -0,0 +1,216 @@
+
+#include "EXTERN.h"
+#include "perl.h"
+#include "XSUB.h"
+/* We're not using anything for which we need ppport.h */
+#include "sieve.h"
+#include "util.h"
+#include "bitarray.h"
+
+MODULE = Math::Prime::Util PACKAGE = Math::Prime::Util
+
+PROTOTYPES: ENABLE
+
+
+void
+prime_precalc(IN UV n)
+
+void
+prime_free()
+
+UV
+prime_count(IN UV n)
+ CODE:
+ if (GIMME_V == G_VOID) {
+ prime_precalc(n);
+ RETVAL = 0;
+ } else {
+ RETVAL = prime_count(n);
+ }
+ OUTPUT:
+ RETVAL
+
+UV
+prime_count_lower(IN UV n)
+
+UV
+prime_count_upper(IN UV n)
+
+UV
+prime_count_approx(IN UV n)
+
+UV
+nth_prime(IN UV n)
+
+UV
+nth_prime_lower(IN UV n)
+
+UV
+nth_prime_upper(IN UV n)
+
+UV
+nth_prime_approx(IN UV n)
+
+int
+is_prime(IN UV n)
+
+UV
+next_prime(IN UV n)
+
+UV
+prev_prime(IN UV n)
+
+
+UV
+_get_prime_cache_size()
+ CODE:
+ RETVAL = get_prime_cache_size();
+ OUTPUT:
+ RETVAL
+
+int
+_maxbits()
+ CODE:
+ RETVAL = BITS_PER_WORD;
+ OUTPUT:
+ RETVAL
+
+
+SV*
+sieve_primes(IN UV low, IN UV high)
+ PREINIT:
+ UV s;
+ const unsigned char* sieve;
+ AV* av = newAV();
+ CODE:
+ if (low <= high) {
+ if (get_prime_cache(high, &sieve) < high) {
+ croak("Could not generate sieve for %ld", high);
+ } else {
+ if ((low <= 2) && (high >= 2)) { av_push(av, newSVuv( 2 )); }
+ if ((low <= 3) && (high >= 3)) { av_push(av, newSVuv( 3 )); }
+ if ((low <= 5) && (high >= 5)) { av_push(av, newSVuv( 5 )); }
+ if (low < 7) { low = 7; }
+ START_DO_FOR_EACH_SIEVE_PRIME( sieve, low, high ) {
+ av_push(av,newSVuv(p));
+ } END_DO_FOR_EACH_SIEVE_PRIME
+ }
+ }
+ RETVAL = newRV_noinc( (SV*) av );
+ OUTPUT:
+ RETVAL
+
+
+SV*
+trial_primes(IN UV low, IN UV high)
+ PREINIT:
+ UV curprime;
+ AV* av = newAV();
+ CODE:
+ if (low <= high) {
+ if (low >= 2) low--; /* Make sure low gets included */
+ curprime = next_trial_prime(low);
+ while (curprime <= high) {
+ av_push(av,newSVuv(curprime));
+ curprime = next_trial_prime(curprime);
+ }
+ }
+ RETVAL = newRV_noinc( (SV*) av );
+ OUTPUT:
+ RETVAL
+
+SV*
+segment_primes(IN UV low, IN UV high, IN UV segment_size = 65536UL)
+ PREINIT:
+ AV* av = newAV();
+ unsigned char* sieve;
+ CODE:
+ if ((low <= 2) && (high >= 2)) { av_push(av, newSVuv( 2 )); }
+ if ((low <= 3) && (high >= 3)) { av_push(av, newSVuv( 3 )); }
+ if ((low <= 5) && (high >= 5)) { av_push(av, newSVuv( 5 )); }
+ if (low < 7) low = 7;
+ if (low <= high) {
+ /* Call the segment siever one or more times */
+ sieve = (unsigned char*) malloc( segment_size );
+ if (sieve == 0)
+ croak("Could not allocate %lu bytes for segment sieve", segment_size);
+ while (low <= high) {
+ UV seghigh = ((high/30 - low/30) < segment_size)
+ ? high
+ : ( (low/30 + segment_size-1)*30+29 );
+ UV startd = low/30;
+ UV endd = seghigh/30;
+ UV ranged = endd - startd + 1;
+ assert(endd >= startd);
+ assert(ranged <= segment_size);
+
+ /* Sieve from startd*30+1 to endd*30+29. */
+ if (sieve_segment(sieve, startd, endd) == 0) {
+ croak("Could not segment sieve from %lu to %lu", startd*30+1, endd*30+29);
+ break;
+ }
+
+ START_DO_FOR_EACH_SIEVE_PRIME( sieve, low-startd*30, seghigh-30*startd )
+ av_push(av,newSVuv( startd*30 + p ));
+ END_DO_FOR_EACH_SIEVE_PRIME
+
+ low = seghigh+2;
+ }
+ free(sieve);
+ }
+ RETVAL = newRV_noinc( (SV*) av );
+ OUTPUT:
+ RETVAL
+
+SV*
+erat_primes(IN UV low, IN UV high)
+ PREINIT:
+ unsigned char* sieve;
+ AV* av = newAV();
+ CODE:
+ if (low <= high) {
+ sieve = sieve_erat30(high);
+ if (sieve == 0) {
+ croak("Could not generate sieve for %lu", high);
+ } else {
+ if ((low <= 2) && (high >= 2)) { av_push(av, newSVuv( 2 )); }
+ if ((low <= 3) && (high >= 3)) { av_push(av, newSVuv( 3 )); }
+ if ((low <= 5) && (high >= 5)) { av_push(av, newSVuv( 5 )); }
+ if (low < 7) { low = 7; }
+ START_DO_FOR_EACH_SIEVE_PRIME( sieve, low, high ) {
+ av_push(av,newSVuv(p));
+ } END_DO_FOR_EACH_SIEVE_PRIME
+ free(sieve);
+ }
+ }
+ RETVAL = newRV_noinc( (SV*) av );
+ OUTPUT:
+ RETVAL
+
+
+SV*
+erat_simple_primes(IN UV low, IN UV high)
+ PREINIT:
+ UV* sieve;
+ UV s;
+ AV* av = newAV();
+ CODE:
+ if (low <= high) {
+ sieve = sieve_erat(high);
+ if (sieve == 0) {
+ croak("Could not generate sieve for %ld", high);
+ } else {
+ if (low <= 2) { av_push(av, newSVuv( 2 )); low = 3; }
+ low = low/2;
+ high = (high-1)/2;
+ for (s = low; s <= high; s++) {
+ if ( ! IS_SET_ARRAY_BIT(sieve, s) ) {
+ av_push(av,newSVuv( 2*s+1 ));
+ }
+ }
+ free(sieve);
+ }
+ }
+ RETVAL = newRV_noinc( (SV*) av );
+ OUTPUT:
+ RETVAL
diff --git a/bitarray.h b/bitarray.h
new file mode 100644
index 0000000..a076296
--- /dev/null
+++ b/bitarray.h
@@ -0,0 +1,54 @@
+#ifndef MPU_BITARRAY_H
+#define MPU_BITARRAY_H
+
+#include "EXTERN.h"
+#include "perl.h"
+
+/* From perl.h, wrapped in PERL_CORE */
+#ifndef U32_CONST
+# if INTSIZE >= 4
+# define U32_CONST(x) ((U32TYPE)x##U)
+# else
+# define U32_CONST(x) ((U32TYPE)x##UL)
+# endif
+#endif
+
+/* From perl.h, wrapped in PERL_CORE */
+#ifndef U64_CONST
+# ifdef HAS_QUAD
+# if INTSIZE >= 8
+# define U64_CONST(x) ((U64TYPE)x##U)
+# elif LONGSIZE >= 8
+# define U64_CONST(x) ((U64TYPE)x##UL)
+# elif QUADKIND == QUAD_IS_LONG_LONG
+# define U64_CONST(x) ((U64TYPE)x##ULL)
+# else /* best guess we can make */
+# define U64_CONST(x) ((U64TYPE)x##UL)
+# endif
+# endif
+#endif
+
+
+#ifdef HAS_QUAD
+ #define BITS_PER_WORD 64
+ #define UVCONST(x) U64_CONST(x)
+#else
+ #define BITS_PER_WORD 32
+ #define UVCONST(x) U32_CONST(x)
+#endif
+
+#define MAXBIT (BITS_PER_WORD-1)
+#define NWORDS(bits) ( ((bits)+BITS_PER_WORD-1) / BITS_PER_WORD )
+#define NBYTES(bits) ( ((bits)+8-1) / 8 )
+
+/* The mod-30 wheel is byte-based so doesn't use these. */
+#define SET_ARRAY_BIT(ar,n) \
+ ar[(n)/BITS_PER_WORD] |= (UVCONST(1) << ((n)%BITS_PER_WORD))
+#define XOR_ARRAY_BIT(ar,n) \
+ ar[(n)/BITS_PER_WORD] ^= (UVCONST(1) << ((n)%BITS_PER_WORD))
+#define CLR_ARRAY_BIT(ar,n) \
+ ar[(n)/BITS_PER_WORD] &= ~(UVCONST(1) << ((n)%BITS_PER_WORD))
+#define IS_SET_ARRAY_BIT(ar,n) \
+ (ar[(n)/BITS_PER_WORD] & (UVCONST(1) << ((n)%BITS_PER_WORD)) )
+
+#endif
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
new file mode 100644
index 0000000..5315dbc
--- /dev/null
+++ b/lib/Math/Prime/Util.pm
@@ -0,0 +1,353 @@
+package Math::Prime::Util;
+use strict;
+use warnings;
+use Carp qw/croak confess/;
+
+BEGIN {
+ $Math::Prime::Util::AUTHORITY = 'cpan:DANAJ';
+ $Math::Prime::Util::VERSION = '0.01';
+}
+
+# parent is cleaner, and in the Perl 5.10.1 / 5.12.0 core, but not earlier.
+# use parent qw( Exporter );
+use base qw( Exporter );
+our @EXPORT_OK = qw(
+ prime_precalc prime_free
+ is_prime
+ primes
+ next_prime prev_prime
+ prime_count prime_count_lower prime_count_upper prime_count_approx
+ nth_prime nth_prime_lower nth_prime_upper nth_prime_approx
+ );
+
+BEGIN {
+ eval {
+ require XSLoader;
+ XSLoader::load(__PACKAGE__, $Math::Prime::Util::VERSION);
+ prime_precalc(0);
+ 1;
+ } or do {
+ # We could insert a Pure Perl implementation here.
+ croak "XS Code not available: $@";
+ }
+}
+
+
+sub primes {
+ my $optref = {}; $optref = shift if ref $_[0] eq 'HASH';
+ croak "no parameters to primes" unless scalar @_ > 0;
+ croak "too many parameters to primes" unless scalar @_ <= 2;
+ my $low = (@_ == 2) ? shift : 2;
+ my $high = shift;
+ my $sref = [];
+
+ # Validate parameters
+ if ( (!defined $low) || (!defined $high) ||
+ #($low < 0) || ($high < 0) ||
+ ($low =~ tr/0123456789//c) || ($high =~ tr/0123456789//c)
+ ) {
+ croak "Parameters must be positive integers";
+ }
+ return $sref if ($low > $high) || ($high < 2);
+
+ my $method = $optref->{'method'};
+ $method = 'Dynamic' unless defined $method;
+
+ if ($method =~ /^(Dyn\w*|Default|Generate)$/i) {
+ # Dynamic -- we should try to do something smart.
+
+ # Tiny range?
+ if (($low+1) >= $high) {
+ $method = 'Trial';
+
+ # Fast for cached sieve?
+ } elsif (($high <= (65536*30)) || ($high <= _get_prime_cache_size)) {
+ $method = 'Sieve';
+
+ # More memory than we should reasonably use for base sieve?
+ } elsif ($high > (32*1024*1024*30)) {
+ $method = 'Segment';
+
+ # Only want half or less of the range low-high ?
+ } elsif ( int($high / ($high-$low)) >= 2 ) {
+ $method = 'Segment';
+
+ } else {
+ $method = 'Sieve';
+ }
+ }
+
+ if ($method =~ /^Trial$/i) { $sref = trial_primes($low, $high); }
+ elsif ($method =~ /^Erat\w*$/i) { $sref = erat_primes($low, $high); }
+ elsif ($method =~ /^Simple\w*$/i) { $sref = erat_simple_primes($low, $high); }
+ elsif ($method =~ /^Seg\w*$/i) { $sref = segment_primes($low, $high); }
+ elsif ($method =~ /^Sieve$/i) { $sref = sieve_primes($low, $high); }
+ else { croak "Unknown prime method: $method"; }
+
+ #return (wantarray) ? @{$sref} : $sref;
+ return $sref;
+}
+
+
+1;
+
+__END__
+
+
+# ABSTRACT: Utilities related to prime numbers, including fast generators / sievers
+
+=pod
+
+=encoding utf8
+
+
+=head1 NAME
+
+Math::Prime::Util - Utilities related to prime numbers, including fast generators / sievers
+
+
+=head1 VERSION
+
+Version 0.01
+
+
+=head1 SYNOPSIS
+
+ use Math::Prime::Util qw/primes/;
+
+ # Get a big array reference of many primes
+ my $aref = primes( 100_000_000 );
+
+ # All the primes between 5k and 10k inclusive
+ my $aref = primes( 5_000, 10_000 );
+
+
+
+=head1 DESCRIPTION
+
+A set of utilities related to prime numbers. These include multiple sieving
+methods, is_prime, prime_count, nth_prime, approximations and bounds for
+the prime_count and nth prime, next_prime and prev_prime, factoring utilities,
+and more.
+
+The default sieving and factoring are intended to be the fastest on CPAN,
+including Math::Prime::XS, Math::Prime::FastSieve, and Math::Factor::XS.
+
+
+
+=head1 FUNCTIONS
+
+=head2 prime_precalc
+
+ prime_precalc( 1_000_000_000 );
+
+Let the module prepare for fast operation up to a specific number. It is not
+necessary to call this, but it gives you more control over when memory is
+allocated and gives faster results for multiple calls in some cases. In the
+current implementation this will calculate a sieve for all numbers up to the
+specified number.
+
+
+=head2 prime_free
+
+ prime_free;
+
+Frees any extra memory the module may have allocated. Like with
+C<prime_precalc>, it is not necessary to call this, but if you're done
+making calls, or want things cleanup up, you can use this.
+
+
+=head2 is_prime
+
+Returns true if the number is prime, false if not.
+
+ print "$n is prime" if is_prime($n);
+
+
+=head2 primes
+
+Returns all the primes between the lower and upper limits (inclusive), with
+a lower limit of C<2> if none is given.
+
+An array reference is returned (with large lists this is much faster and uses
+less memory than returning an array directly).
+
+ my $aref1 = primes( 1_000_000 );
+ my $aref2 = primes( 1_000_000_000_000, 1_000_000_001_000 );
+
+ my @primes = @{ primes( 500 ) };
+
+ print "$_\n" for (@{primes( 20, 100 )});
+
+Sieving will be done if required. The algorithm used will depend on the range
+and whether a sieve result already exists. Possibilities include trial
+division (for ranges with only one expected prime), a Sieve of Eratosthenes
+using wheel factorization, or a segmented sieve.
+
+
+=head2 next_prime
+
+ $n = next_prime($n);
+
+Returns the next prime greater than the input number.
+
+
+=head2 prev_prime
+
+ $n = prev_prime($n);
+
+Returns the prime smaller than the input number. 0 is returned if the
+input is C<2> or lower.
+
+
+=head2 prime_count
+
+ my $number_of_primes = prime_count( 1_000 );
+
+Returns the Prime Count function C<Pi(n)>. The current implementation relies
+on sieving to find the primes within the interval, so will take some time and
+memory. There are slightly faster ways to handle the sieving (e.g. maintain
+a list of counts from C<2 - j> for various C<j>, then do a segmented sieve
+between C<j> and C<n>), and for very large numbers the methods of Meissel,
+Lehmer, or Lagarias-Miller-Odlyzko-Deleglise-Rivat may be appropriate.
+
+
+=head2 prime_count_upper
+
+=head2 prime_count_lower
+
+ my $lower_limit = prime_count_lower($n);
+ die unless prime_count($n) >= $lower_limit;
+
+ my $upper_limit = prime_count_upper($n);
+ die unless prime_count($n) <= $upper_limit;
+
+Returns an upper or lower bound on the number of primes below the input number.
+These are analytical routines, so will take a fixed amount of time and no
+memory. The actual C<prime_count> will always be on or between these numbers.
+
+A common place these would be used is sizing an array to hold the first C<$n>
+primes. It may be desirable to use a bit more memory than is necessary, to
+avoid calling C<prime_count>.
+
+These routines use hand-verified tight limits below a range at least C<2^32>,
+and fall back to the proven Dusart bounds of
+
+ x/logx * (1 + 1/logx + 1.80/log^2x) <= Pi(x)
+
+ x/logx * (1 + 1/logx + 2.51/log^2x) >= Pi(x)
+
+above that range.
+
+
+=head2 prime_count_approx
+
+ print "there are about ",
+ prime_count_approx( 10 ** 18 ),
+ " primes below one quintillion.\n";
+
+Returns an approximation to the C<prime_count> function, without having to
+generate any primes. The results are very close for small numbers, but less
+so with large ranges. The current implementation is 0.00033% too small
+for the example, but takes under a microsecond and no memory to get the
+result.
+
+
+=head2 nth_prime
+
+ say "The ten thousandth prime is ", nth_prime(10_000);
+
+Returns the prime that lies in index C<n> in the array of prime numbers. Put
+another way, this returns the smallest C<p> such that C<Pi(p) E<gt>= n>.
+
+This relies on generating primes, so can require a lot of time and space for
+large inputs.
+
+
+=head2 nth_prime_upper
+
+=head2 nth_prime_lower
+
+ my $lower_limit = nth_prime_lower($n);
+ die unless nth_prime($n) >= $lower_limit;
+
+ my $upper_limit = nth_prime_upper($n);
+ die unless nth_prime($n) <= $upper_limit;
+
+Returns an analytical upper or lower bound on the Nth prime. This will be
+very fast. The lower limit uses the Dusart 1999 bounds for all C<n>, while
+the upper bound uses one of the two Dusart 1999 bounds for C<n E<gt>= 27076>,
+the Robin 1983 bound for C<n E<gt>= 7022>, and the simple bound of
+C<n * (logn + loglogn)> for C<n E<lt> 7022>.
+
+
+=head2 nth_prime_approx
+
+ say "The one trillionth prime is ~ ", nth_prime_approx(10**12);
+
+Returns an approximation to the C<nth_prime> function, without having to
+generate any primes. Uses the Cipolla 1902 approximation with two
+polynomials, plus a correction term for small values to reduce the error.
+
+
+=head1 LIMITATIONS
+
+Some of the functions have not yet transitioned to using a segmented sieve,
+so will use too much memory to be practical when called with very large
+numbers (C<10^11> and up).
+
+I have not completed testing all the functions near the word size limit
+(e.g. C<2^32> for 32-bit machines). Please report any problems you find.
+
+
+=head1 PERFORMANCE
+
+Counting the primes to C<10^10> (10 billion), with time in seconds.
+Pi(10^10) = 455,052,511.
+
+ 1.9 primesieve 3.6 forced to use only a single thread
+ 5.6 Tomás Oliveira e Silva's segmented sieve v2 (Sep 2010)
+ 6.6 primegen (optimized Sieve of Atkin)
+ 11.2 Tomás Oliveira e Silva's segmented sieve v1 (May 2003)
+
+ 15.6 My Sieve of Eratosthenes using a mod-30 wheel
+ 17.2 A slightly modified verion of Terje Mathisen's mod-30 sieve
+ 35.5 Basic Sieve of Eratosthenes on odd numbers
+ 33.4 Sieve of Atkin, from Praxis (not correct)
+ 72.8 Sieve of Atkin, 10-minute fixup of basic algorithm
+ 91.6 Sieve of Atkin, Wikipedia-like
+
+Perl modules, counting the primes to C<800_000_000> (800 million), in seconds:
+
+ 0.9 Math::Prime::Util 0.01
+ 2.9 Math::Prime::FastSieve 0.12
+ 11.7 Math::Prime::XS 0.29
+ 15.0 Bit::Vector 7.2
+ (hours) Math::Primality 0.04
+
+
+=head1 AUTHORS
+
+Dana Jacobsen E<lt>dana at acm.orgE<gt>
+
+
+=head1 ACKNOWLEDGEMENTS
+
+Eratosthenes of Cyrene provided the elegant and simple algorithm for finding
+the primes.
+
+Terje Mathisen, A.R. Quesada, and B. Van Pelt all had useful ideas which I
+used in my wheel sieve.
+
+Tomás Oliveira e Silva has released the source for a very fast segmented sieve.
+The current implementation does not use these ideas, but future versions likely
+will.
+
+
+=head1 COPYRIGHT
+
+Copyright 2011-2012 by Dana Jacobsen E<lt>dana at acm.orgE<gt>
+
+This program is free software; you can redistribute it and/or modify it under the same terms as Perl itself.
+
+=cut
diff --git a/sieve.c b/sieve.c
new file mode 100644
index 0000000..b3765ef
--- /dev/null
+++ b/sieve.c
@@ -0,0 +1,250 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <assert.h>
+#include <limits.h>
+#include <math.h>
+
+#include "sieve.h"
+#include "bitarray.h"
+
+
+static unsigned char* prime_cache_sieve = 0;
+static UV prime_cache_size = 0;
+
+/* Get the maximum sieved value of the cached prime sieve. */
+UV get_prime_cache_size(void) { return prime_cache_size; }
+
+/*
+ * Get the size and a pointer to the cached prime sieve.
+ * Returns the maximum sieved value in the sieve.
+ * Allocates and sieves if needed.
+ *
+ * The sieve holds 30 numbers per byte, using a mod-30 wheel.
+ */
+UV get_prime_cache(UV n, const unsigned char** sieve)
+{
+ if (prime_cache_size < n) {
+
+ if (prime_cache_sieve != 0)
+ free(prime_cache_sieve);
+ prime_cache_sieve = 0;
+ prime_cache_size = 0;
+
+ /* Sieve a bit more than asked, to mitigate thrashing */
+ if (n < (UV_MAX-3840))
+ n += 3840;
+ /* TODO: testing near 2^32-1 */
+
+ prime_cache_sieve = sieve_erat30(n);
+
+ if (prime_cache_sieve != 0)
+ prime_cache_size = n;
+ }
+
+ if (sieve != 0)
+ *sieve = prime_cache_sieve;
+ return prime_cache_size;
+}
+
+
+
+void prime_precalc(UV n)
+{
+ if ( (n == 0) && (prime_cache_sieve == 0) ) {
+ /* On initialization, make a few primes (2-30k using 1k memory) */
+ UV initial_primes_to = 30 * (1024-8);
+ prime_cache_sieve = sieve_erat30(initial_primes_to);
+ if (prime_cache_sieve != 0)
+ prime_cache_size = initial_primes_to;
+ return;
+ }
+
+ get_prime_cache(n, 0); /* Sieve to n */
+}
+
+void prime_free(void)
+{
+ if (prime_cache_sieve != 0)
+ free(prime_cache_sieve);
+ prime_cache_sieve = 0;
+ prime_cache_size = 0;
+
+ prime_precalc(0);
+}
+
+
+
+UV* sieve_erat(UV end)
+{
+ UV* mem;
+ size_t n, s;
+ size_t last = (end+1)/2;
+
+ mem = (UV*) calloc( NWORDS(last), sizeof(UV) );
+ if (mem == 0) {
+ croak("allocation failure in sieve_erat: could not alloc %lu bits", (unsigned long)last);
+ return 0;
+ }
+
+ n = 3;
+ /* TODO: overflow */
+ while ( (n*n) <= end) {
+ for (s = n*n; s <= end; s += 2*n)
+ SET_ARRAY_BIT(mem,s/2);
+ do { n += 2; } while (IS_SET_ARRAY_BIT(mem,n/2));
+ }
+
+ SET_ARRAY_BIT(mem, 1/2); /* 1 is composite */
+
+ return mem;
+}
+
+
+/* Wheel 30 sieve. Ideas from Terje Mathisen and Quesada / Van Pelt. */
+unsigned char* sieve_erat30(UV end)
+{
+ unsigned char* mem;
+ size_t max_buf, buffer_words, limit;
+ UV prime;
+
+ max_buf = (end/30) + ((end%30) != 0);
+ buffer_words = (max_buf + sizeof(UV) - 1) / sizeof(UV);
+ mem = (unsigned char*) calloc( buffer_words, sizeof(UV) );
+ if (mem == 0) {
+ croak("allocation failure in sieve_erat30: could not alloc %lu bytes", (unsigned long)(buffer_words*sizeof(UV)));
+ return 0;
+ }
+
+ /* Shortcut to mark 7. Just an optimization. */
+ if ( (7*7) <= end ) {
+ UV d = 1;
+ while ( (d+6) < max_buf) {
+ mem[d+0] = 0x20; mem[d+1] = 0x10; mem[d+2] = 0x81; mem[d+3] = 0x08;
+ mem[d+4] = 0x04; mem[d+5] = 0x40; mem[d+6] = 0x02; d += 7;
+ }
+ if ( d < max_buf ) mem[d++] = 0x20;
+ if ( d < max_buf ) mem[d++] = 0x10;
+ if ( d < max_buf ) mem[d++] = 0x81;
+ if ( d < max_buf ) mem[d++] = 0x08;
+ if ( d < max_buf ) mem[d++] = 0x04;
+ if ( d < max_buf ) mem[d++] = 0x40;
+ assert(d >= max_buf);
+ }
+ limit = sqrt((double) end); /* prime*prime can overflow */
+ for (prime = 11; prime <= limit; prime = next_prime_in_sieve(mem,prime)) {
+ UV d = (prime*prime)/30;
+ UV m = (prime*prime) - d*30;
+ UV dinc = (2*prime)/30;
+ UV minc = (2*prime) - dinc*30;
+ UV wdinc[8];
+ unsigned char wmask[8];
+ int i;
+
+ /* Find the positions of the next composites we will mark */
+ for (i = 1; i <= 8; i++) {
+ UV dlast = d;
+ do {
+ d += dinc;
+ m += minc;
+ if (m >= 30) { d++; m -= 30; }
+ } while ( masktab30[m] == 0 );
+ wdinc[i-1] = d - dlast;
+ wmask[i%8] = masktab30[m];
+ }
+ d -= prime;
+#if 0
+ assert(d == ((prime*prime)/30));
+ assert(d < max_buf);
+ assert(prime = (wdinc[0]+wdinc[1]+wdinc[2]+wdinc[3]+wdinc[4]+wdinc[5]+wdinc[6]+wdinc[7]));
+#endif
+ i = 0; /* Mark the composites */
+ do {
+ mem[d] |= wmask[i];
+ d += wdinc[i];
+ i = (i+1) & 7;
+ } while (d < max_buf);
+ }
+
+ mem[0] |= masktab30[1]; /* 1 is composite */
+
+ return mem;
+}
+
+
+
+int sieve_segment(unsigned char* mem, UV startd, UV endd)
+{
+ const unsigned char* sieve;
+ UV limit;
+ UV pcsize;
+ UV startp = 30*startd;
+ UV endp = 30*endd+29;
+ UV ranged = endd - startd + 1;
+
+ assert(mem != 0);
+ assert(endd >= startd);
+ memset(mem, 0, ranged);
+
+ limit = sqrt( (double) endp ) + 1;
+ /* printf("segment sieve from %lu to %lu (aux sieve to %lu)\n", startp, endp, limit); */
+ pcsize = get_prime_cache(limit, &sieve);
+ if (pcsize < limit) {
+ croak("Couldn't generate small sieve for segment sieve");
+ return 0;
+ }
+
+ START_DO_FOR_EACH_SIEVE_PRIME(sieve, 7, pcsize)
+ {
+ /* p increments from 7 to at least sqrt(endp) */
+ UV p2 = p*p; /* TODO: overflow */
+ if (p2 >= endp) break;
+ /* Find first multiple of p greater than p*p and larger than startp */
+ if (p2 < startp) {
+ p2 = (startp / p) * p;
+ if (p2 < startp) p2 += p;
+ }
+ /* Bump to next multiple that isn't divisible by 2, 3, or 5 */
+ while (masktab30[p2%30] == 0) { p2 += p; }
+ if (p2 < endp) {
+ /* Sieve from startd to endd starting at p2, stepping p */
+#if 0
+ /* Basic sieve */
+ do {
+ mem[(p2 - startp)/30] |= masktab30[p2%30];
+ do { p2 += 2*p; } while (masktab30[p2%30] == 0);
+ } while (p2 < endp);
+#else
+ UV d = (p2)/30;
+ UV m = (p2) - d*30;
+ UV dinc = (2*p)/30;
+ UV minc = (2*p) - dinc*30;
+ UV wdinc[8];
+ unsigned char wmask[8];
+ int i;
+
+ /* Find the positions of the next composites we will mark */
+ for (i = 1; i <= 8; i++) {
+ UV dlast = d;
+ do {
+ d += dinc;
+ m += minc;
+ if (m >= 30) { d++; m -= 30; }
+ } while ( masktab30[m] == 0 );
+ wdinc[i-1] = d - dlast;
+ wmask[i%8] = masktab30[m];
+ }
+ d -= p;
+ i = 0; /* Mark the composites */
+ do {
+ mem[d-startd] |= wmask[i];
+ d += wdinc[i];
+ i = (i+1) & 7;
+ } while (d <= endd);
+#endif
+ }
+ }
+ END_DO_FOR_EACH_SIEVE_PRIME;
+
+ return 1;
+}
diff --git a/sieve.h b/sieve.h
new file mode 100644
index 0000000..5c08bc4
--- /dev/null
+++ b/sieve.h
@@ -0,0 +1,85 @@
+#ifndef MPU_SIEVE_H
+#define MPU_SIEVE_H
+
+#include "EXTERN.h"
+#include "perl.h"
+
+extern UV get_prime_cache_size(void);
+extern UV get_prime_cache(UV n, const unsigned char** sieve);
+
+extern void prime_precalc(UV x);
+extern void prime_free(void);
+extern UV* sieve_erat(UV end);
+extern unsigned char* sieve_erat30(UV end);
+extern int sieve_segment(unsigned char* mem, UV startd, UV endd);
+
+
+static const UV wheel30[] = {1, 7, 11, 13, 17, 19, 23, 29};
+/* Used for moving between primes */
+static unsigned char nextwheel30[30] = {
+ 1, 7, 7, 7, 7, 7, 7, 11, 11, 11, 11, 13, 13, 17, 17,
+ 17, 17, 19, 19, 23, 23, 23, 23, 29, 29, 29, 29, 29, 29, 1 };
+static unsigned char prevwheel30[30] = {
+ 29, 29, 1, 1, 1, 1, 1, 1, 7, 7, 7, 7, 11, 11, 13,
+ 13, 13, 13, 17, 17, 19, 19, 19, 19, 23, 23, 23, 23, 23, 23 };
+/* The bit mask within a byte */
+static unsigned char masktab30[30] = {
+ 0, 1, 0, 0, 0, 0, 0, 2, 0, 0, 0, 4, 0, 8, 0,
+ 0, 0, 16, 0, 32, 0, 0, 0, 64, 0, 0, 0, 0, 0,128 };
+/* Add this to a number and you'll ensure you're on a wheel location */
+static unsigned char distancewheel30[30] = {
+ 1, 0, 5, 4, 3, 2, 1, 0, 3, 2, 1, 0, 1, 0, 3,
+ 2, 1, 0, 1, 0, 3, 2, 1, 0, 5, 4, 3, 2, 1, 0 };
+
+static int is_prime_in_sieve(const unsigned char* sieve, UV p) {
+ UV d = p/30;
+ UV m = p - d*30;
+ /* If m isn't part of the wheel, we return 0 */
+ return ( (masktab30[m] != 0) && ((sieve[d] & masktab30[m]) == 0) );
+}
+
+/* Warning -- can go off the end of the sieve */
+static UV next_prime_in_sieve(const unsigned char* sieve, UV p) {
+ UV d, m;
+ if (p < 7)
+ return (p < 2) ? 2 : (p < 3) ? 3 : (p < 5) ? 5 : 7;
+ d = p/30;
+ m = p - d*30;
+ do {
+ if (m==29) { d++; m = 1; while (sieve[d] == 0xFF) d++; }
+ else { m = nextwheel30[m]; }
+ } while (sieve[d] & masktab30[m]);
+ return(d*30+m);
+}
+static UV prev_prime_in_sieve(const unsigned char* sieve, UV p) {
+ UV d, m;
+ if (p <= 7)
+ return (p <= 2) ? 0 : (p <= 3) ? 2 : (p <= 5) ? 3 : 5;
+ d = p/30;
+ m = p - d*30;
+ do {
+ m = prevwheel30[m]; if (m==29) { if (d == 0) return 0; d--; }
+ } while (sieve[d] & masktab30[m]);
+ return(d*30+m);
+}
+
+/* Useful macros for the wheel-30 sieve array */
+#define START_DO_FOR_EACH_SIEVE_PRIME(sieve, a, b) \
+ { \
+ UV p = a; \
+ UV l_ = b; \
+ UV d_ = p/30; \
+ UV m_ = p-d_*30; \
+ m_ += distancewheel30[m_]; \
+ p = d_*30 + m_; \
+ while ( p <= l_ ) { \
+ if ((sieve[d_] & masktab30[m_]) == 0)
+
+#define END_DO_FOR_EACH_SIEVE_PRIME \
+ m_ = nextwheel30[m_]; if (m_ == 1) { d_++; } \
+ p = d_*30+m_; \
+ } \
+ }
+
+
+#endif
diff --git a/t/01-load.t b/t/01-load.t
new file mode 100644
index 0000000..5a3eeda
--- /dev/null
+++ b/t/01-load.t
@@ -0,0 +1,7 @@
+#!/usr/bin/env perl
+use strict;
+use warnings;
+
+use Test::More tests => 1;
+
+require_ok 'Math::Prime::Util';
diff --git a/t/02-can.t b/t/02-can.t
new file mode 100644
index 0000000..34748ff
--- /dev/null
+++ b/t/02-can.t
@@ -0,0 +1,16 @@
+#!/usr/bin/env perl
+use strict;
+use warnings;
+use Math::Prime::Util;
+
+use Test::More tests => 1;
+
+my @functions = qw(
+ prime_precalc prime_free
+ is_prime
+ primes
+ next_prime prev_prime
+ prime_count prime_count_lower prime_count_upper prime_count_approx
+ nth_prime nth_prime_lower nth_prime_upper nth_prime_approx
+ );
+can_ok( 'Math::Prime::Util', @functions);
diff --git a/t/03-init.t b/t/03-init.t
new file mode 100644
index 0000000..aad681b
--- /dev/null
+++ b/t/03-init.t
@@ -0,0 +1,22 @@
+#!/usr/bin/env perl
+use strict;
+use warnings;
+use Math::Prime::Util qw/prime_precalc prime_free/;
+
+use Test::More tests => 3;
+
+
+# How do we test these? Use some private functions and assume things about
+# internal behavior. Not the best idea.
+
+can_ok( 'Math::Prime::Util', '_get_prime_cache_size' );
+
+my $init_size = Math::Prime::Util::_get_prime_cache_size;
+
+prime_precalc(10_000_000);
+
+cmp_ok( Math::Prime::Util::_get_prime_cache_size, '>', $init_size, "Internal space grew after large precalc" );
+
+prime_free;
+
+is( Math::Prime::Util::_get_prime_cache_size, $init_size, "Internal space went back to original size after free" );
diff --git a/t/10-isprime.t b/t/10-isprime.t
new file mode 100644
index 0000000..3e96b1d
--- /dev/null
+++ b/t/10-isprime.t
@@ -0,0 +1,112 @@
+#!/usr/bin/env perl
+use strict;
+use warnings;
+
+use Test::More;
+use Math::Prime::Util qw/is_prime/;
+
+my $use64 = Math::Prime::Util::_maxbits > 32;
+my $extra = defined $ENV{RELEASE_TESTING} && $ENV{RELEASE_TESTING};
+
+plan tests => 6 + 19 + 3573 + (5 + 29 + 22 + 23 + 16) + 15 + 27
+ + ($use64 ? 5+1 : 0)
+ + ($extra ? 6 : 0)
+ + (($extra && $use64) ? 19 : 0);
+
+# Some of these tests were inspired by Math::Primality's tests
+
+ok( is_prime(2), '2 is prime');
+ok(!is_prime(1), '1 is not prime');
+ok(!is_prime(0), '0 is not prime');
+ok(!is_prime(-1), '-1 is not prime');
+ok(!is_prime(-2), '-2 is not prime');
+ok(!is_prime(20), '20 is not prime');
+
+# powers of 2
+foreach my $k (2 .. 20) {
+ my $k2 = 2**$k;
+ ok(!is_prime($k2), "2**$k=$k2 is not prime");
+}
+
+my @small_primes = qw/
+2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71
+73 79 83 89 97 101 103 107 109 113 127 131 137 139 149 151 157 163 167 173
+179 181 191 193 197 199 211 223 227 229 233 239 241 251 257 263 269 271 277 281
+283 293 307 311 313 317 331 337 347 349 353 359 367 373 379 383 389 397 401 409
+419 421 431 433 439 443 449 457 461 463 467 479 487 491 499 503 509 521 523 541
+547 557 563 569 571 577 587 593 599 601 607 613 617 619 631 641 643 647 653 659
+661 673 677 683 691 701 709 719 727 733 739 743 751 757 761 769 773 787 797 809
+811 821 823 827 829 839 853 857 859 863 877 881 883 887 907 911 919 929 937 941
+947 953 967 971 977 983 991 997 1009 1013 1019 1021 1031 1033 1039 1049 1051 1061 1063 1069
+1087 1091 1093 1097 1103 1109 1117 1123 1129 1151 1153 1163 1171 1181 1187 1193 1201 1213 1217 1223
+1229 1231 1237 1249 1259 1277 1279 1283 1289 1291 1297 1301 1303 1307 1319 1321 1327 1361 1367 1373
+1381 1399 1409 1423 1427 1429 1433 1439 1447 1451 1453 1459 1471 1481 1483 1487 1489 1493 1499 1511
+1523 1531 1543 1549 1553 1559 1567 1571 1579 1583 1597 1601 1607 1609 1613 1619 1621 1627 1637 1657
+1663 1667 1669 1693 1697 1699 1709 1721 1723 1733 1741 1747 1753 1759 1777 1783 1787 1789 1801 1811
+1823 1831 1847 1861 1867 1871 1873 1877 1879 1889 1901 1907 1913 1931 1933 1949 1951 1973 1979 1987
+1993 1997 1999 2003 2011 2017 2027 2029 2039 2053 2063 2069 2081 2083 2087 2089 2099 2111 2113 2129
+2131 2137 2141 2143 2153 2161 2179 2203 2207 2213 2221 2237 2239 2243 2251 2267 2269 2273 2281 2287
+2293 2297 2309 2311 2333 2339 2341 2347 2351 2357 2371 2377 2381 2383 2389 2393 2399 2411 2417 2423
+2437 2441 2447 2459 2467 2473 2477 2503 2521 2531 2539 2543 2549 2551 2557 2579 2591 2593 2609 2617
+2621 2633 2647 2657 2659 2663 2671 2677 2683 2687 2689 2693 2699 2707 2711 2713 2719 2729 2731 2741
+2749 2753 2767 2777 2789 2791 2797 2801 2803 2819 2833 2837 2843 2851 2857 2861 2879 2887 2897 2903
+2909 2917 2927 2939 2953 2957 2963 2969 2971 2999 3001 3011 3019 3023 3037 3041 3049 3061 3067 3079
+3083 3089 3109 3119 3121 3137 3163 3167 3169 3181 3187 3191 3203 3209 3217 3221 3229 3251 3253 3257
+3259 3271 3299 3301 3307 3313 3319 3323 3329 3331 3343 3347 3359 3361 3371 3373 3389 3391 3407 3413
+3433 3449 3457 3461 3463 3467 3469 3491 3499 3511 3517 3527 3529 3533 3539 3541 3547 3557 3559 3571
+/;
+my %small_primes;
+map { $small_primes{$_} = 1; } @small_primes;
+
+foreach my $n (0 .. 3572) {
+ if (defined $small_primes{$n}) {
+ ok(is_prime($n), "$n is prime");
+ } else {
+ ok(!is_prime($n), "$n is not prime");
+ }
+}
+
+map { ok(!is_prime($_), "A006945 number $_ is not prime") }
+ qw/9 2047 1373653 25326001 3215031751/;
+
+map { ok(!is_prime($_), "A006945 number $_ is not prime") }
+ qw/2152302898747 3474749660383 341550071728321 341550071728321 3825123056546413051/ if $use64;
+
+map { ok(!is_prime($_), "Carmichael Number $_ is not prime") }
+ qw/561 1105 1729 2465 2821 6601 8911 10585 15841 29341 41041 46657 52633
+ 62745 63973 75361 101101 340561 488881 852841 1857241 6733693
+ 9439201 17236801 23382529 34657141 56052361 146843929 216821881/;
+
+map { ok(!is_prime($_), "Pseudoprime (base 2) $_ is not prime" ) }
+ qw/341 561 645 1105 1387 1729 1905 2047 2465 2701 2821 3277 4033 4369 4371
+ 4681 5461 6601 7957 8321 52633 88357/;
+
+map { ok(!is_prime($_), "Pseudoprime (base 3) $_ is not prime" ) }
+ qw/121 703 1891 3281 8401 8911 10585 12403 16531 18721 19345 23521 31621
+ 44287 47197 55969 63139 74593 79003 82513 87913 88573 97567/;
+
+map { ok(!is_prime($_), "Pseudoprime (base 5) $_ is not prime" ) }
+ qw/781 1541 5461 5611 7813 13021 14981 15751 24211 25351 29539 38081
+ 40501 44801 53971 79381/;
+
+map { ok(is_prime($_), "Primegap start $_ is prime" ) }
+ qw/2 3 7 23 89 113 523 887 1129 1327 9551 15683 19609 31397 155921/;
+
+map { ok(is_prime($_), "Primegap end $_ is prime" ) }
+ qw/5 11 29 97 127 541 907 1151 1361 9587 15727 19661 31469 156007 360749
+ 370373 492227 1349651 1357333 2010881 4652507 17051887 20831533 47326913
+ 122164969 189695893 191913031/;
+
+map { ok(is_prime($_), "Primegap end $_ is prime" ) }
+ qw/10726905041/ if $use64;
+
+map { ok(is_prime($_), "Primegap end $_ is prime" ) }
+ qw/387096383 436273291 1294268779 1453168433 2300942869 3842611109/ if $extra;
+
+map { ok(is_prime($_), "Primegap end $_ is prime" ) }
+ qw/4302407713 20678048681 22367085353
+ 25056082543 42652618807 127976334671 182226896239 241160624143
+ 297501075799 303371455241 304599508537 416608695821 461690510011
+ 614487453523 738832927927 1346294310749 1408695493609 1968188556461
+ 2614941710599/
+ if $use64 && $extra;
diff --git a/t/11-primes.t b/t/11-primes.t
new file mode 100644
index 0000000..8fe1ac8
--- /dev/null
+++ b/t/11-primes.t
@@ -0,0 +1,96 @@
+#!/usr/bin/env perl
+use strict;
+use warnings;
+
+use Test::More;
+use Math::Prime::Util qw/primes prime_count/;
+
+my $use64 = Math::Prime::Util::_maxbits > 32;
+my $extra = defined $ENV{RELEASE_TESTING} && $ENV{RELEASE_TESTING};
+
+plan tests => 12 + 1 + 19 + ($use64 ? 1 : 0) + 1;
+
+my @small_primes = qw/
+2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71
+73 79 83 89 97 101 103 107 109 113 127 131 137 139 149 151 157 163 167 173
+179 181 191 193 197 199 211 223 227 229 233 239 241 251 257 263 269 271 277 281
+283 293 307 311 313 317 331 337 347 349 353 359 367 373 379 383 389 397 401 409
+419 421 431 433 439 443 449 457 461 463 467 479 487 491 499 503 509 521 523 541
+547 557 563 569 571 577 587 593 599 601 607 613 617 619 631 641 643 647 653 659
+661 673 677 683 691 701 709 719 727 733 739 743 751 757 761 769 773 787 797 809
+811 821 823 827 829 839 853 857 859 863 877 881 883 887 907 911 919 929 937 941
+947 953 967 971 977 983 991 997 1009 1013 1019 1021 1031 1033 1039 1049 1051 1061 1063 1069
+1087 1091 1093 1097 1103 1109 1117 1123 1129 1151 1153 1163 1171 1181 1187 1193 1201 1213 1217 1223
+1229 1231 1237 1249 1259 1277 1279 1283 1289 1291 1297 1301 1303 1307 1319 1321 1327 1361 1367 1373
+1381 1399 1409 1423 1427 1429 1433 1439 1447 1451 1453 1459 1471 1481 1483 1487 1489 1493 1499 1511
+1523 1531 1543 1549 1553 1559 1567 1571 1579 1583 1597 1601 1607 1609 1613 1619 1621 1627 1637 1657
+1663 1667 1669 1693 1697 1699 1709 1721 1723 1733 1741 1747 1753 1759 1777 1783 1787 1789 1801 1811
+1823 1831 1847 1861 1867 1871 1873 1877 1879 1889 1901 1907 1913 1931 1933 1949 1951 1973 1979 1987
+1993 1997 1999 2003 2011 2017 2027 2029 2039 2053 2063 2069 2081 2083 2087 2089 2099 2111 2113 2129
+2131 2137 2141 2143 2153 2161 2179 2203 2207 2213 2221 2237 2239 2243 2251 2267 2269 2273 2281 2287
+2293 2297 2309 2311 2333 2339 2341 2347 2351 2357 2371 2377 2381 2383 2389 2393 2399 2411 2417 2423
+2437 2441 2447 2459 2467 2473 2477 2503 2521 2531 2539 2543 2549 2551 2557 2579 2591 2593 2609 2617
+2621 2633 2647 2657 2659 2663 2671 2677 2683 2687 2689 2693 2699 2707 2711 2713 2719 2729 2731 2741
+2749 2753 2767 2777 2789 2791 2797 2801 2803 2819 2833 2837 2843 2851 2857 2861 2879 2887 2897 2903
+2909 2917 2927 2939 2953 2957 2963 2969 2971 2999 3001 3011 3019 3023 3037 3041 3049 3061 3067 3079
+3083 3089 3109 3119 3121 3137 3163 3167 3169 3181 3187 3191 3203 3209 3217 3221 3229 3251 3253 3257
+3259 3271 3299 3301 3307 3313 3319 3323 3329 3331 3343 3347 3359 3361 3371 3373 3389 3391 3407 3413
+3433 3449 3457 3461 3463 3467 3469 3491 3499 3511 3517 3527 3529 3533 3539 3541 3547 3557 3559 3571
+/;
+
+my %small_single = (
+ 0 => [],
+ 1 => [],
+ 2 => [2],
+ 3 => [2, 3],
+ 4 => [2, 3],
+ 5 => [2, 3, 5],
+ 6 => [2, 3, 5],
+ 7 => [2, 3, 5, 7],
+ 11 => [2, 3, 5, 7, 11],
+ 18 => [2, 3, 5, 7, 11, 13, 17],
+ 19 => [2, 3, 5, 7, 11, 13, 17, 19],
+ 20 => [2, 3, 5, 7, 11, 13, 17, 19],
+);
+
+while (my($high, $expect) = each (%small_single)) {
+ is_deeply( primes($high), $expect, "primes($high) should return [@{$expect}]");
+}
+
+is_deeply( primes(0, 3572), \@small_primes, "Primes between 0 and 3572" );
+
+my %small_range = (
+ "3 to 9" => [3,5,7],
+ "2 to 20" => [2,3,5,7,11,13,17,19],
+ "30 to 70" => [31,37,41,43,47,53,59,61,67],
+ "70 to 30" => [],
+ "20 to 2" => [],
+ "2 to 2" => [2],
+ "3 to 3" => [3],
+ "2 to 3" => [2,3],
+ "2 to 5" => [2,3,5],
+ "3 to 6" => [3,5],
+ "3 to 7" => [3,5,7],
+ "4 to 8" => [5,7],
+ "2010733 to 2010881" => [2010733,2010881],
+ "2010734 to 2010880" => [],
+ "3088 to 3164" => [3089,3109,3119,3121,3137,3163],
+ "3089 to 3163" => [3089,3109,3119,3121,3137,3163],
+ "3090 to 3162" => [3109,3119,3121,3137],
+ "3842610773 to 3842611109" => [3842610773,3842611109],
+ "3842610774 to 3842611108" => [],
+);
+
+while (my($range, $expect) = each (%small_range)) {
+ my($low,$high) = $range =~ /(\d+) to (\d+)/;
+ is_deeply( primes($low, $high), $expect, "primes($low,$high) should return [@{$expect}]");
+}
+
+if ($use64) {
+ is_deeply( primes(1_693_182_318_746_371, 1_693_182_318_747_671),
+ [qw/1693182318746371 1693182318747503 1693182318747523
+ 1693182318747553 1693182318747583 1693182318747613
+ 1693182318747631 1693182318747637/], "Primes between 1_693_182_318_746_371 and 1_693_182_318_747_671");
+}
+
+is( scalar @{primes(474973,838390)}, prime_count(838390) - prime_count(474973), "count primes within a range" );
diff --git a/t/12-nextprime.t b/t/12-nextprime.t
new file mode 100644
index 0000000..27fb925
--- /dev/null
+++ b/t/12-nextprime.t
@@ -0,0 +1,64 @@
+#!/usr/bin/env perl
+use strict;
+use warnings;
+
+use Test::More;
+use Math::Prime::Util qw/next_prime prev_prime/;
+
+my $use64 = Math::Prime::Util::_maxbits > 32;
+my $extra = defined $ENV{RELEASE_TESTING} && $ENV{RELEASE_TESTING};
+
+plan tests => 499*2 + 3*2 + 6;
+
+my @small_primes = qw/
+2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71
+73 79 83 89 97 101 103 107 109 113 127 131 137 139 149 151 157 163 167 173
+179 181 191 193 197 199 211 223 227 229 233 239 241 251 257 263 269 271 277 281
+283 293 307 311 313 317 331 337 347 349 353 359 367 373 379 383 389 397 401 409
+419 421 431 433 439 443 449 457 461 463 467 479 487 491 499 503 509 521 523 541
+547 557 563 569 571 577 587 593 599 601 607 613 617 619 631 641 643 647 653 659
+661 673 677 683 691 701 709 719 727 733 739 743 751 757 761 769 773 787 797 809
+811 821 823 827 829 839 853 857 859 863 877 881 883 887 907 911 919 929 937 941
+947 953 967 971 977 983 991 997 1009 1013 1019 1021 1031 1033 1039 1049 1051 1061 1063 1069
+1087 1091 1093 1097 1103 1109 1117 1123 1129 1151 1153 1163 1171 1181 1187 1193 1201 1213 1217 1223
+1229 1231 1237 1249 1259 1277 1279 1283 1289 1291 1297 1301 1303 1307 1319 1321 1327 1361 1367 1373
+1381 1399 1409 1423 1427 1429 1433 1439 1447 1451 1453 1459 1471 1481 1483 1487 1489 1493 1499 1511
+1523 1531 1543 1549 1553 1559 1567 1571 1579 1583 1597 1601 1607 1609 1613 1619 1621 1627 1637 1657
+1663 1667 1669 1693 1697 1699 1709 1721 1723 1733 1741 1747 1753 1759 1777 1783 1787 1789 1801 1811
+1823 1831 1847 1861 1867 1871 1873 1877 1879 1889 1901 1907 1913 1931 1933 1949 1951 1973 1979 1987
+1993 1997 1999 2003 2011 2017 2027 2029 2039 2053 2063 2069 2081 2083 2087 2089 2099 2111 2113 2129
+2131 2137 2141 2143 2153 2161 2179 2203 2207 2213 2221 2237 2239 2243 2251 2267 2269 2273 2281 2287
+2293 2297 2309 2311 2333 2339 2341 2347 2351 2357 2371 2377 2381 2383 2389 2393 2399 2411 2417 2423
+2437 2441 2447 2459 2467 2473 2477 2503 2521 2531 2539 2543 2549 2551 2557 2579 2591 2593 2609 2617
+2621 2633 2647 2657 2659 2663 2671 2677 2683 2687 2689 2693 2699 2707 2711 2713 2719 2729 2731 2741
+2749 2753 2767 2777 2789 2791 2797 2801 2803 2819 2833 2837 2843 2851 2857 2861 2879 2887 2897 2903
+2909 2917 2927 2939 2953 2957 2963 2969 2971 2999 3001 3011 3019 3023 3037 3041 3049 3061 3067 3079
+3083 3089 3109 3119 3121 3137 3163 3167 3169 3181 3187 3191 3203 3209 3217 3221 3229 3251 3253 3257
+3259 3271 3299 3301 3307 3313 3319 3323 3329 3331 3343 3347 3359 3361 3371 3373 3389 3391 3407 3413
+3433 3449 3457 3461 3463 3467 3469 3491 3499 3511 3517 3527 3529 3533 3539 3541 3547 3557 3559 3571
+/;
+
+for (my $i = 0; $i < (scalar @small_primes) - 1; $i++) {
+ my $n = next_prime($small_primes[$i]);
+ is($n, $small_primes[$i+1], "the next prime after $small_primes[$i] is $small_primes[$i+1] ?= $n");
+ my $p = prev_prime($small_primes[$i+1]);
+ is($p, $small_primes[$i], "the prev prime before $small_primes[$i+1] is $small_primes[$i] ?= $n");
+}
+
+my %primegaps = (
+ 19609 => 52,
+ 360653 => 96,
+ 2010733 => 148,
+);
+
+while (my($base, $range) = each (%primegaps)) {
+ is( next_prime($base), $base+$range, "next prime of $base is $base+$range" );
+ is( prev_prime($base+$range), $base, "prev prime of $base+$range is $base" );
+}
+
+is( next_prime(19608), 19609, "next prime of 19608 is 19609" );
+is( next_prime(19610), 19661, "next prime of 19610 is 19661" );
+is( next_prime(19660), 19661, "next prime of 19660 is 19661" );
+is( prev_prime(19662), 19661, "prev prime of 19662 is 19661" );
+is( prev_prime(19660), 19609, "prev prime of 19660 is 19609" );
+is( prev_prime(19610), 19609, "prev prime of 19610 is 19609" );
diff --git a/t/13-primecount.t b/t/13-primecount.t
new file mode 100644
index 0000000..a0711b7
--- /dev/null
+++ b/t/13-primecount.t
@@ -0,0 +1,52 @@
+#!/usr/bin/env perl
+use strict;
+use warnings;
+
+use Test::More;
+use Math::Prime::Util qw/prime_count prime_count_lower prime_count_upper prime_count_approx/;
+
+my $use64 = Math::Prime::Util::_maxbits > 32;
+my $extra = defined $ENV{RELEASE_TESTING} && $ENV{RELEASE_TESTING};
+
+plan tests => 10*3 + ($extra ? 10 : 7) + ($use64 ? 2*9 : 0);
+
+my %pivals32 = (
+ 1 => 0,
+ 10 => 4,
+ 100 => 25,
+ 1000 => 168,
+ 10000 => 1229,
+ 100000 => 9592,
+ 1000000 => 78498,
+ 10000000 => 664579,
+ 100000000 => 5761455,
+ 1000000000 => 50847534,
+);
+my %pivals64 = (
+ 10000000000 => 455052511,
+ 100000000000 => 4118054813,
+ 1000000000000 => 37607912018,
+ 10000000000000 => 346065536839,
+ 100000000000000 => 3204941750802,
+ 1000000000000000 => 29844570422669,
+ 10000000000000000 => 279238341033925,
+ 100000000000000000 => 2623557157654233,
+1000000000000000000 => 24739954287740860,
+);
+while (my($n, $pin) = each (%pivals32)) {
+ cmp_ok( prime_count_upper($n), '>=', $pin, "Pi($n) <= upper estimate" );
+ cmp_ok( prime_count_lower($n), '<=', $pin, "Pi($n) >= lower estimate" );
+ if ( ($n <= 2000000) || $extra ) {
+ is( prime_count($n), $pin, "Pi($n) = $pin" );
+ }
+ my $approx_range = abs($pin - prime_count_approx($n));
+ my $range_limit = 1100;
+ cmp_ok( $approx_range, '<=', $range_limit, "prime_count_approx($n) within $range_limit");
+}
+if ($use64) {
+ while (my($n, $pin) = each (%pivals64)) {
+ cmp_ok( prime_count_upper($n), '>=', $pin, "Pi($n) <= upper estimate" );
+ cmp_ok( prime_count_lower($n), '<=', $pin, "Pi($n) >= lower estimate" );
+ }
+}
+
diff --git a/t/14-nthprime.t b/t/14-nthprime.t
new file mode 100644
index 0000000..fcc1674
--- /dev/null
+++ b/t/14-nthprime.t
@@ -0,0 +1,71 @@
+#!/usr/bin/env perl
+use strict;
+use warnings;
+
+use Test::More;
+use Math::Prime::Util qw/nth_prime nth_prime_lower nth_prime_upper nth_prime_approx/;
+
+my $use64 = Math::Prime::Util::_maxbits > 32;
+my $extra = defined $ENV{RELEASE_TESTING} && $ENV{RELEASE_TESTING};
+
+plan tests => 7*2 + 9*3 + ($extra ? 9 : 7) + ($use64 ? 4*3 : 0);
+
+my %pivals32 = (
+ 1 => 0,
+ 10 => 4,
+ 100 => 25,
+ 1000 => 168,
+ 10000 => 1229,
+ 100000 => 9592,
+ 1000000 => 78498,
+);
+
+while (my($n, $pin) = each (%pivals32)) {
+ my $next = $pin+1;
+ cmp_ok( nth_prime($pin), '<=', $n, "nth_prime($pin) <= $n");
+ cmp_ok( nth_prime($next), '>=', $n, "nth_prime($next) >= $n");
+}
+
+my %nthprimes32 = (
+ 1 => 2,
+ 10 => 29,
+ 100 => 541,
+ 1000 => 7919,
+ 10000 => 104729,
+ 100000 => 1299709,
+ 1000000 => 15485863,
+ 10000000 => 179424673,
+ 100000000 => 2038074743,
+);
+my %nthprimes64 = (
+ 1000000000 => 22801763489,
+ 10000000000 => 252097800623,
+ 100000000000 => 2760727302517,
+ 1000000000000 => 29996224275833,
+ # TODO: find more
+);
+
+while (my($n, $nth) = each (%nthprimes32)) {
+ cmp_ok( nth_prime_upper($n), '>=', $nth, "nth_prime($n) <= upper estimate" );
+ cmp_ok( nth_prime_lower($n), '<=', $nth, "nth_prime($n) >= lower estimate" );
+
+ if ( ($n <= 2000000) || $extra ) {
+ is( nth_prime($n), $nth, "nth_prime($n) = $nth" );
+ }
+
+ my $approx = nth_prime_approx($n);
+ my $percent_limit = ($n >= 775) ? 1 : 2;
+ cmp_ok( abs($nth - $approx) / $nth, '<=', $percent_limit/100.0, "nth_prime_approx($n) = $approx within $percent_limit\% of $nth");
+}
+
+if ($use64) {
+ while (my($n, $nth) = each (%nthprimes64)) {
+ cmp_ok( nth_prime_upper($n), '>=', $nth, "nth_prime($n) <= upper estimate" );
+ cmp_ok( nth_prime_lower($n), '<=', $nth, "nth_prime($n) >= lower estimate" );
+
+ my $approx = nth_prime_approx($n);
+ my $percent_limit = 0.001;
+ cmp_ok( abs($nth - $approx) / $nth, '<=', $percent_limit/100.0, "nth_prime_approx($n) = $approx within $percent_limit\% of $nth");
+ }
+}
+
diff --git a/t/30-relations.t b/t/30-relations.t
new file mode 100644
index 0000000..d1b907a
--- /dev/null
+++ b/t/30-relations.t
@@ -0,0 +1,36 @@
+#!/usr/bin/env perl
+use strict;
+use warnings;
+
+use Test::More;
+use Math::Prime::Util qw/primes
+ nth_prime nth_prime_lower nth_prime_upper nth_prime_approx
+ prime_count prime_count_lower prime_count_upper prime_count_approx
+ next_prime prev_prime
+ /;
+
+
+my @trials = qw/1 2 3 4 5 6 7 17 57 89 102 1337 8573 84763 784357 1000001 2573622/;
+
+my $use64 = Math::Prime::Util::_maxbits > 32;
+my $extra = defined $ENV{RELEASE_TESTING} && $ENV{RELEASE_TESTING};
+
+plan tests => 5 * scalar @trials;
+
+my $last = 0;
+foreach my $n (@trials) {
+
+ is( prime_count($n), scalar @{primes($n)}, "Prime count and scalar primes agree for $n" );
+
+ is( prime_count($n) - prime_count($last),
+ scalar @{primes( $last+1, $n )},
+ "scalar primes($last+1,$n) = prime_count($n) - prime_count($last)" );
+
+ is( prime_count(nth_prime($n)), $n, "Pi(pn)) = n for $n");
+
+ is( nth_prime(prime_count($n)+1), next_prime($n), "p(Pi(n)+1) = next_prime(n) for $n" );
+
+ is( nth_prime(prime_count($n)), prev_prime($n+1), "p(Pi(n)) = prev_prime(n) for $n" );
+
+ $last = $n;
+}
diff --git a/t/90-release-perlcritic.t b/t/90-release-perlcritic.t
new file mode 100644
index 0000000..b98bd9f
--- /dev/null
+++ b/t/90-release-perlcritic.t
@@ -0,0 +1,24 @@
+#!/usr/bin/perl
+use strict;
+use warnings;
+use Test::More;
+
+BEGIN {
+ unless ($ENV{RELEASE_TESTING}) {
+ plan( skip_all => 'these tests are for release candidate testing' );
+ }
+}
+
+#---------------------------------------------------------------------
+
+
+eval { require Test::Perl::Critic; };
+plan skip_all => "Test::Perl::Critic required for testing PBP compliance" if $@;
+
+Test::Perl::Critic->import(
+ -verbose => 10,
+ -severity => 'gentle', # default
+ -force => 0, # default (allow ## no critic)
+ );
+
+all_critic_ok();
diff --git a/t/91-release-pod-syntax.t b/t/91-release-pod-syntax.t
new file mode 100644
index 0000000..2b0118e
--- /dev/null
+++ b/t/91-release-pod-syntax.t
@@ -0,0 +1,19 @@
+#!/usr/bin/perl
+use strict;
+use warnings;
+
+BEGIN {
+ unless ($ENV{RELEASE_TESTING}) {
+ require Test::More;
+ Test::More::plan(skip_all => 'these tests are for release candidate testing');
+ }
+}
+
+#---------------------------------------------------------------------
+
+
+use Test::More;
+eval "use Test::Pod 1.41";
+plan skip_all => "Test::Pod 1.41 required for testing POD" if $@;
+
+all_pod_files_ok();
diff --git a/t/92-release-pod-coverage.t b/t/92-release-pod-coverage.t
new file mode 100644
index 0000000..4590319
--- /dev/null
+++ b/t/92-release-pod-coverage.t
@@ -0,0 +1,25 @@
+#!/usr/bin/env perl
+use strict;
+use warnings;
+
+BEGIN {
+ unless ($ENV{RELEASE_TESTING}) {
+ require Test::More;
+ Test::More::plan(skip_all => 'these tests are for release candidate testing');
+ }
+}
+
+#---------------------------------------------------------------------
+
+
+use Test::More;
+eval "use Test::Pod::Coverage 1.08";
+plan skip_all => "Test::Pod::Coverage 1.08 required for testing POD coverage"
+ if $@;
+
+my @modules = Test::Pod::Coverage::all_modules();
+plan tests => scalar @modules;
+
+foreach my $m (@modules) {
+ pod_coverage_ok( $m, { also_private => [ qr/^(erat|erat_simple|segment|trial|sieve)_primes$/ ] } );
+}
diff --git a/util.c b/util.c
new file mode 100644
index 0000000..a775239
--- /dev/null
+++ b/util.c
@@ -0,0 +1,507 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <assert.h>
+#include <limits.h>
+#include <math.h>
+
+#include "util.h"
+#include "sieve.h"
+#include "bitarray.h"
+
+
+static const unsigned char byte_zeros[256] =
+ {8,7,7,6,7,6,6,5,7,6,6,5,6,5,5,4,7,6,6,5,6,5,5,4,6,5,5,4,5,4,4,3,
+ 7,6,6,5,6,5,5,4,6,5,5,4,5,4,4,3,6,5,5,4,5,4,4,3,5,4,4,3,4,3,3,2,
+ 7,6,6,5,6,5,5,4,6,5,5,4,5,4,4,3,6,5,5,4,5,4,4,3,5,4,4,3,4,3,3,2,
+ 6,5,5,4,5,4,4,3,5,4,4,3,4,3,3,2,5,4,4,3,4,3,3,2,4,3,3,2,3,2,2,1,
+ 7,6,6,5,6,5,5,4,6,5,5,4,5,4,4,3,6,5,5,4,5,4,4,3,5,4,4,3,4,3,3,2,
+ 6,5,5,4,5,4,4,3,5,4,4,3,4,3,3,2,5,4,4,3,4,3,3,2,4,3,3,2,3,2,2,1,
+ 6,5,5,4,5,4,4,3,5,4,4,3,4,3,3,2,5,4,4,3,4,3,3,2,4,3,3,2,3,2,2,1,
+ 5,4,4,3,4,3,3,2,4,3,3,2,3,2,2,1,4,3,3,2,3,2,2,1,3,2,2,1,2,1,1,0};
+static UV count_zero_bits(const unsigned char* m, UV nbytes)
+{
+ UV count = 0;
+ while (nbytes--)
+ count += byte_zeros[*m++];
+ return count;
+}
+
+
+static int _is_prime7(UV x)
+{
+ UV q, i;
+ i = 7;
+ while (1) { /* trial division, skipping multiples of 2/3/5 */
+ q = x/i; if (q<i) return 1; if (x==(q*i)) return 0; i += 4;
+ q = x/i; if (q<i) return 1; if (x==(q*i)) return 0; i += 2;
+ q = x/i; if (q<i) return 1; if (x==(q*i)) return 0; i += 4;
+ q = x/i; if (q<i) return 1; if (x==(q*i)) return 0; i += 2;
+ q = x/i; if (q<i) return 1; if (x==(q*i)) return 0; i += 4;
+ q = x/i; if (q<i) return 1; if (x==(q*i)) return 0; i += 6;
+ q = x/i; if (q<i) return 1; if (x==(q*i)) return 0; i += 2;
+ q = x/i; if (q<i) return 1; if (x==(q*i)) return 0; i += 6;
+ }
+ return 1;
+}
+
+
+/* Marked bits for each n, indicating if the number is prime */
+static const unsigned char prime_is_small[] =
+ {0xac,0x28,0x8a,0xa0,0x20,0x8a,0x20,0x28,0x88,0x82,0x08,0x02,0xa2,0x28,0x02,
+ 0x80,0x08,0x0a,0xa0,0x20,0x88,0x20,0x28,0x80,0xa2,0x00,0x08,0x80,0x28,0x82,
+ 0x02,0x08,0x82,0xa0,0x20,0x0a,0x20,0x00,0x88,0x22,0x00,0x08,0x02,0x28,0x82,
+ 0x80,0x20,0x88,0x20,0x20,0x02,0x02,0x28,0x80,0x82,0x08,0x02,0xa2,0x08,0x80,
+ 0x80,0x08,0x88,0x20,0x00,0x0a,0x00,0x20,0x08,0x20,0x08,0x0a,0x02,0x08,0x82,
+ 0x82,0x20,0x0a,0x80,0x00,0x8a,0x20,0x28,0x00,0x22,0x08,0x08,0x20,0x20,0x80,
+ 0x80,0x20,0x88,0x80,0x20,0x02,0x22,0x00,0x08,0x20,0x00,0x0a,0xa0,0x28,0x80,
+ 0x00,0x20,0x8a,0x00,0x20,0x8a,0x00,0x00,0x88,0x80,0x00,0x02,0x22,0x08,0x02};
+#define NPRIME_IS_SMALL (sizeof(prime_is_small)/sizeof(prime_is_small[0]))
+
+int is_prime(UV n)
+{
+ UV d, m;
+ unsigned char mtab;
+ const unsigned char* sieve;
+
+ if ( n < (NPRIME_IS_SMALL*8))
+ return ((prime_is_small[n/8] >> (n%8)) & 1);
+
+ d = n/30;
+ m = n - d*30;
+ mtab = masktab30[m]; /* Bitmask in mod30 wheel */
+
+ /* Return 0 if a multiple of 2, 3, or 5 */
+ if (mtab == 0)
+ return 0;
+
+ if (n <= get_prime_cache(0, &sieve))
+ return ((sieve[d] & mtab) == 0);
+
+ /* Trial division, mod 30 */
+ return _is_prime7(n);
+}
+
+
+static const unsigned char prime_next_small[] =
+ {2,2,3,5,5,7,7,11,11,11,11,13,13,17,17,17,17,19,19,23,23,23,23,
+ 29,29,29,29,29,29,31,31,37,37,37,37,37,37,41,41,41,41,43,43,47,
+ 47,47,47,53,53,53,53,53,53,59,59,59,59,59,59,61,61,67,67,67,67,67,67,71};
+#define NPRIME_NEXT_SMALL (sizeof(prime_next_small)/sizeof(prime_next_small[0]))
+
+UV next_trial_prime(UV n)
+{
+ UV d,m;
+
+ if (n < NPRIME_NEXT_SMALL)
+ return prime_next_small[n];
+
+ d = n/30;
+ m = n - d*30;
+ m = nextwheel30[m]; if (m == 1) d++;
+ while (!_is_prime7(d*30+m)) {
+ m = nextwheel30[m]; if (m == 1) d++;
+ }
+ return(d*30+m);
+}
+
+
+UV next_prime(UV n)
+{
+ UV d, m;
+ const unsigned char* sieve;
+ UV sieve_size;
+
+ if (n < NPRIME_NEXT_SMALL)
+ return prime_next_small[n];
+
+ sieve_size = get_prime_cache(0, &sieve);
+ if (n < sieve_size) {
+ START_DO_FOR_EACH_SIEVE_PRIME(sieve, n+1, sieve_size)
+ return p;
+ END_DO_FOR_EACH_SIEVE_PRIME;
+ /* Not found, so must be larger than the cache size */
+ n = sieve_size;
+ }
+
+ d = n/30;
+ m = n - d*30;
+ m = nextwheel30[m]; if (m == 1) d++;
+ while (!_is_prime7(d*30+m)) {
+ m = nextwheel30[m]; if (m == 1) d++;
+ }
+ return(d*30+m);
+}
+
+
+UV prev_prime(UV n)
+{
+ UV d, m;
+ const unsigned char* sieve;
+ UV sieve_size;
+
+ /* TODO: small prev prime */
+ if (n <= 7)
+ return (n <= 2) ? 0 : (n <= 3) ? 2 : (n <= 5) ? 3 : 5;
+
+ d = n/30;
+ m = n - d*30;
+
+ sieve_size = get_prime_cache(0, &sieve);
+ if (n < sieve_size) {
+ do {
+ m = prevwheel30[m]; if (m==29) { if (d == 0) return 0; d--; }
+ } while (sieve[d] & masktab30[m]);
+ } else {
+ do {
+ m = prevwheel30[m]; if (m==29) { if (d == 0) return 0; d--; }
+ } while (!_is_prime7(d*30+m));
+ }
+ return(d*30+m);
+}
+
+
+
+
+/*
+ * The pi(x) prime count functions. prime_count(x) gives an exact number,
+ * but requires determining all the primes up to x, so will be much slower.
+ *
+ * prime_count_lower(x) and prime_count_upper(x) give lower and upper limits,
+ * which will bound the exact value. These bounds should be fairly tight.
+ *
+ * pi_upper(x) - pi(x) pi_lower(x) - pi(x)
+ * < 10 for x < 5_371 < 10 for x < 9_437
+ * < 50 for x < 295_816 < 50 for x < 136_993
+ * < 100 for x < 1_761_655 < 100 for x < 909_911
+ * < 200 for x < 9_987_821 < 200 for x < 8_787_901
+ * < 400 for x < 34_762_891 < 400 for x < 30_332_723
+ * < 1000 for x < 372_748_528 < 1000 for x < 233_000_533
+ * < 5000 for x < 1_882_595_905 < 5000 for x < over 4300M
+ *
+ * The average of the upper and lower bounds is within 9 for all x < 15809, and
+ * within 50 for all x < 1_763_367.
+ *
+ * It is common to use the following Chebyshev inequality for x >= 17:
+ * 1*x/logx <-> 1.25506*x/logx
+ * but this gives terribly loose bounds.
+ *
+ * Rosser and Schoenfeld's bound for x >= 67 of
+ * x/(logx-1/2) <-> x/(logx-3/2)
+ * is much tighter. These bounds can be tightened even more.
+ *
+ * The formulas of Dusart for higher x are better yet. I recommend the paper
+ * by Burde for further information. Dusart's thesis is also a good resource.
+ *
+ * I have tweaked the bounds formulas for small (under 4000M) numbers so they
+ * are tighter. These bounds are verified via trial. The Dusart bounds
+ * (1.8 and 2.51) are used for larger numbers since those are proven.
+ *
+ */
+
+static const unsigned char prime_count_small[] =
+ {0,0,1,2,2,3,3,4,4,4,4,5,5,6,6,6,6,7,7,8,8,8,8,9,9,9,9,9,9,10,10,
+ 11,11,11,11,11,11,12,12,12,12,13,13,14,14,14,14,15,15,15,15,15,15,
+ 16,16,16,16,16,16,17,17,18,18,18,18,18,18,19};
+#define NPRIME_COUNT_SMALL (sizeof(prime_count_small)/sizeof(prime_count_small[0]))
+
+static const float F1 = 1.0;
+UV prime_count_lower(UV x)
+{
+ float fx, flogx;
+ float a = 1.80;
+
+ if (x < NPRIME_COUNT_SMALL)
+ return prime_count_small[x];
+
+ fx = (float)x;
+ flogx = logf(x);
+
+ if (x < 599)
+ return (UV) (fx / (flogx-0.7));
+
+ if (x < 2700) { a = 0.30; }
+ else if (x < 5500) { a = 0.90; }
+ else if (x < 19400) { a = 1.30; }
+ else if (x < 32299) { a = 1.60; }
+ else if (x < 176000) { a = 1.80; }
+ else if (x < 315000) { a = 2.10; }
+ else if (x < 1100000) { a = 2.20; }
+ else if (x < 4500000) { a = 2.31; }
+ else if (x <233000000) { a = 2.36; }
+ else if (x <240000000) { a = 2.32; }
+ else if (x <UVCONST(0xFFFFFFFF)) { a = 2.32; }
+
+ return (UV) ( (fx/flogx) * (F1 + F1/flogx + a/(flogx*flogx)) );
+}
+
+
+UV prime_count_upper(UV x)
+{
+ float fx, flogx;
+ float a = 2.51;
+
+ if (x < NPRIME_COUNT_SMALL)
+ return prime_count_small[x];
+
+ fx = (float)x;
+ flogx = logf(x);
+
+ /* This function is unduly complicated. */
+
+ if (x < 1621) return (UV) (fx / (flogx-1.048) + F1);
+ if (x < 5000) return (UV) (fx / (flogx-1.071) + F1);
+ if (x < 15900) return (UV) (fx / (flogx-1.098) + F1);
+
+ if (x < 24000) { a = 2.30; }
+ else if (x < 59000) { a = 2.48; }
+ else if (x < 350000) { a = 2.52; }
+ else if (x < 355991) { a = 2.54; }
+ else if (x < 356000) { a = 2.51; }
+ else if (x < 3550000) { a = 2.50; }
+ else if (x < 3560000) { a = 2.49; }
+ else if (x < 5000000) { a = 2.48; }
+ else if (x < 8000000) { a = 2.47; }
+ else if (x < 13000000) { a = 2.46; }
+ else if (x < 18000000) { a = 2.45; }
+ else if (x < 31000000) { a = 2.44; }
+ else if (x < 41000000) { a = 2.43; }
+ else if (x < 48000000) { a = 2.42; }
+ else if (x <119000000) { a = 2.41; }
+ else if (x <182000000) { a = 2.40; }
+ else if (x <192000000) { a = 2.395; }
+ else if (x <213000000) { a = 2.390; }
+ else if (x <271000000) { a = 2.385; }
+ else if (x <322000000) { a = 2.380; }
+ else if (x <400000000) { a = 2.375; }
+ else if (x <510000000) { a = 2.370; }
+ else if (x <682000000) { a = 2.367; }
+ else if (x <UVCONST(0xFFFFFFFF)) { a = 2.362; }
+
+ /*
+ * An alternate idea:
+ * float alog[23] = { 2.30,2.30,2.30,2.30,2.30,2.30,2.30 ,2.30,2.30,2.30,
+ * 2.47,2.49,2.53,2.50,2.49,2.49,2.456,2.44,2.40,2.370,
+ * 2.362,2.362,2.362,2.362};
+ * float clog[23] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
+ * 3, 1, 2, 1, 3, 2, 5, -6, 1, 1,
+ * 1, 1, 1, 1};
+ * if ((int)flogx < 23) {
+ * a = alog[(int)flogx];
+ * return ((UV) ( (fx/flogx) * (F1 + F1/flogx + a/(flogx*flogx)) ) + clog[(int)flogx] + 0.01);
+ * }
+ *
+ * Another thought is to use more terms in the Li(x) expansion along with
+ * a subtraction [Li(x) > Pi(x) for x < 10^316 or so, so for our 64-bit
+ * version we should be fine].
+ */
+
+ return (UV) ( (fx/flogx) * (F1 + F1/flogx + a/(flogx*flogx)) + F1 );
+}
+
+
+UV prime_count_approx(UV x)
+{
+ /* Placeholder for fancy algorithms, like Tomás Oliveira e Silva's:
+ * http://trac.sagemath.org/sage_trac/ticket/8135
+ * or an approximation to Li(x) plus a delta.
+ */
+ UV lower = prime_count_lower(x);
+ UV upper = prime_count_upper(x);
+ return ((lower + upper) / 2);
+}
+
+
+UV prime_count(UV n)
+{
+ const unsigned char* sieve;
+ static UV last_bytes = 0;
+ static UV last_count = 3;
+ UV s, bytes;
+ UV count = 3;
+
+ if (n < NPRIME_COUNT_SMALL)
+ return prime_count_small[n];
+
+ /* Get the cached sieve. */
+ if (get_prime_cache(n, &sieve) < n) {
+ croak("Couldn't generate sieve for prime_count");
+ return 0;
+ }
+
+ /* Simple:
+ * START_DO_FOR_EACH_SIEVE_PRIME(sieve, 7, n)
+ * count++;
+ * END_DO_FOR_EACH_SIEVE_PRIME;
+ */
+
+ bytes = n / 30;
+ s = 0;
+
+ /* Start from last word position if we can. This is a big speedup when
+ * calling prime_count many times with successively larger numbers. */
+ if (bytes >= last_bytes) {
+ s = last_bytes;
+ count = last_count;
+ }
+
+ count += count_zero_bits(sieve+s, bytes-s);
+
+ last_bytes = bytes;
+ last_count = count;
+
+ START_DO_FOR_EACH_SIEVE_PRIME(sieve, 30*bytes, n)
+ count++;
+ END_DO_FOR_EACH_SIEVE_PRIME;
+
+ return count;
+}
+
+
+
+static const unsigned short primes_small[] =
+ {0,2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,
+ 101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,
+ 193,197,199,211,223,227,229,233,239,241,251,257,263,269,271,277,281,283,
+ 293,307,311,313,317,331,337,347,349,353,359,367,373,379,383,389,397,401,
+ 409,419,421,431,433,439,443,449,457,461,463,467,479,487,491,499};
+#define NPRIMES_SMALL (sizeof(primes_small)/sizeof(primes_small[0]))
+
+/* The nth prime will be more than this number */
+UV nth_prime_lower(UV n)
+{
+ double fn = (double) n;
+ double flogn, flog2n, lower;
+
+ if (n < NPRIMES_SMALL)
+ return (n==0) ? 0 : primes_small[n]-1;
+
+ flogn = log(n);
+ flog2n = log(flogn);
+
+ /* Dusart 1999, for all n >= 2 */
+ lower = fn * (flogn + flog2n - 1.0 + ((flog2n-2.25)/flogn));
+
+ if (lower > (double)UV_MAX)
+ return 0;
+
+ return (UV) lower;
+}
+
+
+/* The nth prime will be less than this number */
+UV nth_prime_upper(UV n)
+{
+ double fn = (double) n;
+ double flogn, flog2n, upper;
+
+ if (n < NPRIMES_SMALL)
+ return primes_small[n]+1;
+
+ flogn = log(n);
+ flog2n = log(flogn);
+
+ if (n >= 39017)
+ upper = fn * ( flogn + flog2n - 0.9484 ); /* Dusart 1999 */
+ else if (n >= 27076)
+ upper = fn * (flogn + flog2n - 1.0 + ((flog2n-1.80)/flogn)); /*Dusart 1999*/
+ else if (n >= 7022)
+ upper = fn * ( flogn + 0.9385 * flog2n ); /* Robin 1983 */
+ else
+ upper = fn * ( flogn + flog2n );
+
+ if (upper > (double)UV_MAX)
+ return 0;
+
+ return (UV) ceil(upper);
+}
+
+
+UV nth_prime_approx(UV n)
+{
+ double fn, flogn, flog2n, order, approx;
+
+ if (n < NPRIMES_SMALL)
+ return primes_small[n];
+
+ /* This isn't too bad:
+ * return ((nth_prime_lower(n) + nth_prime_upper(n)) / 2);
+ */
+
+ fn = (double) n;
+ flogn = log(n);
+ flog2n = log(flogn);
+
+ /* Cipolla 1902:
+ * m=0 fn * ( flogn + flog2n - 1 );
+ * m=1 + ((flog2n - 2)/flogn) );
+ * m=2 - (((flog2n*flog2n) - 6*flog2n + 11) / (2*flogn*flogn))
+ * + O((flog2n/flogn)^3)
+ */
+ approx = fn * ( flogn + flog2n - 1
+ + ((flog2n - 2)/flogn)
+ - (((flog2n*flog2n) - 6*flog2n + 11) / (2*flogn*flogn))
+ );
+
+ /* Apply a correction to help keep small inputs close. */
+ order = flog2n/flogn;
+ order = order*order*order * fn;
+
+ if (n < 259) approx += 10.4 * order;
+ else if (n < 775) approx += 7.52 * order;
+ else if (n < 1271) approx += 5.6 * order;
+ else if (n < 2000) approx += 5.2 * order;
+ else if (n < 4000) approx += 4.3 * order;
+ else if (n < 12000) approx += 3.0 * order;
+ else if (n <150000) approx += 2.1 * order;
+
+ return (UV) rint(approx);
+}
+
+
+UV nth_prime(UV n)
+{
+ const unsigned char* sieve;
+ UV upper_limit, start, count, s, bytes_left;
+
+ if (n < NPRIMES_SMALL)
+ return primes_small[n];
+
+if (n == 100000) return 1299709;
+if (n == 1000000) return 15485863;
+if (n == 10000000) return 179424673;
+if (n == 100000000) return 2038074743;
+if (n == 1000000000) return 22801763489;
+if (n == 10000000000) return 252097800623;
+
+ upper_limit = nth_prime_upper(n);
+ if (upper_limit == 0) {
+ croak("nth_prime(%lu) would overflow", (unsigned long)n);
+ return 0;
+ }
+ /* The nth prime is guaranteed to be within this range */
+ if (get_prime_cache(upper_limit, &sieve) < upper_limit) {
+ croak("Couldn't generate sieve for nth(%lu) [sieve to %lu]", (unsigned long)n, (unsigned long)upper_limit);
+ return 0;
+ }
+
+ count = 3;
+ start = 7;
+ s = 0;
+ bytes_left = (n-count) / 8;
+ while ( bytes_left > 0 ) {
+ /* There is at minimum one byte we can count (and probably many more) */
+ count += count_zero_bits(sieve+s, bytes_left);
+ assert(count <= n);
+ s += bytes_left;
+ bytes_left = (n-count) / 8;
+ }
+ if (s > 0)
+ start = s * 30;
+
+ START_DO_FOR_EACH_SIEVE_PRIME(sieve, start, upper_limit)
+ if (++count == n) return p;
+ END_DO_FOR_EACH_SIEVE_PRIME;
+ croak("nth_prime failed for %lu, not found in range %lu - %lu", (unsigned long)n, (unsigned long) start, (unsigned long)upper_limit);
+ return 0;
+}
diff --git a/util.h b/util.h
new file mode 100644
index 0000000..5c81e2c
--- /dev/null
+++ b/util.h
@@ -0,0 +1,22 @@
+#ifndef MPU_UTIL_H
+#define MPU_UTIL_H
+
+#include "EXTERN.h"
+#include "perl.h"
+
+extern int is_prime(UV x);
+extern UV next_trial_prime(UV x);
+extern UV next_prime(UV x);
+extern UV prev_prime(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 prime_count(UV x);
+
+extern UV nth_prime_lower(UV n);
+extern UV nth_prime_upper(UV x);
+extern UV nth_prime_approx(UV x);
+extern UV nth_prime(UV x);
+
+#endif
--
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