[libmath-prime-util-perl] 06/06: Changes for v0.04
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:44:21 UTC 2015
This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.04
in repository libmath-prime-util-perl.
commit 18c272dd6534ee68ebdc07173d0e209476e161d2
Author: Dana Jacobsen <dana at acm.org>
Date: Fri Jun 8 13:58:10 2012 -0600
Changes for v0.04
---
Changes | 1 +
MANIFEST | 7 +++-
TODO | 4 +-
XS.xs | 4 +-
factor.c | 60 +++++++++++++++++++-----------
lib/Math/Prime/Util.pm | 40 +++++++++++++-------
sieve.c | 3 +-
sieve.h | 2 +-
t/02-can.t | 6 ++-
t/03-init.t | 18 ++++-----
t/16-randomprime.t | 99 ++++++++++++++++++++++++++++++++++++++++++++++++++
util.c | 3 +-
12 files changed, 189 insertions(+), 58 deletions(-)
diff --git a/Changes b/Changes
index 64a99b7..9277d67 100644
--- a/Changes
+++ b/Changes
@@ -7,6 +7,7 @@ Revision history for Perl extension Math::Prime::Util.
is_prime run up to 10x faster (which impacts next_prime, factoring, etc.)
- replace all assert with croak indicating an internal error.
- Add random_prime and random_ndigit_prime
+ - renamed prime_free to prime_memfree.
0.03 6 June 2012
- Speed up factoring.
diff --git a/MANIFEST b/MANIFEST
index f12eb59..80c5c05 100644
--- a/MANIFEST
+++ b/MANIFEST
@@ -15,10 +15,12 @@ sieve.c
util.h
util.c
examples/bench-nthprime.pl
-examples/bench-factor-extra.t
+examples/bench-factor.pl
+examples/bench-factor-extra.pl
examples/bench-factor-semiprime.pl
examples/bench-is-prime.pl
-examples/bench-miller-rabin.t
+examples/bench-miller-rabin.pl
+examples/bench-random-prime.pl
examples/test-factoring.pl
t/01-load.t
t/02-can.t
@@ -29,6 +31,7 @@ t/12-nextprime.t
t/13-primecount.t
t/14-nthprime.t
t/15-probprime.t
+t/16-randomprime.t
t/30-relations.t
t/50-factoring.t
t/90-release-perlcritic.t
diff --git a/TODO b/TODO
index e3aa794..c0c7240 100644
--- a/TODO
+++ b/TODO
@@ -17,8 +17,6 @@
- Faster SQUFOF
-- prime count bounds are <=, nth_prime bounds are <. Pick one.
-
-- test file for random_prime
+- nth_prime bounds are < in the code, but <= in the documentation.
- speed up random_prime for large numbers
diff --git a/XS.xs b/XS.xs
index 1d23c6e..0188258 100644
--- a/XS.xs
+++ b/XS.xs
@@ -18,7 +18,7 @@ void
prime_precalc(IN UV n)
void
-prime_free()
+prime_memfree()
UV
prime_count(IN UV n)
@@ -148,7 +148,7 @@ segment_primes(IN UV low, IN UV high, IN UV segment_size = 65536UL)
UV seghigh = (seghigh_d == high_d) ? high : (seghigh_d*30+29);
UV segbase = low_d * 30;
/* printf(" startd = %"UVuf" endd = %"UVuf"\n", startd, endd); */
-
+
MPUassert( seghigh_d >= low_d, "segment_primes highd < lowd");
MPUassert( range_d <= segment_size, "segment_primes range > segment size");
diff --git a/factor.c b/factor.c
index fef8b55..b02ff1a 100644
--- a/factor.c
+++ b/factor.c
@@ -1,8 +1,7 @@
-//#include <stdio.h>
-//#include <stdlib.h>
-//#include <string.h>
-//#include <limits.h>
-//#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
#include "ptypes.h"
#include "factor.h"
@@ -101,25 +100,44 @@ int trial_factor(UV n, UV *factors, UV maxtrial)
/* UV is the largest integral type available (that we know of). */
- /* if n is smaller than this, you can multiply without overflow */
- #define HALF_WORD (UVCONST(1) << (BITS_PER_WORD/2))
-
- static UV _mulmod(UV a, UV b, UV m) {
- UV r = 0;
- while (b > 0) {
- if (b & 1) {
- if (r == 0) {
- r = a;
- } else {
- r = m - r;
- r = (a >= r) ? a-r : m-r+a;
+ #if defined(__GNUC__) && (defined(__i386__) || defined(__x86_64__))
+ /* Inline assembly -- basically as fast as a regular (a*b)%m */
+ /* Arguably we don't need all the ifs to avoid it */
+ static UV _mulmod(UV a, UV b, UV c) {
+ UV d; /* to hold the result of a*b mod c */
+ /* calculates a*b mod c, stores result in d */
+ asm ("mov %1, %%rax;" /* put a into rax */
+ "mul %2;" /* mul a*b -> rdx:rax */
+ "div %3;" /* (a*b)/c -> quot in rax remainder in rdx */
+ "mov %%rdx, %0;" /* store result in d */
+ :"=r"(d) /* output */
+ :"r"(a), "r"(b), "r"(c) /* input */
+ :"%rax", "%rdx" /* clobbered registers */
+ );
+ return d;
+ }
+ #else
+ /* Do it by hand */
+ static UV _mulmod(UV a, UV b, UV m) {
+ UV r = 0;
+ while (b > 0) {
+ if (b & 1) {
+ if (r == 0) {
+ r = a;
+ } else {
+ r = m - r;
+ r = (a >= r) ? a-r : m-r+a;
+ }
}
+ a = (a > (m-a)) ? (a-m)+a : a+a;
+ b >>= 1;
}
- a = (a > (m-a)) ? (a-m)+a : a+a;
- b >>= 1;
+ return r;
}
- return r;
- }
+ #endif
+
+ /* if n is smaller than this, you can multiply without overflow */
+ #define HALF_WORD (UVCONST(1) << (BITS_PER_WORD/2))
#define addmod(n,a,m) ((((m)-(n)) > (a)) ? ((n)+(a)) : ((n)+(a)-(m)))
#define mulmod(a,b,m) (((a)|(b)) < HALF_WORD) ? ((a)*(b))%(m) : _mulmod(a,b,m)
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 01dfb85..217f8a4 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -12,7 +12,7 @@ BEGIN {
# use parent qw( Exporter );
use base qw( Exporter );
our @EXPORT_OK = qw(
- prime_precalc prime_free
+ prime_precalc prime_memfree
is_prime is_prob_prime miller_rabin
primes
next_prime prev_prime
@@ -51,6 +51,8 @@ sub primes {
if ( (!defined $low) || (!defined $high) ||
($low =~ tr/0123456789//c) || ($high =~ tr/0123456789//c)
) {
+ $low = 'undef' unless defined $low;
+ $high = 'undef' unless defined $high;
croak "Parameters [ $low $high ] must be positive integers";
}
@@ -103,6 +105,8 @@ sub random_prime {
if ( (!defined $low) || (!defined $high) ||
($low =~ tr/0123456789//c) || ($high =~ tr/0123456789//c)
) {
+ $low = 'undef' unless defined $low;
+ $high = 'undef' unless defined $high;
croak "Parameters [ $low $high ] must be positive integers";
}
croak "Parameters [ $low $high ] not in range 0-$_maxparam" unless $low <= $_maxparam && $high <= $_maxparam;
@@ -131,12 +135,12 @@ sub random_prime {
if ($range <= 4294967295) {
do {
$prime = $low + int(rand($range));
- } while ( !($prime%2) && !($prime%3) && !is_prime($prime) );
+ } while ( !($prime%2) || !($prime%3) || !is_prime($prime) );
} else {
do {
my $rand = ( (int(rand(4294967295)) << 32) + int(rand(4294967295)) ) % $range;
$prime = $low + $rand;
- } while ( !($prime%2) && !($prime%3) && !is_prime($prime) );
+ } while ( !($prime%2) || !($prime%3) || !is_prime($prime) );
}
}
return $prime;
@@ -167,7 +171,7 @@ sub new {
sub DESTROY {
my $self = shift;
confess "instances count mismatch" unless $memfree_instances > 0;
- Math::Prime::Util::prime_free if --$memfree_instances == 0;
+ Math::Prime::Util::prime_memfree if --$memfree_instances == 0;
}
package Math::Prime::Util;
@@ -212,11 +216,14 @@ Version 0.04
my @primes = @{primes( 500 )};
- # is_prime returns 0 for composite, 1 for prime
+ # is_prime returns 0 for composite, 2 for prime
say "$n is prime" if is_prime($n);
+ # is_prob_prime returns 0 for composite, 2 for prime, and 1 for maybe prime
+ say "$n is ", qw(composite maybe_prime? prime)[is_prob_prime($n)];
- # step to the next prime
+
+ # step to the next prime (returns 0 if the next one is more than ~0)
$n = next_prime($n);
# step back (returns 0 if given input less than 2)
@@ -253,7 +260,7 @@ Version 0.04
prime_precalc( 1_000_000_000 );
# Free any memory used by the module.
- prime_free;
+ prime_memfree;
# Alternate way to free. When this leaves scope, memory is freed.
my $mf = Math::Prime::Util::MemFree->new;
@@ -477,9 +484,12 @@ The current algorithm does a random index selection for small numbers, which
is deterministic. For larger numbers, this can be very slow, so the
obvious Monte Carlo method is used, where random numbers in the range are
selected until one is prime. That also gets slow as the number of digits
-increases, but not something that impacts us in 64-bit. A GMP version would
-likely need to consider using C<next_prime>/C<prev_prime> -- giving up a
-small bit of uniformity for reasonable performance.
+increases, but not something that impacts us in 64-bit.
+
+If you want cryptographically secure primes, I suggest looking at
+L<Crypt::Primes> or something similar. The current L<Math::Prime::Util>
+module does not use strong randomness, and its primes are ridiculously small
+by cryptographic standards.
=head2 random_ndigit_prime
@@ -506,9 +516,9 @@ current implementation this will calculate a sieve for all numbers up to the
specified number.
-=head2 prime_free
+=head2 prime_memfree
- prime_free;
+ prime_memfree;
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
@@ -541,8 +551,9 @@ always result in the input value, though those are the only cases where
the returned factors are not prime.
The current algorithm is to use trial division for small numbers, while large
-numbers go through a sequence of small trials, SQUFOF, and Pollard's Rho.
-This process is repeated for each non-prime factor.
+numbers go through a sequence of small trials, SQUFOF, Pollard's Rho, and
+finally trial division for any survivors. This process is repeated for
+each non-prime factor.
=head2 trial_factor
@@ -626,6 +637,7 @@ 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
+ 2.2 yafu 1.31
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)
diff --git a/sieve.c b/sieve.c
index e5003d3..6ac1af6 100644
--- a/sieve.c
+++ b/sieve.c
@@ -1,7 +1,6 @@
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
-#include <limits.h>
#include <math.h>
#include "sieve.h"
@@ -66,7 +65,7 @@ void prime_precalc(UV n)
/* TODO: should we prealloc the segment here? */
}
-void prime_free(void)
+void prime_memfree(void)
{
if (prime_cache_sieve != 0)
free(prime_cache_sieve);
diff --git a/sieve.h b/sieve.h
index f381639..e693b5b 100644
--- a/sieve.h
+++ b/sieve.h
@@ -8,7 +8,7 @@ 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 void prime_memfree(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);
diff --git a/t/02-can.t b/t/02-can.t
index 34748ff..f3e0aea 100644
--- a/t/02-can.t
+++ b/t/02-can.t
@@ -6,11 +6,13 @@ use Math::Prime::Util;
use Test::More tests => 1;
my @functions = qw(
- prime_precalc prime_free
- is_prime
+ prime_precalc prime_memfree
+ is_prime is_prob_prime miller_rabin
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
+ random_prime random_ndigit_prime
+ factor
);
can_ok( 'Math::Prime::Util', @functions);
diff --git a/t/03-init.t b/t/03-init.t
index feb3940..079b07a 100644
--- a/t/03-init.t
+++ b/t/03-init.t
@@ -1,7 +1,7 @@
#!/usr/bin/env perl
use strict;
use warnings;
-use Math::Prime::Util qw/prime_precalc prime_free/;
+use Math::Prime::Util qw/prime_precalc prime_memfree/;
use Test::More tests => 3 + 3 + 3 + 6;
@@ -19,9 +19,9 @@ prime_precalc($bigsize);
cmp_ok( Math::Prime::Util::_get_prime_cache_size, '>', $init_size, "Internal space grew after large precalc" );
-prime_free;
+prime_memfree;
-is( Math::Prime::Util::_get_prime_cache_size, $init_size, "Internal space went back to original size after free" );
+is( Math::Prime::Util::_get_prime_cache_size, $init_size, "Internal space went back to original size after memfree" );
# Now do the object way.
@@ -48,12 +48,12 @@ is( Math::Prime::Util::_get_prime_cache_size, $init_size, "Memory released after
is( Math::Prime::Util::_get_prime_cache_size, $init_size, "Memory released after last MemFree object goes out of scope");
# Show how an eval death can leak
-eval { prime_precalc($bigsize); 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, "Memory freed after eval");
+eval { prime_precalc($bigsize); cmp_ok( Math::Prime::Util::_get_prime_cache_size, '>', $init_size, "Internal space grew after large precalc" ); prime_memfree; };
+is( Math::Prime::Util::_get_prime_cache_size, $init_size, "Memory freed after successful eval");
-eval { prime_precalc($bigsize); cmp_ok( Math::Prime::Util::_get_prime_cache_size, '>', $init_size, "Internal space grew after large precalc" ); die; prime_free; };
+eval { prime_precalc($bigsize); cmp_ok( Math::Prime::Util::_get_prime_cache_size, '>', $init_size, "Internal space grew after large precalc" ); die; prime_memfree; };
isnt( Math::Prime::Util::_get_prime_cache_size, $init_size, "Memory normally not freed after eval die");
-prime_free;
+prime_memfree;
-eval { my $mf = Math::Prime::Util::MemFree->new; prime_precalc($bigsize); cmp_ok( Math::Prime::Util::_get_prime_cache_size, '>', $init_size, "Internal space grew after large precalc" ); die; prime_free; };
-is( Math::Prime::Util::_get_prime_cache_size, $init_size, "Memory is not freed after eval die using object scoper");
+eval { my $mf = Math::Prime::Util::MemFree->new; prime_precalc($bigsize); cmp_ok( Math::Prime::Util::_get_prime_cache_size, '>', $init_size, "Internal space grew after large precalc" ); die; };
+is( Math::Prime::Util::_get_prime_cache_size, $init_size, "Memory is freed after eval die using object scoper");
diff --git a/t/16-randomprime.t b/t/16-randomprime.t
new file mode 100644
index 0000000..bb51aba
--- /dev/null
+++ b/t/16-randomprime.t
@@ -0,0 +1,99 @@
+#!/usr/bin/env perl
+use strict;
+use warnings;
+
+use Test::More;
+use Math::Prime::Util qw/random_prime random_ndigit_prime is_prime/;
+
+my $use64 = Math::Prime::Util::_maxbits > 32;
+my $extra = defined $ENV{RELEASE_TESTING} && $ENV{RELEASE_TESTING};
+
+plan tests => 13 + 6*1 + 11*3 + 10*2 + 12*2;
+
+my $infinity = ~0 * ~0;
+
+ok(!eval { random_prime(undef); }, "random_prime(undef)");
+ok(!eval { random_prime(-3); }, "random_prime(-3)");
+ok(!eval { random_prime("a"); }, "random_prime(a)");
+ok(!eval { random_prime(undef,undef); }, "random_prime(undef,undef)");
+ok(!eval { random_prime(2,undef); }, "random_prime(2,undef)");
+ok(!eval { random_prime(2,"a"); }, "random_prime(2,a)");
+ok(!eval { random_prime(undef,0); }, "random_prime(undef,0)");
+ok(!eval { random_prime(0,undef); }, "random_prime(0,undef)");
+ok(!eval { random_prime(2,undef); }, "random_prime(2,undef)");
+ok(!eval { random_prime(2,-4); }, "random_prime(2,-4)");
+ok(!eval { random_prime(2,$infinity); }, "random_prime(2,+infinity)");
+ok(!eval { random_prime($infinity); }, "random_prime(+infinity)");
+ok(!eval { random_prime(-$infinity); }, "random_prime(-infinity)");
+
+my %range_edge_empty = (
+ "0 to 0" => [],
+ "0 to 1" => [],
+ "2 to 1" => [],
+ "3 to 2" => [],
+ "1294268492 to 1294268778" => [],
+ "3842610774 to 3842611108" => [],
+);
+while (my($range, $expect) = each (%range_edge_empty)) {
+ my($low,$high) = $range =~ /(\d+) to (\d+)/;
+ is( random_prime($low,$high), undef, "primes($low,$high) should return undef" );
+}
+
+my %range_edge = (
+ "0 to 2" => [2,2],
+ "2 to 2" => [2,2],
+ "2 to 3" => [2,3],
+ "3 to 5" => [3,5],
+ "10 to 20" => [11,19],
+ "8 to 12" => [11,11],
+ "10 to 12" => [11,11],
+ "16706143 to 16706143" => [16706143,16706143],
+ "16706142 to 16706144" => [16706143,16706143],
+ "3842610773 to 3842611109" => [3842610773,3842611109],
+ "3842610772 to 3842611110" => [3842610773,3842611109],
+);
+while (my($range, $expect) = each (%range_edge)) {
+ my($low,$high) = $range =~ /(\d+) to (\d+)/;
+ my $got = random_prime($low,$high);
+ ok( is_prime($got), "Prime in range $low-$high is indeed prime" );
+ cmp_ok( $got, '>=', $expect->[0], "random_prime($low,$high) >= $expect->[0]");
+ cmp_ok( $got, '<=', $expect->[1], "random_prime($low,$high) >= $expect->[1]");
+}
+
+my %ranges = (
+ "2 to 20" => [2,19],
+ "3 to 7" => [3,7],
+ "20 to 100" => [23,97],
+ "5678 to 9876" => [5683,9871],
+ "27767 to 88493" => [27767,88493],
+ "27764 to 88498" => [27767,88493],
+ "27764 to 88493" => [27767,88493],
+ "27767 to 88498" => [27767,88493],
+ "17051687 to 17051899" => [17051687,17051899],
+ "17051688 to 17051898" => [17051707,17051887],
+);
+while (my($range, $expect) = each (%ranges)) {
+ my($low,$high) = $range =~ /(\d+) to (\d+)/;
+ my $isprime = 1;
+ my $inrange = 1;
+ for (1 .. 1000) {
+ my $got = random_prime($low,$high);
+ $isprime *= is_prime($got) ? 1 : 0;
+ $inrange *= (($got >= $expect->[0]) && ($got <= $expect->[1])) ? 1 : 0;
+ }
+ ok($isprime, "All returned values for $low-$high were prime" );
+ ok($inrange, "All returned values for $low-$high were in the range" );
+}
+
+my @to = (2, 3, 4, 5, 6, 7, 8, 9, 100, 1000, 1000000, 4294967295);
+foreach my $high (@to) {
+ my $isprime = 1;
+ my $inrange = 1;
+ for (1 .. 1000) {
+ my $got = random_prime($high);
+ $isprime *= is_prime($got) ? 1 : 0;
+ $inrange *= (($got >= 2) && ($got <= $high)) ? 1 : 0;
+ }
+ ok($isprime, "All returned values for $high were prime" );
+ ok($inrange, "All returned values for $high were in the range" );
+}
diff --git a/util.c b/util.c
index 9b47e74..f969af5 100644
--- a/util.c
+++ b/util.c
@@ -1,7 +1,6 @@
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
-#include <limits.h>
#include <math.h>
#include "util.h"
@@ -537,7 +536,7 @@ UV nth_prime_approx(UV n)
fn = (double) n;
flogn = log(n);
flog2n = log(flogn);
-
+
/* Cipolla 1902:
* m=0 fn * ( flogn + flog2n - 1 );
* m=1 + ((flog2n - 2)/flogn) );
--
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