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

_______________________________________________
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