Commit 274aad6d161080633221f243ec282cbb8f79565a - libmath-prime-util-perl

Imported Upstream version 0.20 gregor herrmann 11 years ago

24 changed file(s) with 792 addition(s) and 212 deletion(s). Raw diff Collapse all Expand all

-3

.travis.yml less more

13	13	- cpanm Math::MPFR
14	14
15	15	env:
	16	-
	17	- MPU_NO_GMP=1
	18	- MPU_NO_MPFR=1
	19	- MPU_NO_GMP=1 MPU_NO_MPFR=1
16	20	- MPU_NO_XS=1 MPU_NO_GMP=1 MPU_NO_MPFR=1
17		- MPU_NO_MPFR=1
18		- MPU_NO_GMP=1
19		- MPU_NO_GMP=1 MPU_NO_MPFR=1
20	21
21	22
22	23	install:

+34

-0

Changes less more

0	0	Revision history for Perl extension Math::Prime::Util.
	1
	2	0.20 3 February 2012
	3
	4	- Speedup for PP AKS, and turn off test on 32-bit machines.
	5
	6	- Replaced fast sqrt detection in PP.pm with a slightly slower version.
	7	The bloom filter doesn't work right in 32-bit Perl. These two changes
	8	should speed up testing on some machines by a huge amount.
	9
	10	- Fix is_perfect_power in XS AKS.
	11
	12	0.19 1 February 2012
	13
	14	- Update MR bases with newest from http://miller-rabin.appspot.com/.
	15
	16	- Fixed some issues when using bignum and Calc BigInt backend, and bignum
	17	and Perl 5.6.
	18
	19	- Added tests for bigint is_provable_prime.
	20
	21	- Added a few tests to give better coverage.
	22
	23	- Adjust some validation subroutines to cut down on overhead.
	24
	25	0.18 14 January 2012
	26
	27	- Add random_strong_prime.
	28
	29	- Fix builds with Solaris 9 and older.
	30
	31	- Add some debug info to perhaps find out why old ActiveState Perls are
	32	dying in Math::BigInt::Calc, as if they were using really old versions
	33	that run out of memory trying to calculate '2 ** 66'.
	34	http://code.activestate.com/ppm/Math-Prime-Util/
1	35
2	36	0.17 20 December 2012
3	37

-0

MANIFEST less more

38	38	examples/bench-pp-count.pl
39	39	examples/bench-pp-isprime.pl
40	40	examples/bench-pp-sieve.pl
	41	examples/bench-mp-nextprime.pl
	42	examples/bench-mp-psrp.pl
	43	examples/bench-mp-prime_count.pl
41	44	examples/test-factor-yafu.pl
42	45	examples/test-factor-mpxs.pl
43	46	examples/test-nextprime-yafu.pl

45	48	examples/test-holf.pl
46	49	examples/test-nthapprox.pl
47	50	examples/test-pcapprox.pl
	51	examples/test-primality-small.pl
48	52	examples/sophie_germain.pl
49	53	examples/twin_primes.pl
50	54	examples/find_mr_bases.pl

76	80	t/31-threading.t
77	81	t/50-factoring.t
78	82	t/51-primearray.t
	83	t/70-rt-bignum.t
79	84	t/80-pp.t
80	85	t/81-bignum.t
81	86	t/90-release-perlcritic.t

+11

-3

META.json less more

3	3	"Dana A Jacobsen <dana@acm.org>"
4	4	],
5	5	"dynamic_config" : 1,
6		"generated_by" : "ExtUtils::MakeMaker version 6.6302, CPAN::Meta::Converter version 2.120921",
	6	"generated_by" : "ExtUtils::MakeMaker version 6.64, CPAN::Meta::Converter version 2.120921",
7	7	"license" : [
8	8	"perl_5"
9	9	],

32	32	},
33	33	"runtime" : {
34	34	"recommends" : {
35		"Math::MPFR" : "0"
	35	"Math::BigInt::GMP" : "0",
	36	"Math::MPFR" : "2.03",
	37	"Math::Prime::Util::GMP" : "0.06"
36	38	},
37	39	"requires" : {
38	40	"Carp" : "0",

48	50	}
49	51	},
50	52	"release_status" : "stable",
51		"version" : "0.17"
	53	"resources" : {
	54	"homepage" : "https://github.com/danaj/Math-Prime-Util",
	55	"repository" : {
	56	"url" : "https://github.com/danaj/Math-Prime-Util"
	57	}
	58	},
	59	"version" : "0.20"
52	60	}

-3

META.yml less more

7	7	configure_requires:
8	8	ExtUtils::MakeMaker: 0
9	9	dynamic_config: 1
10		generated_by: 'ExtUtils::MakeMaker version 6.6302, CPAN::Meta::Converter version 2.120921'
	10	generated_by: 'ExtUtils::MakeMaker version 6.64, CPAN::Meta::Converter version 2.120921'
11	11	license: perl
12	12	meta-spec:
13	13	url: http://module-build.sourceforge.net/META-spec-v1.4.html

18	18	- t
19	19	- inc
20	20	recommends:
21		Math::MPFR: 0
	21	Math::BigInt::GMP: 0
	22	Math::MPFR: 2.03
	23	Math::Prime::Util::GMP: 0.06
22	24	requires:
23	25	Carp: 0
24	26	Config: 0

29	31	XSLoader: 0.01
30	32	base: 0
31	33	perl: 5.006002
32		version: 0.17
	34	resources:
	35	homepage: https://github.com/danaj/Math-Prime-Util
	36	repository: https://github.com/danaj/Math-Prime-Util
	37	version: 0.20

+11

-5

Makefile.PL less more

19	19
20	20	BUILD_REQUIRES=>{
21	21	'Test::More' => '0.45',
22		'bignum' => '0.22', # Used for bigint tests
	22	'bignum' => '0.22', # 'use bigint' in tests
23	23	},
24	24	PREREQ_PM => {
25	25	'Exporter' => '5.562',

32	32	'Math::BigFloat' => '1.59',
33	33	},
34	34	META_MERGE => {
35		recommends => { 'Math::Prime::Util::GMP' => 0.06, },
36		recommends => { 'Math::BigInt::GMP' => 0, },
37		recommends => { 'Math::MPFR' => 0, },
38		},
	35	resources => {
	36	homepage => 'https://github.com/danaj/Math-Prime-Util',
	37	repository => 'https://github.com/danaj/Math-Prime-Util',
	38	},
	39	recommends => {
	40	'Math::Prime::Util::GMP' => 0.06,
	41	'Math::BigInt::GMP' => 0,
	42	'Math::MPFR' => 2.03,
	43	},
	44	},
39	45
40	46	MIN_PERL_VERSION => 5.006002,
41	47	);

-1

README less more

0		Math::Prime::Util version 0.17
	0	Math::Prime::Util version 0.20
1	1
2	2	A set of utilities related to prime numbers. These include multiple sieving
3	3	methods, is_prime, prime_count, nth_prime, approximations and bounds for

-1

TODO less more

33	33
34	34	- Tests for:
35	35
36		is_provable_prime
37	36	RiemannZeta
38	37	RiemannR
39	38

44	43	tinymt32, seed it using the system rand (multiple calls, just use 8-bits
45	44	each), then use it to get 32-bit irands. This would only be used if they
46	45	didn't give us a RNG (so they don't care about strict crypto).
	46
	47	- Add PE 35 Circular primes and PE 291 Panatoipal primes to bin/primes.pl

+43

-8

aks.c less more

1	1	#include <stdlib.h>
2	2	#include <string.h>
3	3	#include <math.h>
	4	#include <float.h>
4	5
5	6	/*
6	7	* The AKS v6 algorithm, for native integers. Based on the AKS v6 paper.

36	37	return log2n;
37	38	}
38	39
39		/* See Bach and Sorenson (1993) for much better algorithm */
40		static int is_perfect_power(UV x) {
	40	/* Bach and Sorenson (1993) would be better */
	41	static int is_perfect_power(UV n) {
41	42	UV b, last;
42		if ((x & (x-1)) == 0) return 1; /* powers of 2 */
43		b = sqrt(x); if (bb == x) return 1; / perfect square */
44		last = log2floor(x) + 1;
45		for (b = 3; b < last; b = _XS_next_prime(b)) {
46		UV root = pow(x, 1.0 / (double)b);
47		if (pow(root, b) == x) return 1;
	43	if ((n <= 3) \|\| (n == UV_MAX)) return 0;
	44	if ((n & (n-1)) == 0) return 1; /* powers of 2 */
	45	last = log2floor(n-1) + 1;
	46	#if (BITS_PER_WORD == 32) \|\| (DBL_DIG > 19)
	47	if (1) {
	48	#elif DBL_DIG == 10
	49	if (n < UVCONST(10000000000)) {
	50	#elif DBL_DIG == 15
	51	if (n < UVCONST(1000000000000000)) {
	52	#else
	53	if ( n < (UV) pow(10, DBL_DIG) ) {
	54	#endif
	55	/* Simple floating point method. Fast, but need enough mantissa. */
	56	b = sqrt(n)+0.5; if (bb == n) return 1; / perfect square */
	57	for (b = 3; b < last; b = _XS_next_prime(b)) {
	58	UV root = pow(n, 1.0 / (double)b) + 0.5;
	59	if ( ((UV)(pow(root, b)+0.5)) == n) return 1;
	60	}
	61	} else {
	62	/* Dietzfelbinger, algorithm 2.3.5 (without optimized exponential) */
	63	for (b = 2; b <= last; b++) {
	64	UV a = 1;
	65	UV c = n;
	66	while (c >= HALF_WORD) c = (1+c)>>1;
	67	while ((c-a) >= 2) {
	68	UV m, maxm, p, i;
	69	m = (a+c) >> 1;
	70	maxm = UV_MAX / m;
	71	p = m;
	72	for (i = 2; i <= b; i++) {
	73	if (p > maxm) { p = n+1; break; }
	74	p *= m;
	75	}
	76	if (p == n) return 1;
	77	if (p < n)
	78	a = m;
	79	else
	80	c = m;
	81	}
	82	}
48	83	}
49	84	return 0;
50	85	}

+28

-0

examples/bench-mp-nextprime.pl less more

	0	#!/usr/bin/env perl
	1	use strict;
	2	use warnings;
	3	use Math::Prime::Util;
	4	use Math::Prime::Util::GMP;
	5	use Math::Primality;
	6	use Benchmark qw/:all/;
	7	my $count = shift \|\| -2;
	8	srand(29); # So we have repeatable results
	9
	10	test_at_digits($_, 1000) for (5, 15, 25, 50, 200);
	11
	12	sub test_at_digits {
	13	my($digits, $numbers) = @_;
	14	die "Digits must be > 0" unless $digits > 0;
	15
	16	my $start = Math::Prime::Util::random_ndigit_prime($digits) - 3;
	17	my $end = $start;
	18	$end = Math::Prime::Util::GMP::next_prime($end) for 1 .. $numbers;
	19
	20	print "next_prime x $numbers starting at $start\n";
	21
	22	cmpthese($count,{
	23	'MP' => sub { my $n = $start; $n = Math::Primality::next_prime($n) for 1..$numbers; die "MP ended with $n instead of $end" unless $n == $end; },
	24	'MPU' => sub { my $n = $start; $n = Math::Prime::Util::next_prime($n) for 1..$numbers; die "MPU ended with $n instead of $end" unless $n == $end; },
	25	'MPU GMP' => sub { my $n = $start; $n = Math::Prime::Util::GMP::next_prime($n) for 1..$numbers; die "MPU GMP ended with $n instead of $end" unless $n == $end; },
	26	});
	27	}

+21

-0

examples/bench-mp-prime_count.pl less more

	0	#!/usr/bin/env perl
	1	use strict;
	2	use warnings;
	3	use Math::Prime::Util;
	4	use Math::Prime::Util::GMP;
	5	use Math::Primality;
	6	use Benchmark qw/:all/;
	7	my $count = shift \|\| -2;
	8
	9	#my($n, $exp) = (100000,9592);
	10	my($n, $exp) = (1000000,78498);
	11	#my($n, $exp) = (10000000,664579);
	12	cmpthese($count,{
	13	'MP' =>sub { die unless $exp == Math::Primality::prime_count($n); },
	14	'MPU default' =>sub { die unless $exp == Math::Prime::Util::prime_count($n); },
	15	'MPU XS Sieve' =>sub { die unless $exp == Math::Prime::Util::_XS_prime_count($n); },
	16	'MPU XS Lehmer'=>sub { die unless $exp == Math::Prime::Util::_XS_lehmer_pi($n); },
	17	'MPU PP Sieve' =>sub { die unless $exp == Math::Prime::Util::PP::_sieve_prime_count($n); },
	18	'MPU PP Lehmer'=>sub { die unless $exp == Math::Prime::Util::PP::_lehmer_pi($n); },
	19	'MPU GMP Trial'=>sub { die unless $exp == Math::Prime::Util::GMP::prime_count(2,$n); },
	20	});

+29

-0

examples/bench-mp-psrp.pl less more

	0	#!/usr/bin/env perl
	1	use strict;
	2	use warnings;
	3	use Math::Prime::Util;
	4	use Math::Prime::Util::GMP;
	5	use Math::Primality;
	6	use Benchmark qw/:all/;
	7	use List::Util qw/min max/;
	8	my $count = shift \|\| -2;
	9	srand(29); # So we have repeatable results
	10
	11	test_at_digits($_, 1000) for (5, 15, 25, 50, 200);
	12
	13	sub test_at_digits {
	14	my($digits, $numbers) = @_;
	15	die "Digits must be > 0" unless $digits > 0;
	16
	17	# We get a mix of primes and non-primes.
	18	my @nums = map { Math::Prime::Util::random_ndigit_prime($digits)+2 } 1 .. $numbers;
	19	print "is_strong_pseudoprime for $numbers random $digits-digit numbers",
	20	" (", min(@nums), " - ", max(@nums), ")\n";
	21
	22	cmpthese($count,{
	23	'MP' =>sub {Math::Primality::is_strong_pseudoprime($_,3) for @nums;},
	24	'MPU' =>sub {Math::Prime::Util::is_strong_pseudoprime($_,3) for @nums;},
	25	'MPU PP' =>sub {Math::Prime::Util::PP::miller_rabin($_,3) for @nums;},
	26	'MPU GMP' =>sub {Math::Prime::Util::GMP::is_strong_pseudoprime($_,3) for @nums;},
	27	});
	28	}

-1

examples/find_mr_bases.pl less more

30	30
31	31	# Parallel:
32	32	my $maxn :shared;
33		my $start = 1; # People have mined below ~ 2^52
	33	my $start = int(259+241); # People have mined below 2^55
34	34	$maxn = 2047;
35	35	my $nextn = 2049;
36	36	my @threads;

82	82	2012-07-02 base 64390572806844 good up to 161701
83	83	2012-10-15 base 1769236083487960 good up to 192001
84	84	2012-10-17 base 1948244569546278 good up to 212321
	85	2013-01-14 base 34933608779780163 good up to 218245

+42

-0

examples/test-primality-small.pl less more

	0	#!/usr/bin/env perl
	1	use strict;
	2	use warnings;
	3	$\| = 1; # fast pipes
	4
	5	# Make sure the is_prob_prime functionality is working for small inputs.
	6	# Good for making sure the first few M-R bases are set up correctly.
	7
	8	use Math::Prime::Util qw/is_prob_prime/;
	9	use Math::Prime::Util::PrimeArray;
	10
	11	my @primes; tie @primes, 'Math::Prime::Util::PrimeArray';
	12
	13	# Test just primes
	14	if (0) {
	15	foreach my $i (1 .. 10000000) {
	16	my $n = shift @primes;
	17	die unless is_prob_prime($n);
	18	#print "." unless $i % 16384;
	19	print "$i $n\n" unless $i % 262144;
	20	}
	21	}
	22
	23	# Test every number up to the 100Mth prime (about 2000M)
	24	if (1) {
	25	die "2 should be prime" unless is_prob_prime(2);
	26	shift @primes;
	27	my $n = shift @primes;
	28	foreach my $i (2 .. 100_000_000) {
	29	die "$n should be prime" unless is_prob_prime($n);
	30	print "$i $n\n" unless $i % 262144;
	31	my $next = shift @primes;
	32	my $diff = ($next - $n) >> 1;
	33	if ($diff > 1) {
	34	foreach my $d (1 .. $diff-1) {
	35	my $cn = $n + 2*$d;
	36	die "$cn should be composite" if is_prob_prime($cn);
	37	}
	38	}
	39	$n = $next;
	40	}
	41	}

+18

-18

factor.c less more

207	207	}
208	208	#else
209	209	#if 1
210		/* Better bases from: http://miller-rabin.appspot.com/ */
211		if (n < UVCONST(212321)) {
212		bases[0] = UVCONST(1948244569546278);
	210	/* Better bases from http://miller-rabin.appspot.com/, 30 Jan 2013 */
	211	if (n < UVCONST(272161)) {
	212	bases[0] = UVCONST(62769592775616394);
213	213	nbases = 1;
214		} else if (n < UVCONST(360018361)) {
215		bases[0] = UVCONST( 1143370 );
216		bases[1] = UVCONST( 2350307676 );
	214	} else if (n < UVCONST(466758181)) {
	215	bases[0] = UVCONST( 91869414 );
	216	bases[1] = UVCONST( 6346128598129234 );
217	217	nbases = 2;
218		} else if (n < UVCONST(105936894253)) {
	218	} else if (n < UVCONST(109134866497)) {
219	219	bases[0] = 2;
220		bases[1] = UVCONST( 1005905886 );
221		bases[2] = UVCONST( 1340600841 );
	220	bases[1] = UVCONST( 45650740 );
	221	bases[2] = UVCONST( 3722628058 );
222	222	nbases = 3;
223		} else if (n < UVCONST(31858317218647)) {
	223	} else if (n < UVCONST(47636622961201)) {
224	224	bases[0] = 2;
225		bases[1] = UVCONST( 642735 );
226		bases[2] = UVCONST( 553174392 );
227		bases[3] = UVCONST( 3046413974 );
	225	bases[1] = UVCONST( 2570940 );
	226	bases[2] = UVCONST( 211991001 );
	227	bases[3] = UVCONST( 3749873356 );
228	228	nbases = 4;
229		} else if (n < UVCONST(3071837692357849)) {
	229	} else if (n < UVCONST(3770579582154547)) {
230	230	bases[0] = 2;
231		bases[1] = UVCONST( 75088 );
232		bases[2] = UVCONST( 642735 );
233		bases[3] = UVCONST( 203659041 );
234		bases[4] = UVCONST( 3613982119 );
	231	bases[1] = UVCONST( 2570940 );
	232	bases[2] = UVCONST( 880937 );
	233	bases[3] = UVCONST( 610386380 );
	234	bases[4] = UVCONST( 4130785767 );
235	235	nbases = 5;
236	236	} else {
237	237	bases[0] = 2;

+71

-46

lib/Math/Prime/Util/PP.pm less more

4	4
5	5	BEGIN {
6	6	$Math::Prime::Util::PP::AUTHORITY = 'cpan:DANAJ';
7		$Math::Prime::Util::PP::VERSION = '0.14';
	7	$Math::Prime::Util::PP::VERSION = '0.20';
8	8	}
9	9
10	10	# The Pure Perl versions of all the Math::Prime::Util routines.

64	64	sub _validate_positive_integer {
65	65	my($n, $min, $max) = @_;
66	66	croak "Parameter must be defined" if !defined $n;
67		croak "Parameter '$n' must be a positive integer" if $n =~ tr/0123456789//c;
	67	croak "Parameter '$n' must be a positive integer"
	68	if ref($n) ne 'Math::BigInt' && $n =~ tr/0123456789//c;
68	69	croak "Parameter '$n' must be >= $min" if defined $min && $n < $min;
69	70	croak "Parameter '$n' must be <= $max" if defined $max && $n > $max;
70	71	if ($n <= ~0) {

73	74	} elsif (ref($n) ne 'Math::BigInt') {
74	75	croak "Parameter '$n' outside of integer range" if !defined $bigint::VERSION;
75	76	$_[0] = Math::BigInt->new("$n"); # Make $n a proper bigint object
	77	$_[0]->upgrade(undef) if $_[0]->upgrade(); # Stop BigFloat upgrade
	78	} else {
	79	$_[0]->upgrade(undef) if $_[0]->upgrade(); # Stop BigFloat upgrade
76	80	}
77	81	# One of these will be true:
78	82	# 1) $n <= max and $n is not a bigint

664	668
665	669	sub _is_perfect_power {
666	670	my $n = shift;
667		my $log2n = _log2($n);
	671	return 0 if $n <= 3 \|\| $n != int($n);
	672	return 1 if ($n & ($n-1)) == 0; # Power of 2
668	673	$n = Math::BigInt->new("$n") unless ref($n) eq 'Math::BigInt';
669		for my $e (@{primes($log2n)}) {
	674	# Perl 5.6.2 chokes on this, so do it via as_bin
	675	# my $log2n = 0; { my $num = $n; $log2n++ while $num >>= 1; }
	676	my $log2n = length($n->as_bin) - 2;
	677	for (my $e = 2; $e <= $log2n; $e = next_prime($e)) {
670	678	return 1 if $n->copy()->broot($e)->bpow($e) == $n;
671	679	}
672	680	0;

843	851
844	852	# Check for perfect square
845	853	if (ref($n) eq 'Math::BigInt') {
846		my $mcheck = int(($n & 127)->bstr);
847		if (($mcheck0x8bc40d7d) & ($mcheck0xa1e2f5d1) & 0x14020a) {
848		# ~82% of non-squares were rejected by the bloom filter
	854	my $mc = int(($n & 31)->bstr);
	855	if ($mc==0\|\|$mc==1\|\|$mc==4\|\|$mc==9\|\|$mc==16\|\|$mc==17\|\|$mc==25) {
849	856	my $sq = $n->copy->bsqrt->bfloor;
850	857	$sq->bmul($sq);
851	858	return 0 if $sq == $n;
852	859	}
853	860	} else {
854		my $mcheck = $n & 127;
855		if (($mcheck0x8bc40d7d) & ($mcheck0xa1e2f5d1) & 0x14020a) {
	861	my $mc = $n & 31;
	862	if ($mc==0\|\|$mc==1\|\|$mc==4\|\|$mc==9\|\|$mc==16\|\|$mc==17\|\|$mc==25) {
856	863	my $sq = int(sqrt($n));
857	864	return 0 if ($sq*$sq) == $n;
858	865	}

870	877	# It's now time to perform the Lucas pseudoprimality test using $D.
871	878
872	879	if (ref($n) ne 'Math::BigInt') {
873		require Math::BigInt;
	880	if (!defined $MATH::BigInt::VERSION) {
	881	eval { require Math::BigInt; Math::BigInt->import(try=>'GMP,Pari'); 1; }
	882	or do { croak "Cannot load Math::BigInt "; }
	883	}
874	884	$n = Math::BigInt->new("$n");
875	885	}
876	886

973	983	for (my $ix = 0; $ix <= $px_degree; $ix++) {
974	984	my $px_at_ix = $px->[$ix];
975	985	next unless $px_at_ix;
976		foreach my $iy (@indices_y) {
977		my $py_at_iy = $py->[$iy];
978		my $rindex = ($ix + $iy) % $r; # reduce mod X^r-1
979		if (!defined $res[$rindex]) {
980		$res[$rindex] = $_poly_bignum ? Math::BigInt->bzero : 0
981		}
982		$res[$rindex] = ($res[$rindex] + ($py_at_iy * $px_at_ix)) % $n;
	986	if ($_poly_bignum) {
	987	foreach my $iy (@indices_y) {
	988	my $py_px = $py->[$iy] * $px_at_ix;
	989	my $rindex = ($ix + $iy) % $r; # reduce mod X^r-1
	990	$res[$rindex] = Math::BigInt->bzero unless defined $res[$rindex];
	991	$res[$rindex]->badd($py_px)->bmod($n);
	992	}
	993	} else {
	994	foreach my $iy (@indices_y) {
	995	my $py_px = $py->[$iy] * $px_at_ix;
	996	my $rindex = ($ix + $iy) % $r; # reduce mod X^r-1
	997	$res[$rindex] = 0 unless defined $res[$rindex];
	998	$res[$rindex] = ($res[$rindex] + $py_px) % $n;
	999	}
983	1000	}
984	1001	}
985	1002	# In case we had upper terms go to zero after modulo, reduce the degree.

1012	1029	sub is_aks_prime {
1013	1030	my $n = shift;
1014	1031
1015		if (ref($n) ne 'Math::BigInt') {
1016		if (!defined $MATH::BigInt::VERSION) {
1017		eval { require Math::BigInt; Math::BigInt->import(try=>'GMP,Pari'); 1; }
1018		or do { croak "Cannot load Math::BigInt "; }
1019		}
1020		if (!defined $MATH::BigFloat::VERSION) {
1021		eval { require Math::BigFloat; Math::BigFloat->import(); 1; }
1022		or do { croak "Cannot load Math::BigFloat "; }
1023		}
1024		$n = Math::BigInt->new("$n");
1025		}
	1032	if (!defined $MATH::BigInt::VERSION) {
	1033	eval { require Math::BigInt; Math::BigInt->import(try=>'GMP,Pari'); 1; }
	1034	or do { croak "Cannot load Math::BigInt "; }
	1035	}
	1036	if (!defined $MATH::BigFloat::VERSION) {
	1037	eval { require Math::BigFloat; Math::BigFloat->import(); 1; }
	1038	or do { croak "Cannot load Math::BigFloat "; }
	1039	}
	1040	$n = Math::BigInt->new("$n") unless ref($n) eq 'Math::BigInt';
1026	1041
1027	1042	return 0 if $n < 2;
1028	1043	return 0 if _is_perfect_power($n);

1071	1086	return ($_[0]) if $_[0] < 4;
1072	1087
1073	1088	my @factors;
1074		while ( !($_[0] % 2) ) { push @factors, 2; $_[0] = int($_[0] / 2); }
1075		while ( !($_[0] % 3) ) { push @factors, 3; $_[0] = int($_[0] / 3); }
1076		while ( !($_[0] % 5) ) { push @factors, 5; $_[0] = int($_[0] / 5); }
1077
1078		# Stop using bignum if we can
1079		$_[0] = int($_[0]->bstr) if ref($_[0]) eq 'Math::BigInt' && $_[0] <= ~0;
	1089	if (ref($_[0]) ne 'Math::BigInt') {
	1090	while ( !($_[0] % 2) ) { push @factors, 2; $_[0] = int($_[0] / 2); }
	1091	while ( !($_[0] % 3) ) { push @factors, 3; $_[0] = int($_[0] / 3); }
	1092	while ( !($_[0] % 5) ) { push @factors, 5; $_[0] = int($_[0] / 5); }
	1093	} else {
	1094	if (Math::BigInt::bgcd($_[0], 235) != 1) {
	1095	while ( $_[0]->is_even) { push @factors, 2; $_[0]->brsft(1); }
	1096	foreach my $div (3, 5) {
	1097	my ($q, $r) = $_[0]->copy->bdiv($div);
	1098	while ($r->is_zero) {
	1099	push @factors, $div;
	1100	$_[0] = $q;
	1101	($q, $r) = $_[0]->copy->bdiv($div);
	1102	}
	1103	}
	1104	}
	1105	$_[0] = int($_[0]->bstr) if $_[0] <= ~0;
	1106	}
1080	1107
1081	1108	if ( ($_[0] > 1) && _is_prime7($_[0]) ) {
1082	1109	push @factors, $_[0];

1209	1236	if ($f == $n) {
1210	1237	last if $inloop++; # We've been here before
1211	1238	} elsif ($f != 1) {
1212		my $f2 = $n->copy->bdiv($f);
	1239	my $f2 = $n->copy->bdiv($f)->as_int;
1213	1240	push @factors, $f;
1214	1241	push @factors, $f2;
1215	1242	croak "internal error in prho" unless ($f * $f2) == $n;

1282	1309	$Xi->bmul($Xi); $Xi->badd($a); $Xi->bmod($n);
1283	1310	my $f = Math::BigInt::bgcd( ($Xi > $Xm) ? $Xi-$Xm : $Xm-$Xi, $n);
1284	1311	if ( ($f != 1) && ($f != $n) ) {
1285		my $f2 = $n->copy->bdiv($f);
	1312	my $f2 = $n->copy->bdiv($f)->as_int;
1286	1313	push @factors, $f;
1287	1314	push @factors, $f2;
1288	1315	croak "internal error in pbrent" unless ($f * $f2) == $n;

1342	1369	$kf = $n if $kf == 0;
1343	1370	my $f = Math::BigInt::bgcd( $kf-1, $n );
1344	1371	if ( ($f != 1) && ($f != $n) ) {
1345		my $f2 = $n->copy->bdiv($f);
	1372	my $f2 = $n->copy->bdiv($f)->as_int;
1346	1373	push @factors, $f;
1347	1374	push @factors, $f2;
1348	1375	croak "internal error in pminus1" unless ($f * $f2) == $n;

1379	1406	if ( ref($n) eq 'Math::BigInt' ) {
1380	1407	for my $i ($startrounds .. $rounds) {
1381	1408	my $ni = $n->copy->bmul($i);
1382		my $s = $ni->copy->bsqrt->bfloor;
	1409	my $s = $ni->copy->bsqrt->bfloor->as_int;
1383	1410	$s->binc if ($s * $s) != $ni;
1384	1411	my $m = $s->copy->bmul($s)->bmod($n);
1385	1412	# Check for perfect square
1386		my $mcheck = int(($m & 127)->bstr);
1387		next if (($mcheck0x8bc40d7d) & ($mcheck0xa1e2f5d1) & 0x14020a);
1388		# ... 82% of non-squares were rejected by the bloom filter
1389		my $f = $m->copy->bsqrt->bfloor;
	1413	my $mc = int(($m & 31)->bstr);
	1414	next unless $mc==0\|\|$mc==1\|\|$mc==4\|\|$mc==9\|\|$mc==16\|\|$mc==17\|\|$mc==25;
	1415	my $f = $m->copy->bsqrt->bfloor->as_int;
1390	1416	next unless ($f*$f) == $m;
1391	1417	$f = Math::BigInt::bgcd( ($s > $f) ? $s-$f : $f-$s, $n);
1392	1418	last if $f == 1 \|\| $f == $n; # Should never happen
1393		my $f2 = $n->copy->bdiv($f);
	1419	my $f2 = $n->copy->bdiv($f)->as_int;
1394	1420	push @factors, $f;
1395	1421	push @factors, $f2;
1396	1422	croak "internal error in HOLF" unless ($f * $f2) == $n;

1403	1429	$s++ if ($s * $s) != ($n * $i);
1404	1430	my $m = ($s < $_half_word) ? ($s*$s) % $n : _mulmod($s, $s, $n);
1405	1431	# Check for perfect square
1406		my $mcheck = $m & 127;
1407		next if (($mcheck0x8bc40d7d) & ($mcheck0xa1e2f5d1) & 0x14020a);
1408		# ... 82% of non-squares were rejected by the bloom filter
	1432	my $mc = $m & 31;
	1433	next unless $mc==0\|\|$mc==1\|\|$mc==4\|\|$mc==9\|\|$mc==16\|\|$mc==17\|\|$mc==25;
1409	1434	my $f = int(sqrt($m));
1410	1435	next unless $f*$f == $m;
1411	1436	$f = _gcd_ui( ($s > $f) ? $s - $f : $f - $s, $n);

+297

-109

lib/Math/Prime/Util.pm less more

4	4
5	5	BEGIN {
6	6	$Math::Prime::Util::AUTHORITY = 'cpan:DANAJ';
7		$Math::Prime::Util::VERSION = '0.17';
	7	$Math::Prime::Util::VERSION = '0.20';
8	8	}
9	9
10	10	# parent is cleaner, and in the Perl 5.10.1 / 5.12.0 core, but not earlier.

21	21	next_prime prev_prime
22	22	prime_count prime_count_lower prime_count_upper prime_count_approx
23	23	nth_prime nth_prime_lower nth_prime_upper nth_prime_approx
24		random_prime random_ndigit_prime random_nbit_prime random_maurer_prime
	24	random_prime random_ndigit_prime random_nbit_prime
	25	random_strong_prime random_maurer_prime
25	26	primorial pn_primorial
26	27	factor all_factors
27	28	moebius euler_phi jordan_totient

82	83	*fermat_factor = \&Math::Prime::Util::PP::fermat_factor;
83	84	*holf_factor = \&Math::Prime::Util::PP::holf_factor;
84	85	*squfof_factor = \&Math::Prime::Util::PP::squfof_factor;
	86	*rsqufof_factor = \&Math::Prime::Util::PP::squfof_factor;
85	87	*pbrent_factor = \&Math::Prime::Util::PP::pbrent_factor;
86	88	*prho_factor = \&Math::Prime::Util::PP::prho_factor;
87	89	*pminus1_factor = \&Math::Prime::Util::PP::pminus1_factor;

206	208	1;
207	209	}
208	210
	211	my $_bigint_small;
209	212	sub _validate_positive_integer {
210	213	my($n, $min, $max) = @_;
211	214	croak "Parameter must be defined" if !defined $n;
212		croak "Parameter '$n' must be a positive integer" if $n =~ tr/0123456789//c;
	215	if (ref($n) eq 'Math::BigInt') {
	216	croak "Parameter '$n' must be a positive integer" unless $n->sign() eq '+';
	217	} else {
	218	croak "Parameter '$n' must be a positive integer"
	219	if $n eq '' \|\| $n =~ tr/0123456789//c;
	220	}
213	221	croak "Parameter '$n' must be >= $min" if defined $min && $n < $min;
214	222	croak "Parameter '$n' must be <= $max" if defined $max && $n > $max;
215		# The second term is used instead of '<=' to fix strings like ~0+delta.
216		# The third works around a rare BigInt bug (e.g. 23 > 18446744073709551615 !!)
217		if ($n < $_Config{'maxparam'} \|\| int($n) eq $_Config{'maxparam'} \|\| "$n" < $_Config{'maxparam'}) {
218		$_[0] = $_[0]->as_number() if ref($_[0]) eq 'Math::BigFloat';
219		$_[0] = int($_[0]->bstr) if ref($_[0]) eq 'Math::BigInt';
220		} elsif (ref($n) ne 'Math::BigInt') {
221		croak "Parameter '$n' outside of integer range" if !defined $bigint::VERSION;
222		$_[0] = Math::BigInt->new("$n"); # Make $n a proper bigint object
	223
	224	if (ref($_[0])) {
	225	$_[0] = Math::BigInt->new("$_[0]") unless ref($_[0]) eq 'Math::BigInt';
	226	$_bigint_small = Math::BigInt->new("$_Config{'maxparam'}")
	227	unless defined $_bigint_small;
	228	if ($_[0]->bacmp($_bigint_small) <= 0) {
	229	$_[0] = int($_[0]->bstr);
	230	} else {
	231	$_[0]->upgrade(undef) if $_[0]->upgrade(); # Stop BigFloat upgrade
	232	}
	233	} else {
	234	# The second term is used instead of '<=' to fix strings like ~0+delta.
	235	if ( ! ($n < $_Config{'maxparam'} \|\| int($n) eq $_Config{'maxparam'}) ) {
	236	# We were handed a string representing a big number.
	237	croak "Parameter '$n' outside of integer range" if !defined $bigint::VERSION;
	238	$_[0] = Math::BigInt->new("$n"); # Make $n a proper bigint object
	239	$_[0]->upgrade(undef) if $_[0]->upgrade(); # Stop BigFloat upgrade
	240	}
223	241	}
224	242	# One of these will be true:
225	243	# 1) $n <= max and $n is not a bigint

445	463	my $p1 = primorial(Math::BigInt->new(2052));
446	464	my $p2 = primorial(Math::BigInt->new(6028));
447	465	my $p3 = primorial(Math::BigInt->new($_big_gcd_top));
448		$_big_gcd[0] = $p0 / 223092870;
449		$_big_gcd[1] = $p1 / $p0;
450		$_big_gcd[2] = $p2 / $p1;
451		$_big_gcd[3] = $p3 / $p2;
	466	$_big_gcd[0] = $p0->bdiv(223092870)->bfloor->as_int;
	467	$_big_gcd[1] = $p1->bdiv($p0)->bfloor->as_int;
	468	$_big_gcd[2] = $p2->bdiv($p1)->bfloor->as_int;
	469	$_big_gcd[3] = $p3->bdiv($p2)->bfloor->as_int;
452	470	}
453	471
454	472	# Returns a function that will get a uniform random number between 0 and

517	535	$U = ($irandf->() << 32) + $irandf->();
518	536	} while $U >= (18446744073709551615 - $remainder);
519	537	}
520		#return $U % $range;
521		$U %= $range;
522		die "random failure: $max $U\n" if $U > $max;
523		return $U;
	538	return $U % $range;
524	539	};
525	540	return $randf;
526	541	}

550	565	}
551	566
552	567	$low-- if $low == 2; # Low of 2 becomes 1 for our program.
553		croak "Invalid _random_prime parameters" if ($low % 2) == 0 \|\| ($high % 2) == 0;
	568	confess "Invalid _random_prime parameters: $low, $high" if ($low % 2) == 0 \|\| ($high % 2) == 0;
554	569
555	570	# We're going to look at the odd numbers only.
556	571	#my $range = $high - $low + 1;

637	652	my($nbins, $rem);
638	653	($nbins, $rem) = $oddrange->copy->bdiv( "$rand_part_size" );
639	654	$nbins++ if $rem > 0;
	655	$nbins = $nbins->as_int();
640	656	($binsize,$rem) = $oddrange->copy->bdiv($nbins);
641	657	$binsize++ if $rem > 0;
642		$nparts = $oddrange->copy->bdiv($binsize);
	658	$binsize = $binsize->as_int();
	659	$nparts = $oddrange->copy->bdiv($binsize)->as_int();
643	660	$low = $high->copy->bzero->badd($low) if ref($low) ne 'Math::BigInt';
644	661	} else {
645	662	my $nbins = int($oddrange / $rand_part_size);

790	807	_validate_positive_integer($bits, 2);
791	808
792	809	if (!defined $_random_nbit_ranges[$bits]) {
793		my $bigbits = $bits > $_Config{'maxbits'};
	810	my $bigbits = $bits > $_Config{'maxbits'}; # \|\| ($] < 5.8 && $bits > 49);
794	811	if ($bigbits) {
795	812	if (!defined $Math::BigInt::VERSION) {
796	813	eval { require Math::BigInt; Math::BigInt->import(try=>'GMP,Pari'); 1; }

823	840	# could go even higher if we used is_provable_prime or looked for is_prime
824	841	# returning 2. This should be reasonably fast to ~128 bits with MPU::GMP.
825	842	my $p0 = $_Config{'maxbits'};
	843	$p0 = 32 if $] < 5.8 && $p0 > 32;
826	844
827	845	return random_nbit_prime($k) if $k <= $p0;
828	846

857	875	# I've seen +0, +1, and +2 here. Maurer uses +0. Menezes uses +1.
858	876	my $q = random_maurer_prime( ($r * $k)->bfloor + 1 );
859	877	$q = Math::BigInt->new("$q") unless ref($q) eq 'Math::BigInt';
860		my $I = Math::BigInt->new(2)->bpow($k-2)->bdiv($q)->bfloor;
	878	my $I = Math::BigInt->new(2)->bpow($k-2)->bdiv($q)->bfloor->as_int();
861	879	print "B = $B r = $r k = $k q = $q I = $I\n" if $verbose && $verbose != 3;
862	880
863	881	# Big GCD's are hugely fast with GMP or Pari, but super slow with Calc.

939	957	}
940	958	croak "Failure in random_maurer_prime, could not find a prime\n";
941	959	} # End of random_maurer_prime
	960
	961	# Gordon's algorithm for generating a strong prime.
	962	sub random_strong_prime {
	963	my($t) = @_;
	964	_validate_positive_integer($t, 128);
	965
	966	if (!defined $Math::BigInt::VERSION) {
	967	eval { require Math::BigInt; Math::BigInt->import(try=>'GMP,Pari'); 1; }
	968	or do { croak "Cannot load Math::BigInt"; };
	969	}
	970	my $irandf = _get_rand_func();
	971
	972	my $l = (($t+1) >> 1) - 2;
	973	my $lp = int($t/2) - 20;
	974	my $lpp = $l - 20;
	975	while (1) {
	976	my $qp = random_nbit_prime($lp);
	977	my $qpp = random_nbit_prime($lpp);
	978	$qp = Math::BigInt->new("$qp") unless ref($qp) eq 'Math::BigInt';
	979	$qpp = Math::BigInt->new("$qpp") unless ref($qpp) eq 'Math::BigInt';
	980	my ($il, $rem) = Math::BigInt->new(2)->bpow($l-1)->bsub(1)->bdiv(2*$qpp);
	981	$il++ if $rem > 0;
	982	$il = $il->as_int();
	983	my $iu = Math::BigInt->new(2)->bpow($l)->bsub(2)->bdiv(2*$qpp)->as_int();
	984	my $istart = $il + $irandf->($iu - $il);
	985	for (my $i = $istart; $i <= $iu; $i++) { # Search for q
	986	my $q = 2 * $i * $qpp + 1;
	987	next unless is_prob_prime($q);
	988	my $pp = $qp->copy->bmodpow($q-2, $q)->bmul(2)->bmul($qp)->bsub(1);
	989	my ($jl, $rem) = Math::BigInt->new(2)->bpow($t-1)->bsub($pp)->bdiv(2$q$qp);
	990	$jl++ if $rem > 0;
	991	$jl = $jl->as_int();
	992	my $ju = Math::BigInt->new(2)->bpow($t)->bsub(1)->bsub($pp)->bdiv(2$q$qp)->as_int();
	993	my $jstart = $jl + $irandf->($ju - $jl);
	994	for (my $j = $jstart; $j <= $ju; $j++) { # Search for p
	995	my $p = $pp + 2 * $j * $q * $qp;
	996	return $p if is_prob_prime($p);
	997	}
	998	}
	999	}
	1000	}
942	1001
943	1002	} # end of the random prime section
944	1003

1156	1215	return 0 if defined $n && $n < 2;
1157	1216	_validate_positive_integer($n);
1158	1217
1159		return _XS_is_prime($n) if $n <= $_XS_MAXVAL;
	1218	return _XS_is_prime($n) if ref($n) ne 'Math::BigInt' && $n <= $_XS_MAXVAL;
1160	1219	return Math::Prime::Util::GMP::is_prime($n) if $_HAVE_GMP;
1161	1220	return is_prob_prime($n);
1162	1221	}

1178	1237	_validate_positive_integer($n);
1179	1238
1180	1239	# If we have XS and n is either small or bigint is unknown, then use XS.
1181		return _XS_next_prime($n) if $n <= $_XS_MAXVAL
	1240	return _XS_next_prime($n) if ref($n) ne 'Math::BigInt' && $n <= $_XS_MAXVAL
1182	1241	&& (!defined $bigint::VERSION \|\| $n < $_Config{'maxprime'} );
1183	1242
1184	1243	# Try to stick to the plan with respect to maximum return values.

1197	1256	my($n) = @_;
1198	1257	_validate_positive_integer($n);
1199	1258
1200		return _XS_prev_prime($n) if $n <= $_XS_MAXVAL;
	1259	return _XS_prev_prime($n) if ref($n) ne 'Math::BigInt' && $n <= $_XS_MAXVAL;
1201	1260	if ($_HAVE_GMP) {
1202	1261	# If $n is a bigint object, try to make the return value the same
1203	1262	return (ref($n) eq 'Math::BigInt')

1254	1313	my($n) = @_;
1255	1314	_validate_positive_integer($n);
1256	1315
1257		return _XS_factor($n) if $n <= $_XS_MAXVAL;
	1316	return _XS_factor($n) if ref($n) ne 'Math::BigInt' && $n <= $_XS_MAXVAL;
1258	1317
1259	1318	if ($_HAVE_GMP) {
1260	1319	my @factors = Math::Prime::Util::GMP::factor($n);

1271	1330	my($n) = shift;
1272	1331	_validate_positive_integer($n);
1273	1332	# validate bases?
1274		return _XS_miller_rabin($n, @_) if $n <= $_XS_MAXVAL;
	1333	return _XS_miller_rabin($n, @_) if ref($n) ne 'Math::BigInt' && $n <= $_XS_MAXVAL;
1275	1334	return Math::Prime::Util::GMP::is_strong_pseudoprime($n, @_) if $_HAVE_GMP;
1276	1335	return Math::Prime::Util::PP::miller_rabin($n, @_);
1277	1336	}

1317	1376	return 0 if defined $n && $n < 2;
1318	1377	_validate_positive_integer($n);
1319	1378
1320		return _XS_is_prob_prime($n) if $n <= $_XS_MAXVAL;
	1379	return _XS_is_prob_prime($n) if ref($n) ne 'Math::BigInt' && $n <= $_XS_MAXVAL;
1321	1380	return Math::Prime::Util::GMP::is_prob_prime($n) if $_HAVE_GMP;
1322	1381
1323	1382	return 2 if $n == 2 \|\| $n == 3 \|\| $n == 5 \|\| $n == 7;

1359	1418	_validate_positive_integer($n);
1360	1419
1361	1420	# Shortcut some of the calls.
1362		return _XS_is_prime($n) if $n <= $_XS_MAXVAL;
	1421	return _XS_is_prime($n) if ref($n) ne 'Math::BigInt' && $n <= $_XS_MAXVAL;
1363	1422	return Math::Prime::Util::GMP::is_provable_prime($n) if $_HAVE_GMP
1364	1423	&& defined &Math::Prime::Util::GMP::is_provable_prime;
1365	1424

1392	1451	next if $ap->copy->bmodpow($nm1, $n) != 1;
1393	1452	# 2. a^((n-1)/f) != 1 mod n for all f.
1394	1453	next if (scalar grep { $_ == 1 }
1395		map { $ap->copy->bmodpow($nm1/$_,$n); }
	1454	map { $ap->copy->bmodpow(int($nm1/$_),$n); }
1396	1455	@factors) > 0;
1397	1456	return 2;
1398	1457	}

1746	1805
1747	1806	=head1 VERSION
1748	1807
1749		Version 0.17
	1808	Version 0.20
1750	1809
1751	1810
1752	1811	=head1 SYNOPSIS

1849	1908	my $rand_prime = random_prime(100, 10000); # random prime within a range
1850	1909	my $rand_prime = random_ndigit_prime(6); # random 6-digit prime
1851	1910	my $rand_prime = random_nbit_prime(128); # random 128-bit prime
	1911	my $rand_prime = random_strong_prime(256); # random 256-bit strong prime
1852	1912	my $rand_prime = random_maurer_prime(256); # random 256-bit provable prime
1853	1913
1854	1914

1918	1978
1919	1979	=over 4
1920	1980
1921		=item Install L<Math::Prime::Util::GMP>, as that will vastly increase the
1922		speed of many of the functions. This does require the
1923		L<GMP\|gttp://gmplib.org> library be installed on your system, but this
1924		increasingly comes pre-installed or easily available using the OS vendor
1925		package installation tool.
1926
1927		=item Install and use L<Math::BigInt::GMP> or L<Math::BigInt::Pari>, then use
1928		C<use bigint try =E<gt> 'GMP,Pari'> in your script, or on the command
1929		line C<-Mbigint=lib,GMP>. Large modular exponentiation is much faster
1930		using the GMP or Pari backends, as are the math and approximation
1931		functions when called with very large inputs.
1932
1933		=item Install L<Math::MPFR> if you use the Ei, li, Zeta, or R functions. If
1934		that module can be loaded, these functions will run much faster on
1935		bignum inputs, and are able to provide higher accuracy.
1936
1937		=item Having run these functions on many versions of Perl, if you're using
1938		anything older than Perl 5.14, I would recommend you upgrade if you
1939		are using bignums a lot. There are some brittle behaviors on
1940		5.12.4 and earlier with bignums.
	1981	=item *
	1982
	1983	Install L<Math::Prime::Util::GMP>, as that will vastly increase the speed
	1984	of many of the functions. This does require the L<GMP\|gttp://gmplib.org>
	1985	library be installed on your system, but this increasingly comes
	1986	pre-installed or easily available using the OS vendor package installation tool.
	1987
	1988	=item *
	1989
	1990	Install and use L<Math::BigInt::GMP> or L<Math::BigInt::Pari>, then use
	1991	C<use bigint try =E<gt> 'GMP,Pari'> in your script, or on the command line
	1992	C<-Mbigint=lib,GMP>. Large modular exponentiation is much faster using the
	1993	GMP or Pari backends, as are the math and approximation functions when
	1994	called with very large inputs.
	1995
	1996	=item *
	1997
	1998	Install L<Math::MPFR> if you use the Ei, li, Zeta, or R functions. If that
	1999	module can be loaded, these functions will run much faster on bignum inputs,
	2000	and are able to provide higher accuracy.
	2001
	2002	=item *
	2003
	2004	Having run these functions on many versions of Perl, if you're using anything
	2005	older than Perl 5.14, I would recommend you upgrade if you are using bignums
	2006	a lot. There are some brittle behaviors on 5.12.4 and earlier with bignums.
1941	2007
1942	2008	=back
1943	2009

2444	2510	L<Math::BigInt::GMP>.
2445	2511
2446	2512
	2513	=head2 random_strong_prime
	2514
	2515	my $bigprime = random_strong_prime(512);
	2516
	2517	Constructs an n-bit strong prime using Gordon's algorithm. We consider a
	2518	strong prime I<p> to be one where
	2519
	2520	=over
	2521
	2522	=item * I<p> is large. This function requires at least 128 bits.
	2523
	2524	=item * I<p-1> has a large prime factor I<r>.
	2525
	2526	=item * I<p+1> has a large prime factor I<s>
	2527
	2528	=item * I<r-1> has a large prime factor I<t>
	2529
	2530	=back
	2531
	2532	Using a strong prime in cryptography guards against easy factoring with
	2533	algorithms like Pollard's Rho. Rivest and Silverman (1999) present a case
	2534	that using strong primes is unnecessary, and most modern cryptographic systems
	2535	agree. First, the smoothness does not affect more modern factoring methods
	2536	such as ECM. Second, modern factoring methods like GNFS are far faster than
	2537	either method so make the point moot. Third, due to key size growth and
	2538	advances in factoring and attacks, for practical purposes, using large random
	2539	primes offer security equivalent to using strong primes.
	2540
	2541
2447	2542	=head2 random_maurer_prime
2448	2543
2449	2544	my $bigprime = random_maurer_prime(512);

2825	2920	# perl -MMath::Pari=:int,PARI,nextprime -E 'my $start = PARI "100000000000000000000"; my $end = $start+1000; my $p=nextprime($start); while ($p <= $end) { say $p; $p = nextprime($p+1); }'
2826	2921
2827	2922
	2923	Project Euler, problem 3 (Largest prime factor):
	2924
	2925	use Math::Prime::Util qw/factor/;
	2926	use bigint; # Only necessary for 32-bit machines.
	2927	say "", (factor(600851475143))[-1]
	2928
	2929	Project Euler, problem 7 (10001st prime):
	2930
	2931	use Math::Prime::Util qw/nth_prime/;
	2932	say nth_prime(10_001);
	2933
	2934	Project Euler, problem 10 (summation of primes):
	2935
	2936	use Math::Prime::Util qw/primes/;
	2937	my $sum = 0;
	2938	$sum += $_ for @{primes(2_000_000)};
	2939	say $sum;
	2940
	2941	Project Euler, problem 21 (Amicable numbers):
	2942
	2943	use Math::Prime::Util qw/divisor_sum/;
	2944	sub dsum { my $n = shift; divisor_sum($n, sub {$_[0]}) - $n; }
	2945	my $sum = 0;
	2946	foreach my $a (1..10000) {
	2947	my $b = dsum($a);
	2948	$sum += $a + $b if $b > $a && dsum($b) == $a;
	2949	}
	2950	say $sum;
	2951
	2952	Project Euler, problem 41 (Pandigital prime), brute force command line:
	2953
	2954	perl -MMath::Prime::Util=:all -E 'my @p = grep { /1/&&/2/&&/3/&&/4/&&/5/&&/6/&&/7/} @{primes(1000000,9999999)}; say $p[-1];
	2955
	2956	Project Euler, problem 47 (Distinct primes factors):
	2957
	2958	use Math::Prime::Util qw/factor/;
	2959	use List::MoreUtils qw/distinct/;
	2960	sub nfactors { scalar distinct factor(shift); }
	2961	my $n = pn_primorial(4); # Start with the first 4-factor number
	2962	$n++ while (nfactors($n) != 4 \|\| nfactors($n+1) != 4 \|\| nfactors($n+2) != 4 \|\| nfactors($n+3) != 4);
	2963	say $n;
	2964
	2965	Project Euler, problem 69, stupid brute force solution (about 5 seconds):
	2966
	2967	use Math::Prime::Util qw/euler_phi/;
	2968	my ($n, $max) = (0,0);
	2969	do {
	2970	my $ndivphi = $_/euler_phi($_);
	2971	($n, $max) = ($_, $ndivphi) if $ndivphi > $max;
	2972	} for 1..1000000;
	2973	say "$n $max";
	2974
	2975	Here's the right way to do PE problem 69 (under 0.03s):
	2976
	2977	use Math::Prime::Util qw/pn_primorial/;
	2978	my $n = 0;
	2979	$n++ while pn_primorial($n+1) < 1000000;
	2980	say pn_primorial($n);'
	2981
	2982
2828	2983	=head1 LIMITATIONS
2829	2984
2830	2985	I have not completed testing all the functions near the word size limit

2879	3034	--------- -------------------------- ------- -----------
2880	3035	0.03 Math::Prime::Util 0.12 using Lehmer's method
2881	3036	0.28 Math::Prime::Util 0.17 segmented mod-30 sieve
2882		0.52 Math::Prime::Util::PP 0.14 Perl (Lehmer's method)
	3037	0.47 Math::Prime::Util::PP 0.14 Perl (Lehmer's method)
2883	3038	0.9 Math::Prime::Util 0.01 mod-30 sieve
2884	3039	2.9 Math::Prime::FastSieve 0.12 decent odd-number sieve
2885	3040	11.7 Math::Prime::XS 0.29 "" but needs a count API
2886	3041	15.0 Bit::Vector 7.2
2887		57.3 Math::Prime::Util::PP 0.14 Perl (fastest I know of)
	3042	48.9 Math::Prime::Util::PP 0.14 Perl (fastest I know of)
2888	3043	170.0 Faster Perl sieve (net) 2012-01 array of odds
2889	3044	548.1 RosettaCode sieve (net) 2012-06 simplistic Perl
	3045	3048.1 Math::Primality 0.08 Perl + Math::GMPz
2890	3046	~11000 Math::Primality 0.04 Perl + Math::GMPz
2891	3047	>20000 Math::Big 1.12 Perl, > 26GB RAM used
2892	3048

2919	3075
2920	3076	=over 4
2921	3077
2922		=item L<Math::Prime::Util> looks in the sieve for a fast bit lookup if that
2923		exists (default up to 30,000 but it can be expanded, e.g.
2924		C<prime_precalc>), uses trial division for numbers higher than this but
2925		not too large (0.1M on 64-bit machines, 100M on 32-bit machines), a
2926		deterministic set of Miller-Rabin tests for 64-bit and smaller numbers,
2927		and a BPSW test for bigints.
2928
2929		=item L<Math::Prime::XS> does trial divisions, which is wonderful if the input
2930		has a small factor (or is small itself). But if given a large prime it
2931		can take orders of magnitude longer. It does not support bigints.
2932
2933		=item L<Math::Prime::FastSieve> only works in a sieved range, which is really
2934		fast if you can do it (M::P::U will do the same if you call
2935		C<prime_precalc>). Larger inputs just need too much time and memory
2936		for the sieve.
2937
2938		=item L<Math::Primality> uses GMP for all work. Under ~32-bits it uses 2 or 3
2939		MR tests, while above 4759123141 it performs a BPSW test. This is is
2940		fantastic for bigints over 2^64, but it is significantly slower than
2941		native precision tests. With 64-bit numbers it is generally an order of
2942		magnitude or more slower than any of the others. Once bigints are being
2943		used, its performance is quite good. It is faster than this module unless
2944		L<Math::Prime::Util::GMP> has been installed, in which case this module
2945		is just a little bit faster.
2946
2947		=item L<Math::Pari> has some very effective code, but it has some overhead to
2948		get to it from Perl. That means for small numbers it is relatively slow:
2949		an order of magnitude slower than M::P::XS and M::P::Util (though arguably
2950		this is only important for benchmarking since "slow" is ~2 microseconds).
2951		Large numbers transition over to smarter tests so don't slow down much.
2952		The C<ispseudoprime(n,0)> function will perform the BPSW test and is
2953		fast even for large inputs.
	3078	=item L<Math::Prime::Util>
	3079
	3080	looks in the sieve for a fast bit lookup if that exists (default up to 30,000
	3081	but it can be expanded, e.g. C<prime_precalc>), uses trial division for
	3082	numbers higher than this but not too large (0.1M on 64-bit machines,
	3083	100M on 32-bit machines), a deterministic set of Miller-Rabin tests for
	3084	64-bit and smaller numbers, and a BPSW test for bigints.
	3085
	3086	=item L<Math::Prime::XS>
	3087
	3088	does trial divisions, which is wonderful if the input has a small factor
	3089	(or is small itself). But if given a large prime it can take orders of
	3090	magnitude longer. It does not support bigints.
	3091
	3092	=item L<Math::Prime::FastSieve>
	3093
	3094	only works in a sieved range, which is really fast if you can do it
	3095	(M::P::U will do the same if you call C<prime_precalc>). Larger inputs
	3096	just need too much time and memory for the sieve.
	3097
	3098	=item L<Math::Primality>
	3099
	3100	uses GMP for all work. Under ~32-bits it uses 2 or 3 MR tests, while
	3101	above 4759123141 it performs a BPSW test. This is is fantastic for
	3102	bigints over 2^64, but it is significantly slower than native precision
	3103	tests. With 64-bit numbers it is generally an order of magnitude or more
	3104	slower than any of the others. Once bigints are being used, its
	3105	performance is quite good. It is faster than this module unless
	3106	L<Math::Prime::Util::GMP> has been installed, in which case this module
	3107	is just a little bit faster.
	3108
	3109	=item L<Math::Pari>
	3110
	3111	has some very effective code, but it has some overhead to get to it from
	3112	Perl. That means for small numbers it is relatively slow: an order of
	3113	magnitude slower than M::P::XS and M::P::Util (though arguably this is
	3114	only important for benchmarking since "slow" is ~2 microseconds). Large
	3115	numbers transition over to smarter tests so don't slow down much. The
	3116	C<ispseudoprime(n,0)> function will perform the BPSW test and is fast
	3117	even for large inputs.
2954	3118
2955	3119	=back
2956	3120

3018	3182
3019	3183	=over 4
3020	3184
3021		=item Pierre Dusart, "Estimates of Some Functions Over Primes without R.H.", preprint, 2010. L<http://arxiv.org/abs/1002.0442/>
3022
3023		=item Pierre Dusart, "Autour de la fonction qui compte le nombre de nombres premiers", PhD thesis, 1998. In French, but the mathematics is readable and highly recommended reading if you're interesting in prime number bounds. L<http://www.unilim.fr/laco/theses/1998/T1998_01.html>
3024
3025		=item Gabriel Mincu, "An Asymptotic Expansion", Journal of Inequalities in Pure and Applied Mathematics, v4, n2, 2003. A very readable account of Cipolla's 1902 nth prime approximation. L<http://www.emis.de/journals/JIPAM/images/153_02_JIPAM/153_02.pdf>
3026
3027		=item David M. Smith, "Multiple-Precision Exponential Integral and Related Functions".
3028
3029		=item Vincent Pegoraro and Philipp Slusallek, "On the Evaluation of the Complex-Valued Exponential Integral".
3030
3031		=item William H. Press et al., "Numerical Recipes", 3rd edition.
3032
3033		=item W. J. Cody and Henry C. Thacher, Jr., "Rational Chevyshev Approximations for the Exponential Integral E_1(x)".
3034
3035		=item W. J. Cody, K. E. Hillstrom, and Henry C. Thacher Jr., "Chebyshev Approximations for the Riemann Zeta Function", Mathematics of Computation, v25, n115, pp 537-547, July 1971.
3036
3037		=item Ueli M. Maurer, "Fast Generation of Prime Numbers and Secure Public-Key Cryptographic Parameters", 1995. L<http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.26.2151>
3038
3039		=item Pierre-Alain Fouque and Mehdi Tibouchi, "Close to Uniform Prime Number Generation With Fewer Random Bits", 2011. L<http://eprint.iacr.org/2011/481>
3040
3041		=item Douglas A. Stoll and Patrick Demichel , "The impact of ζ(s) complex zeros on π(x) for x E<lt> 10^{10^{13}}", Mathematics of Computation, v80, n276, pp 2381-2394, October 2011. L<http://www.ams.org/journals/mcom/2011-80-276/S0025-5718-2011-02477-4/home.html>
3042
3043		=item L<OEIS: Primorial\|http://oeis.org/wiki/Primorial>.
	3185	=item *
	3186
	3187	Pierre Dusart, "Estimates of Some Functions Over Primes without R.H.", preprint, 2010. L<http://arxiv.org/abs/1002.0442/>
	3188
	3189	=item *
	3190
	3191	Pierre Dusart, "Autour de la fonction qui compte le nombre de nombres premiers", PhD thesis, 1998. In French, but the mathematics is readable and highly recommended reading if you're interesting in prime number bounds. L<http://www.unilim.fr/laco/theses/1998/T1998_01.html>
	3192
	3193	=item *
	3194
	3195	Gabriel Mincu, "An Asymptotic Expansion", Journal of Inequalities in Pure and Applied Mathematics, v4, n2, 2003. A very readable account of Cipolla's 1902 nth prime approximation. L<http://www.emis.de/journals/JIPAM/images/153_02_JIPAM/153_02.pdf>
	3196
	3197	=item *
	3198
	3199	David M. Smith, "Multiple-Precision Exponential Integral and Related Functions".
	3200
	3201	=item *
	3202
	3203	Vincent Pegoraro and Philipp Slusallek, "On the Evaluation of the Complex-Valued Exponential Integral".
	3204
	3205	=item *
	3206
	3207	William H. Press et al., "Numerical Recipes", 3rd edition.
	3208
	3209	=item *
	3210
	3211	W. J. Cody and Henry C. Thacher, Jr., "Rational Chevyshev Approximations for the Exponential Integral E_1(x)".
	3212
	3213	=item *
	3214
	3215	W. J. Cody, K. E. Hillstrom, and Henry C. Thacher Jr., "Chebyshev Approximations for the Riemann Zeta Function", Mathematics of Computation, v25, n115, pp 537-547, July 1971.
	3216
	3217	=item *
	3218
	3219	Ueli M. Maurer, "Fast Generation of Prime Numbers and Secure Public-Key Cryptographic Parameters", 1995. L<http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.26.2151>
	3220
	3221	=item *
	3222
	3223	Pierre-Alain Fouque and Mehdi Tibouchi, "Close to Uniform Prime Number Generation With Fewer Random Bits", 2011. L<http://eprint.iacr.org/2011/481>
	3224
	3225	=item *
	3226
	3227	Douglas A. Stoll and Patrick Demichel , "The impact of ζ(s) complex zeros on π(x) for x E<lt> 10^{10^{13}}", Mathematics of Computation, v80, n276, pp 2381-2394, October 2011. L<http://www.ams.org/journals/mcom/2011-80-276/S0025-5718-2011-02477-4/home.html>
	3228
	3229	=item *
	3230
	3231	L<OEIS: Primorial\|http://oeis.org/wiki/Primorial>.
3044	3232
3045	3233	=back
3046	3234

+10

-4

ptypes.h less more

0	0	#ifndef MPU_PTYPES_H
1	1	#define MPU_PTYPES_H
2	2
3		#ifndef _MSC_VER
4		#define __STDC_LIMIT_MACROS
5		#include <stdint.h>
6		#else
	3	#ifdef _MSC_VER
7	4	/* No stdint.h for MS C, so we lose the chance to possibly optimize
8	5	* some operations on 64-bit machines running a 32-bit Perl. It's probably
9	6	* a rare enough case that we don't need to be too concerned. If we do want,

13	10	* Thanks to Sisyphus for bringing the MSC issue to my attention (and even
14	11	* submitting a working patch!).
15	12	*/
	13	#elif defined(__sun) \|\| defined(__sun__)
	14	/* stdint.h is only in Solaris 10+. */
	15	#if defined(__SunOS_5_10) \|\| defined(__SunOS_5_11) \|\| defined(__SunOS_5_12)
	16	#define __STDC_LIMIT_MACROS
	17	#include <stdint.h>
	18	#endif
	19	#else
	20	#define __STDC_LIMIT_MACROS
	21	#include <stdint.h>
16	22	#endif
17	23
18	24	#include "EXTERN.h"

+14

-1

t/13-primecount.t less more

4	4	use Test::More;
5	5	use Math::Prime::Util qw/prime_count prime_count_lower prime_count_upper prime_count_approx/;
6	6
	7	my $isxs = Math::Prime::Util::prime_get_config->{'xs'};
7	8	my $use64 = Math::Prime::Util::prime_get_config->{'maxbits'} > 32;
8	9	my $extra = defined $ENV{RELEASE_TESTING} && $ENV{RELEASE_TESTING};
9	10

81	82	+ scalar(keys %pivals_small)
82	83	+ $use64 * 3 * scalar(keys %pivals64)
83	84	+ scalar(keys %intervals)
84		+ 1;
	85	+ 1
	86	+ 6; # prime count specific methods
85	87
86	88	ok( eval { prime_count(13); 1; }, "prime_count in void context");
87	89

147	149	# 109726486, // prime count 2^32 interval starting at 10^17
148	150	# 103626726, // prime count 2^32 interval starting at 10^18
149	151	# 98169972}; // prime count 2^32 interval starting at 10^19
	152
	153	# Make sure each specific algorithm isn't broken.
	154	SKIP: {
	155	skip "Not XS -- skipping direct primecount tests", 4 unless $isxs;
	156	is(Math::Prime::Util::_XS_lehmer_pi (3456789), 247352, "XS Lehmer count");
	157	is(Math::Prime::Util::_XS_meissel_pi (3456789), 247352, "XS Meissel count");
	158	is(Math::Prime::Util::_XS_legendre_pi(3456789), 247352, "XS Legendre count");
	159	is(Math::Prime::Util::_XS_prime_count(3456789), 247352, "XS sieve count");
	160	}
	161	is(Math::Prime::Util::PP::_lehmer_pi (3456789), 247352, "PP Lehmer count");
	162	is(Math::Prime::Util::PP::_sieve_prime_count(3456789), 247352, "PP sieve count");

-1

t/22-aks-prime.t less more

34	34	#diag "Unfortunately these tests are very slow.";
35	35
36	36	SKIP: {
37		skip "Skipping PP AKS on 32-bit -- just too slow.", 1 if $ispp && !$use64;
	37	# If we're pure Perl, then this is definitely too slow.
	38	# Arguably we should check to see if they have the GMP code.
	39	skip "Skipping PP AKS on PP -- just too slow.", 1 if $ispp;
	40	# If we have 64-bit available in the compiler (e.g. uint64_t), this can
	41	# still be quite fast. However for pretty much everyone else, this is
	42	# just far too slow for running in a test suite.
	43	skip "Skipping PP AKS on 32-bit -- just too slow.", 1 if !$use64;
38	44	# The first number that makes it past the sqrt test to actually run.
39	45	is( is_aks_prime(69197), 1, "is_aks_prime(69197) is true" );
40	46	}

-3

t/50-factoring.t less more

29	29	1363 989 779 629 403
30	30	547308031
31	31	808 2727 12625 34643 134431 221897 496213 692759 1228867
32		2463289 3008891 5115953 6961021 8030207 10486123
	32	2231139 2463289 3008891 5115953 6961021 8030207 10486123
33	33	10893343 12327779 701737021
	34	549900
34	35	/;
35	36
36	37	my @testn64 = qw/37607912018 346065536839 600851475143

71	72	0 => [],
72	73	);
73	74
74		plan tests => (2 * scalar @testn) + scalar(keys %all_factors) + 6*7;
	75	plan tests => (2 * scalar @testn) + scalar(keys %all_factors) + 7*8;
75	76
76	77	foreach my $n (@testn) {
77	78	my @f = factor($n);

99	100	extra_factor_test("fermat_factor", sub {Math::Prime::Util::fermat_factor(shift)});
100	101	extra_factor_test("holf_factor", sub {Math::Prime::Util::holf_factor(shift)});
101	102	extra_factor_test("squfof_factor", sub {Math::Prime::Util::squfof_factor(shift)});
	103	extra_factor_test("rsqufof_factor", sub {Math::Prime::Util::rsqufof_factor(shift)});
102	104	extra_factor_test("pbrent_factor", sub {Math::Prime::Util::pbrent_factor(shift)});
103	105	extra_factor_test("prho_factor", sub {Math::Prime::Util::prho_factor(shift)});
104	106	extra_factor_test("pminus1_factor",sub {Math::Prime::Util::pminus1_factor(shift)});

110	112	is_deeply( [ sort {$a<=>$b} $fsub->(1) ], [1], "$fname(1)" );
111	113	is_deeply( [ sort {$a<=>$b} $fsub->(4) ], [2, 2], "$fname(4)" );
112	114	is_deeply( [ sort {$a<=>$b} $fsub->(9) ], [3, 3], "$fname(9)" );
113		is_deeply( [ sort {$a<=>$b} $fsub->(25) ], [5, 5], "$fname(9)" );
	115	is_deeply( [ sort {$a<=>$b} $fsub->(25) ], [5, 5], "$fname(25)" );
114	116	is_deeply( [ sort {$a<=>$b} $fsub->(175) ], [5, 5, 7], "$fname(175)" );
115	117	is_deeply( [ sort {$a<=>$b} $fsub->(403) ], [13, 31], "$fname(403)" );
	118	is_deeply( [ sort {$a<=>$b} $fsub->(549900) ], [2,2,3,3,5,5,13,47], "$fname(549900)" );
116	119	}
117	120

+31

-0

t/70-rt-bignum.t less more

	0	#!/usr/bin/env perl
	1	use strict;
	2	use warnings;
	3
	4	# I found these issues when doing some testing of is_provable_prime. When
	5	# bignum is loaded, we get some strange behavior. There are two fixes for
	6	# it in the code:
	7	# 1) make sure every divide and bdiv is coerced back to an integer.
	8	# 2) turn off upgrade in input validation.
	9	# The second method in theory is all that is needed.
	10
	11	use Math::Prime::Util qw/:all/;
	12	use bignum;
	13
	14	use Test::More tests => 1;
	15
	16	if ($] < 5.008) {
	17	diag "A prototype warning was expected with old, old Perl";
	18	}
	19
	20	my $n = 100199294509778143137521762187425301691197073534078445671945250753109628678272;
	21	# 2 2 2 2 2 2 2 3 7 509 277772399 263650456338779643073784729209358382310353002641378210462709359
	22
	23	my @partial_factor = Math::Prime::Util::PP::prho_factor(100199294509778143137521762187425301691197073534078445671945250753109628678272, 5);
	24
	25	is_deeply( \@partial_factor,
	26	[2,2,2,2,2,2,2,3,7,37276523255125797298185179385202865212498911284999421752955822452793760669],
	27	"PP prho factors correctly with 'use bignum'" );
	28
	29	# The same thing happens in random primes, PP holf factoring,
	30	# PP is_provable_primes, and possibly elsewhere

+74

-1

t/80-pp.t less more

54	54	61 => [ qw/217 341 1261 2701 3661 6541 6697 7613 13213 16213 22177 23653 23959 31417 50117 61777 63139 67721 76301 77421 79381 80041/ ],
55	55	73 => [ qw/205 259 533 1441 1921 2665 3439 5257 15457 23281 24617 26797 27787 28939 34219 39481 44671 45629 64681 67069 76429 79501 93521/ ],
56	56	);
	57	# Test a pseudoprime larger than 2^32.
	58	push @{$pseudoprimes{2}}, 75792980677 if $use64;
57	59	my $num_pseudoprimes = 0;
58	60	foreach my $ppref (values %pseudoprimes) {
59	61	push @composites, @$ppref;

220	222	2*scalar(keys %pivals_small) + scalar(keys %nthprimes_small) +
221	223	4 + scalar(keys %pseudoprimes) +
222	224	scalar(keys %eivals) + scalar(keys %livals) + scalar(keys %rvals) +
223		1 + 1 +
	225	1 + 1 + # factor
	226	8 + 4*3 + # factoring subs
	227	10 + # AKS
224	228	0;
225	229
226	230	use Math::Prime::Util qw/primes prime_count_approx prime_count_lower/;
	231	use Math::BigInt try => 'GMP';
227	232	require_ok 'Math::Prime::Util::PP';
228	233	# This function skips some setup
229	234	undef *primes;

237	242	*prev_prime = \&Math::Prime::Util::PP::prev_prime;
238	243
239	244	*miller_rabin = \&Math::Prime::Util::PP::miller_rabin;
	245	*is_aks_prime = \&Math::Prime::Util::PP::is_aks_prime;
240	246
241	247	*factor = \&Math::Prime::Util::PP::factor;
242	248

395	401	}
396	402	}
397	403	is_deeply( \@gotfactor, \@expfactor, "test factoring for $ntests composites");
	404	}
	405
	406	# The PP factor code does small trials, then loops doing 64k rounds of HOLF
	407	# if the composite is less than a half word, followed by 64k rounds each of
	408	# prho with a = {3,5,7,11,13}. Most numbers are handled by these. The ones
	409	# that aren't end up being too slow for us to put in a test. So we'll try
	410	# running the various factoring methods manually.
	411	{
	412	is_deeply( [ sort {$a<=>$b} Math::Prime::Util::PP::holf_factor(403) ],
	413	[ 13, 31 ],
	414	"holf(403)" );
	415	is_deeply( [ sort {$a<=>$b} Math::Prime::Util::PP::prho_factor(403) ],
	416	[ 13, 31 ],
	417	"prho(403)" );
	418	is_deeply( [ sort {$a<=>$b} Math::Prime::Util::PP::pbrent_factor(403) ],
	419	[ 13, 31 ],
	420	"pbrent(403)" );
	421	is_deeply( [ sort {$a<=>$b} Math::Prime::Util::PP::pminus1_factor(403) ],
	422	[ 13, 31 ],
	423	"pminus1(403)" );
	424	is_deeply( [ sort {$a<=>$b} Math::Prime::Util::PP::prho_factor(851981) ],
	425	[ 13, 65537 ],
	426	"prho(851981)" );
	427	is_deeply( [ sort {$a<=>$b} Math::Prime::Util::PP::pbrent_factor(851981) ],
	428	[ 13, 65537 ],
	429	"pbrent(851981)" );
	430	my $n64 = $use64 ? 55834573561 : Math::BigInt->new("55834573561");
	431	is_deeply( [ sort {$a<=>$b} Math::Prime::Util::PP::prho_factor($n64) ],
	432	[ 13, 4294967197 ],
	433	"prho(55834573561)" );
	434	is_deeply( [ sort {$a<=>$b} Math::Prime::Util::PP::pbrent_factor($n64) ],
	435	[ 13, 4294967197 ],
	436	"pbrent(55834573561)" );
	437	# 1013 4294967197 4294967291
	438	my $nbig = Math::BigInt->new("18686551294184381720251");
	439	my @nfac;
	440	@nfac = sort {$a<=>$b} Math::Prime::Util::PP::holf_factor($nbig);
	441	is(scalar @nfac, 2, "holf finds a factor of 18686551294184381720251");
	442	is($nfac[0] * $nfac[1], $nbig, "holf found a correct factor");
	443	ok($nfac[0] != 1 && $nfac[1] != 1, "holf didn't return a degenerate factor");
	444	@nfac = sort {$a<=>$b} Math::Prime::Util::PP::prho_factor($nbig);
	445	is(scalar @nfac, 2, "prho finds a factor of 18686551294184381720251");
	446	is($nfac[0] * $nfac[1], $nbig, "prho found a correct factor");
	447	ok($nfac[0] != 1 && $nfac[1] != 1, "prho didn't return a degenerate factor");
	448	@nfac = sort {$a<=>$b} Math::Prime::Util::PP::pbrent_factor($nbig);
	449	is(scalar @nfac, 2, "pbrent finds a factor of 18686551294184381720251");
	450	is($nfac[0] * $nfac[1], $nbig, "pbrent found a correct factor");
	451	ok($nfac[0] != 1 && $nfac[1] != 1, "pbrent didn't return a degenerate factor");
	452	@nfac = sort {$a<=>$b} Math::Prime::Util::PP::pminus1_factor($nbig);
	453	is(scalar @nfac, 2, "pminus1 finds a factor of 18686551294184381720251");
	454	is($nfac[0] * $nfac[1], $nbig, "pminus1 found a correct factor");
	455	ok($nfac[0] != 1 && $nfac[1] != 1, "pminus1 didn't return a degenerate factor");
	456	}
	457
	458	##### AKS primality test. Be very careful with performance.
	459	is( is_aks_prime(1), 0, "AKS: 1 is composite (less than 2)" );
	460	is( is_aks_prime(2), 1, "AKS: 2 is prime" );
	461	is( is_aks_prime(3), 1, "AKS: 3 is prime" );
	462	is( is_aks_prime(4), 0, "AKS: 4 is composite" );
	463	is( is_aks_prime(64), 0, "AKS: 64 is composite (perfect power)" );
	464	is( is_aks_prime(65), 0, "AKS: 65 is composite (caught in trial)" );
	465	is( is_aks_prime(23), 1, "AKS: 23 is prime (r >= n)" );
	466	is( is_aks_prime(101), 1, "AKS: 101 is prime (passed anr test)" );
	467	is( is_aks_prime(70747), 0, "AKS: 70747 is composite (n mod r)" );
	468	SKIP: {
	469	skip "Skipping PP AKS test on 32-bit machine", 1 unless $use64 \|\| $extra;
	470	is( is_aks_prime(74513), 0, "AKS: 74513 is composite (failed anr test)" );
398	471	}
399	472
400	473	###############################################################################

+23

-4

t/81-bignum.t less more

24	24	my @composites = qw/
25	25	777777777777777777777777 877777777777777777777777
26	26	87777777777777777777777795475 890745785790123461234805903467891234681234
	27	/;
	28
	29	# Primes where n-1 is easy to factor, so we finish quickly.
	30	my @proveprimes = qw/
	31	65635624165761929287 1162566711635022452267983 77123077103005189615466924501
	32	3991617775553178702574451996736229 273952953553395851092382714516720001799
27	33	/;
28	34
29	35	# pseudoprimes to various small prime bases

58	64	);
59	65
60	66	plan tests => 0
61		+ 2*(@primes + @composites)
	67	+ 2*(@primes + @composites + @proveprimes)
62	68	+ 1 # primes
63	69	+ 2 # next/prev prime
64	70	+ 1 # primecount large base small range

68	74	+ scalar(keys %factors)
69	75	+ scalar(keys %allfactors)
70	76	+ 2 # moebius, euler_phi
71		+ 12 # random primes
	77	+ 15 # random primes
72	78	+ 0;
73	79
74	80	# Using GMP makes these tests run about 2x faster on some machines

95	101	prime_count
96	102	nth_prime
97	103	is_prime
	104	is_provable_prime
98	105	next_prime
99	106	prev_prime
100	107	is_strong_pseudoprime
101	108	random_prime
102	109	random_ndigit_prime
103	110	random_nbit_prime
	111	random_strong_prime
104	112	random_maurer_prime
105	113	/;
106	114	# TODO: is_strong_lucas_pseudoprime

108	116	# LogarithmicIntegral
109	117	# RiemannR
110	118
	119	my $bignumver = $bigint::VERSION;
	120	my $bigintver = $Math::BigInt::VERSION;
111	121	my $bigintlib = Math::BigInt->config()->{lib};
112	122	$bigintlib =~ s/^Math::BigInt:://;
113	123	my $mpugmpver = Math::Prime::Util::prime_get_config->{gmp}
114	124	? $Math::Prime::Util::GMP::VERSION : "<none>";
115		diag "BigInt library: $bigintlib, MPU::GMP $mpugmpver\n";
	125	diag "BigInt $bignumver/$bigintver, lib: $bigintlib. MPU::GMP $mpugmpver\n";
116	126
117	127
118	128	###############################################################################

124	134	foreach my $n (@composites) {
125	135	ok( !is_prime($n), "$n is not prime" );
126	136	ok( !is_prob_prime($n), "$n is not probably prime");
	137	}
	138	foreach my $n (@proveprimes) {
	139	ok( is_prime($n), "$n is prime" );
	140	ok( is_provable_prime($n), "$n is provably prime" );
127	141	}
128	142
129	143	###############################################################################

205	219	cmp_ok( $randprime, '<', 2**80, "random 80-bit prime isn't too big");
206	220	ok( is_prime($randprime), "random 80-bit prime is prime");
207	221
	222	$randprime = random_strong_prime(256);
	223	cmp_ok( $randprime, '>', 2**255, "random 256-bit strong prime isn't too small");
	224	cmp_ok( $randprime, '<', 2**256, "random 256-bit strong prime isn't too big");
	225	ok( is_prime($randprime), "random 80-bit strong prime is prime");
	226
208	227	SKIP: {
209	228	skip "Your 64-bit Perl is broken, skipping maurer prime", 3 if $broken64;
210	229	$randprime = random_maurer_prime(80);

229	248	prime_set_config(assume_rh=>1);
230	249	my $pclo_rh = prime_count_lower($n);
231	250	my $pcup_rh = prime_count_upper($n);
232		prime_set_config(assume_rh=>0);
	251	prime_set_config(assume_rh => undef);
233	252
234	253	#diag "lower: " . $pclo->bstr() . " " . ($pcap-$pclo)->bstr;
235	254	#diag "rh lower: " . $pclo_rh->bstr() . " " . ($pcap-$pclo_rh)->bstr;