[libmath-prime-util-perl] 04/04: Updates for next release
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:44:14 UTC 2015
This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.02
in repository libmath-prime-util-perl.
commit 40b1a2a86271801d2803f23b431e9d611045bd79
Author: Dana Jacobsen <dana at acm.org>
Date: Tue Jun 5 05:59:44 2012 -0600
Updates for next release
---
Changes | 1 +
factor.c | 20 ++++---
lib/Math/Prime/Util.pm | 82 ++++++++++++++++++++++++----
sieve.c | 5 ++
sieve.h | 8 +--
t/14-nthprime.t | 8 ++-
util.c | 141 ++++++++++++++++++++++++++++++++++++++-----------
util.h | 4 ++
8 files changed, 214 insertions(+), 55 deletions(-)
diff --git a/Changes b/Changes
index 3b02d1d..78d7472 100644
--- a/Changes
+++ b/Changes
@@ -11,6 +11,7 @@ Revision history for Perl extension Math::Prime::Util.
Not memory intensive any more, and faster for large inputs. The time
growth is slightly over linear however, so expect to wait a long
time for 10^12 or more.
+ - nth_prime also transitioned to segmented sieve.
0.01 4 June 2012
- Initial release
diff --git a/factor.c b/factor.c
index 9ed9c11..7fa1474 100644
--- a/factor.c
+++ b/factor.c
@@ -23,8 +23,9 @@
int trial_factor(UV n, UV *factors, UV maxtrial)
{
- UV f, d, m;
- UV limit;
+ static const wheeladvance[30] =
+ {0,6,0,0,0,0,0,4,0,0,0,2,0,4,0,0,0,2,0,4,0,0,0,6,0,0,0,0,0,2};
+ UV f, m, limit;
int nfactors = 0;
if (maxtrial == 0) maxtrial = UV_MAX;
@@ -58,7 +59,6 @@ int trial_factor(UV n, UV *factors, UV maxtrial)
/* wheel 30 */
f = 11;
- d = 0;
m = 11;
while (f <= limit) {
if ( (n%f) == 0 ) {
@@ -70,8 +70,8 @@ int trial_factor(UV n, UV *factors, UV maxtrial)
newlimit = sqrt(n);
if (newlimit < limit) limit = newlimit;
}
- m = nextwheel30[m]; if (m == 1) d++;
- f = d*30 + m;
+ f += wheeladvance[m];
+ m = nextwheel30[m];
}
if (n != 1)
factors[nfactors++] = n;
@@ -91,15 +91,19 @@ static UV gcd_ui(UV x, UV y) {
}
/* n^power + a mod m */
+/* This has serious overflow issues, making the programs that use it dubious */
static UV powmod(UV n, UV power, UV add, UV m) {
UV t = 1;
+ n = n % m;
+ /* t and n will always be < m from now on */
while (power) {
if (power & 1)
- t = ((t % m) * (n % m)) % m;
- n = ((n % m) * (n % m)) % m;
+ t = (t * n) % m;
+ n = (n * n) % m;
power >>= 1;
}
- return (t+add) % m;
+ /* (t+a) % m, noting t is always < m */
+ return ( ((m-t) > add) ? (t+add) : (t+add-m) );
}
/* Knuth volume 2, algorithm C.
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 64747a1..1cedbf8 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -20,6 +20,7 @@ our @EXPORT_OK = qw(
nth_prime nth_prime_lower nth_prime_upper nth_prime_approx
factor
);
+our %EXPORT_TAGS = (all => [ @EXPORT_OK ]);
BEGIN {
eval {
@@ -136,7 +137,10 @@ Version 0.02
=head1 SYNOPSIS
- use Math::Prime::Util qw/primes/;
+ # Normally you would just import the functions you are using.
+ # Nothing is exported by default.
+ use Math::Prime::Util ':all';
+
# Get a big array reference of many primes
my $aref = primes( 100_000_000 );
@@ -144,6 +148,55 @@ Version 0.02
# All the primes between 5k and 10k inclusive
my $aref = primes( 5_000, 10_000 );
+ # If you want them in an array instead
+ my @primes = @{primes( 500 )};
+
+
+ # is_prime returns 0 for composite, 1 for prime
+ say "$n is prime" if is_prime($n);
+
+
+ # step to the next prime
+ $n = next_prime($n);
+
+ # step back (returns 0 if given input less than 2)
+ $n = prev_prime($n);
+
+
+ # Return Pi(n) -- the number of primes E<lt>= n.
+ my $number_of_primes = prime_count( 1_000_000 );
+
+ # Quickly return an approximation to Pi(n)
+ my $approx_number_of_primes = prime_count_approx( 10**17 );
+
+ # Lower and upper bounds. lower <= Pi(n) <= upper for all n
+ die unless prime_count_lower($n) <= prime_count($n);
+ die unless prime_count_upper($n) >= prime_count($n);
+
+
+ # Return p_n, the nth prime
+ say "The ten thousandth prime is ", nth_prime(10_000);
+
+ # Return a quick approximation to the nth prime
+ say "The one trillionth prime is ~ ", nth_prime_approx(10**12);
+
+ # Lower and upper bounds. lower <= nth_prime(n) <= upper for all n
+ die unless nth_prime_lower($n) <= nth_prime($n);
+ die unless nth_prime_upper($n) >= nth_prime($n);
+
+
+ # Get the prime factors of a number
+ @prime_factors = factor( $n );
+
+
+ # Precalculate a sieve, possibly speeding up later work.
+ prime_precalc( 1_000_000_000 );
+
+ # Free any memory used by the module.
+ prime_free;
+
+ # Alternate way to free. When this leaves scope, memory is freed.
+ my $mf = Math::Prime::Util::MemFree->new;
=head1 DESCRIPTION
@@ -266,7 +319,10 @@ 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.
+large inputs. A segmented sieve is used for large inputs, so it is memory
+efficient. On my machine it will return the 203,280,221st prime (the largest
+that fits in 32-bits) in 2.5 seconds. The 10^9th prime takes 15 seconds to
+find, while the 10^10th prime takes nearly four minutes.
=head2 nth_prime_upper
@@ -319,7 +375,7 @@ might be a better choice for complicated uses.
=head2 Math::Prime::Util::MemFree->new
- my $mf = Math::Prime::Util::MemFree->new
+ my $mf = Math::Prime::Util::MemFree->new;
# perform operations. When $mf goes out of scope, memory will be recovered.
This is a more robust way of making sure any cached memory is freed, as it
@@ -390,10 +446,6 @@ typically much faster for inputs in the 19+ digit range.
=head1 LIMITATIONS
-The function C<nth_prime> has 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.
@@ -414,6 +466,7 @@ Pi(10^10) = 455,052,511.
6.6 primegen (optimized Sieve of Atkin)
11.2 Tomás Oliveira e Silva's segmented sieve v1 (May 2003)
+ 5.9 My SoE called in segments
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
@@ -423,11 +476,22 @@ Pi(10^10) = 455,052,511.
Perl modules, counting the primes to C<800_000_000> (800 million), in seconds:
+ Time (s) Module Version
+ -------- -------------------------- ------
+ 0.4 Math::Prime::Util 0.02
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
+ [hours] Math::Primality 0.04
+
+
+I have not done extensive timing on factoring. The difference in speed for
+most inputs between Math::Factor::XS and Math::Prime::Util is small.
+M::F::XS is a bit faster than the current implementation for most ranges,
+except for large inputs with large factors, where using SQUFOF is a big
+advantage (e.g. C<204518747 * 16476429743>, which is ~50x faster with this
+module).
=head1 AUTHORS
@@ -449,7 +513,7 @@ will.
The SQUFOF implementation being used is my modifications to Ben Buhrow's
modifications to Bob Silverman's code. I may experiment with some other
-implementations (Ben Buhrows and Jason Papadopoulos both have published their
+implementations (Ben Buhrows and Jason Papadopoulos both have published
excellent versions in the public domain).
diff --git a/sieve.c b/sieve.c
index 8af0f70..adb8ead 100644
--- a/sieve.c
+++ b/sieve.c
@@ -7,6 +7,7 @@
#include "sieve.h"
#include "bitarray.h"
+#include "util.h" /* For freeing the segment cache */
static unsigned char* prime_cache_sieve = 0;
@@ -61,6 +62,8 @@ void prime_precalc(UV n)
}
get_prime_cache(n, 0); /* Sieve to n */
+
+ /* TODO: should we prealloc the segment here? */
}
void prime_free(void)
@@ -70,6 +73,8 @@ void prime_free(void)
prime_cache_sieve = 0;
prime_cache_size = 0;
+ free_prime_segment();
+
prime_precalc(0);
}
diff --git a/sieve.h b/sieve.h
index 5c08bc4..28044d4 100644
--- a/sieve.h
+++ b/sieve.h
@@ -16,18 +16,18 @@ 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] = {
+static const 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] = {
+static const 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] = {
+static const 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] = {
+static const 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 };
diff --git a/t/14-nthprime.t b/t/14-nthprime.t
index 3b5e48a..1164387 100644
--- a/t/14-nthprime.t
+++ b/t/14-nthprime.t
@@ -8,7 +8,7 @@ use Math::Prime::Util qw/nth_prime nth_prime_lower nth_prime_upper nth_prime_app
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);
+plan tests => 7*2 + 9*3 + ($extra ? 9 : 7) + ($use64 ? 9*3 : 0);
my %pivals32 = (
1 => 0,
@@ -42,7 +42,11 @@ my %nthprimes64 = (
10000000000 => 252097800623,
100000000000 => 2760727302517,
1000000000000 => 29996224275833,
- # TODO: find more
+ 10000000000000 => 323780508946331,
+ 100000000000000 => 3475385758524527,
+ 1000000000000000 => 37124508045065437,
+ 10000000000000000 => 394906913903735329,
+ 100000000000000000 => 4185296581467695669,
);
while (my($n, $nth) = each (%nthprimes32)) {
diff --git a/util.c b/util.c
index e1cea2b..169f38b 100644
--- a/util.c
+++ b/util.c
@@ -9,6 +9,23 @@
#include "sieve.h"
#include "bitarray.h"
+/*
+ * I'm undecided as to whether we want this, or just let the functions alloc
+ * and free it per call.
+ */
+static unsigned char* prime_segment = 0;
+unsigned char* get_prime_segment(void) {
+ if (prime_segment == 0)
+ prime_segment = (unsigned char*) malloc( SEGMENT_CHUNK_SIZE );
+ if (prime_segment == 0)
+ croak("Could not allocate %"UVuf" bytes for segment sieve", SEGMENT_CHUNK_SIZE);
+ return prime_segment;
+}
+void free_prime_segment(void) {
+ if (prime_segment != 0)
+ free(prime_segment);
+ prime_segment = 0;
+}
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,
@@ -349,17 +366,15 @@ UV prime_count(UV n)
} else {
/* We don't have enough primes. Repeatedly segment sieve */
- UV const segment_size = 262144;
+ UV const segment_size = SEGMENT_CHUNK_SIZE;
unsigned char* segment;
- /* TODO: we really shouldn't need this */
+ /* Get this over with once */
prime_precalc( sqrt(n) + 2 );
- segment = (unsigned char*) malloc( segment_size );
- if (segment == 0) {
- croak("Could not allocate %"UVuf" bytes for segment sieve", segment_size);
+ segment = get_prime_segment();
+ if (segment == 0)
return 0;
- }
for (s = 0; s <= bytes; s += segment_size) {
/* We want to sieve one extra byte, to handle the last fragment */
@@ -383,8 +398,6 @@ UV prime_count(UV n)
count++;
END_DO_FOR_EACH_SIEVE_PRIME;
- free(segment);
-
}
return count;
@@ -492,42 +505,106 @@ UV nth_prime_approx(UV n)
}
+/*
+ * Given a sieve of size nbytes, walk it counting zeros (primes) until:
+ *
+ * (1) we counted them all: return the count, which will be less than maxcount.
+ *
+ * (2) we hit maxcount: set position to the index of the maxcount'th prime
+ * and return count (which will be equal to maxcount).
+ */
+static UV count_segment(const unsigned char* sieve, UV nbytes, UV maxcount, UV* pos)
+{
+ UV count = 0;
+ UV bytes_left;
+ UV byte = 0;
+
+ assert(sieve != 0);
+ assert(pos != 0);
+ *pos = 0;
+ if ( (nbytes == 0) || (maxcount == 0) )
+ return 0;
+
+ while ( (byte < nbytes) && (count < maxcount) )
+ count += byte_zeros[sieve[byte++]];
+
+ if (count >= maxcount) { /* One too far -- back up */
+ count -= byte_zeros[sieve[--byte]];
+ }
+
+ assert(count < maxcount);
+
+ if (byte == nbytes)
+ return count;
+
+ /* The result is somewhere in the next byte */
+ START_DO_FOR_EACH_SIEVE_PRIME(sieve, byte*30+1, nbytes*30-1)
+ if (++count == maxcount) { *pos = p; return count; }
+ END_DO_FOR_EACH_SIEVE_PRIME;
+
+ croak("count_segment failed");
+ return 0;
+}
+
UV nth_prime(UV n)
{
- const unsigned char* sieve;
- UV upper_limit, start, count, s, bytes_left;
+ const unsigned char* cache_sieve;
+ unsigned char* segment;
+ UV upper_limit, segbase, segment_size;
+ UV p = 0;
+ UV target = n-3;
+ UV count = 0;
+ /* If very small, return the table entry */
if (n < NPRIMES_SMALL)
return primes_small[n];
+ /* Determine a bound on the nth prime. We know it comes before this. */
upper_limit = nth_prime_upper(n);
if (upper_limit == 0) {
croak("nth_prime(%"UVuf") would overflow", 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(%"UVuf") [sieve to %"UVuf"]", n, upper_limit);
+
+ /* Get the primary cache, and ensure it is at least this large. If the
+ * input is small enough, get a sieve covering the range. Otherwise, we'll
+ * walk segments. Make sure we have enough primes in the cache so the
+ * segmented siever won't have to keep resieving.
+ */
+ if (upper_limit <= (1*1024*1024*30))
+ segment_size = get_prime_cache(upper_limit, &cache_sieve) / 30;
+ else
+ segment_size = get_prime_cache(sqrt(upper_limit), &cache_sieve) / 30;
+
+ /* Count up everything in the cached sieve. */
+ count += count_segment(cache_sieve, segment_size, target, &p);
+ if (count == target)
+ return p;
+
+ /* Start segment sieving. Get memory to sieve into. */
+ segbase = segment_size;
+ segment_size = SEGMENT_CHUNK_SIZE;
+ segment = get_prime_segment();
+ if (segment == 0)
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;
+ while (count < target) {
+ /* Limit the segment size if we know the answer comes earlier */
+ if ( (30*(segbase+segment_size)+29) > upper_limit )
+ segment_size = (upper_limit - segbase*30 + 30) / 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 %"UVuf", not found in range %"UVuf" - %"UVuf, n, start, upper_limit);
- return 0;
+ /* Do the actual sieving in the range */
+ if (sieve_segment(segment, segbase, segbase + segment_size-1) == 0) {
+ croak("Could not segment sieve from %"UVuf" to %"UVuf, 30*segbase+1, 30*(segbase+segment_size)+29);
+ break;
+ }
+
+ /* Count up everything in this segment */
+ count += count_segment(segment, segment_size, target-count, &p);
+
+ if (count < target)
+ segbase += segment_size;
+ }
+ assert(count == target);
+ return ( (segbase*30) + p );
}
diff --git a/util.h b/util.h
index 5c81e2c..35659d0 100644
--- a/util.h
+++ b/util.h
@@ -19,4 +19,8 @@ extern UV nth_prime_upper(UV x);
extern UV nth_prime_approx(UV x);
extern UV nth_prime(UV x);
+#define SEGMENT_CHUNK_SIZE 262144
+extern unsigned char* get_prime_segment(void);
+extern void free_prime_segment(void);
+
#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