[libmath-prime-util-perl] 34/50: Add and enhance examples, add bin/factor.pl
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:45:37 UTC 2015
This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.13
in repository libmath-prime-util-perl.
commit a23c8b742bc6dc82e64f5487fb0bf1834fe2e145
Author: Dana Jacobsen <dana at acm.org>
Date: Thu Nov 8 04:21:43 2012 -0800
Add and enhance examples, add bin/factor.pl
---
.gitignore | 11 +++
Changes | 2 +-
MANIFEST | 1 +
Makefile.PL | 2 +-
bin/factor.pl | 50 +++++++++++++
bin/primes.pl | 11 ++-
examples/bench-factor-extra.pl | 45 +++++++++---
examples/bench-factor-semiprime.pl | 51 +++++++++----
examples/bench-factor.pl | 39 ++++++++--
examples/test-bpsw.pl | 145 ++++++++++++++++++++++++++++++++++++
examples/test-euler-pari.pl | 37 ++++++++++
examples/test-factor-gnufactor.pl | 147 +++++++++++++++++++++++++++++++++++++
examples/test-factor-mpxs.pl | 23 +++++-
examples/test-factor-yafu.pl | 39 ++++++++--
14 files changed, 557 insertions(+), 46 deletions(-)
diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..dfc9838
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,11 @@
+Math-Prime-Util-*.tar.gz
+Makefile
+Makefile.old
+MYMETA.*
+Util.bs
+XS.[co]
+cache.o
+factor.o
+sieve.o
+util.o
+blib*
diff --git a/Changes b/Changes
index d85d848..cf18037 100644
--- a/Changes
+++ b/Changes
@@ -2,7 +2,7 @@ Revision history for Perl extension Math::Prime::Util.
0.12 2 November 2012
- - Add bin/primes.pl
+ - Add bin/primes.pl and bin/factor.pl
- Add functions:
primorial product of primes <= n
diff --git a/MANIFEST b/MANIFEST
index d2db4c1..eef1b35 100644
--- a/MANIFEST
+++ b/MANIFEST
@@ -42,6 +42,7 @@ examples/test-pcapprox.pl
examples/sophie_germain.pl
examples/twin_primes.pl
bin/primes.pl
+bin/factor.pl
t/01-load.t
t/02-can.t
t/03-init.t
diff --git a/Makefile.PL b/Makefile.PL
index a8a29c2..1592bf7 100644
--- a/Makefile.PL
+++ b/Makefile.PL
@@ -14,7 +14,7 @@ WriteMakefile1(
'XS.o',
LIBS => ['-lm'],
- EXE_FILES => ['bin/primes.pl'],
+ EXE_FILES => ['bin/primes.pl', 'bin/factor.pl'],
BUILD_REQUIRES=>{
'Test::More' => '0.45',
diff --git a/bin/factor.pl b/bin/factor.pl
new file mode 100755
index 0000000..2b585a1
--- /dev/null
+++ b/bin/factor.pl
@@ -0,0 +1,50 @@
+#!/usr/bin/env perl
+use strict;
+use warnings;
+use Getopt::Long;
+use bigint try => 'GMP';
+use Math::Prime::Util qw/factor/;
+$| = 1;
+no bigint;
+
+my %opts;
+GetOptions(\%opts,
+ 'version', # turn off MPU::GMP for debugging
+ 'help',
+ ) || die_usage();
+if (exists $opts{'version'}) {
+ my $version_str =
+ "factor.pl version 1.0 using Math::Prime::Util $Math::Prime::Util::VERSION";
+ $version_str .= " and MPU::GMP $Math::Prime::Util::GMP::VERSION"
+ if Math::Prime::Util::prime_get_config->{'gmp'};
+ $version_str .= "\nWritten by Dana Jacobsen.\n";
+ die "$version_str";
+}
+die_usage() if exists $opts{'help'};
+
+if (@ARGV) {
+ foreach my $n (@ARGV) {
+ print "$n: ", join(" ", factor($n)), "\n";
+ }
+} else {
+ while (<>) {
+ chomp;
+ foreach my $n (split / /) {
+ print "$n: ", join(" ", factor($n)), "\n";
+ }
+ }
+}
+sub die_usage {
+ die <<EOU;
+Usage: $0 [options] [number] ...
+
+Print the prime factors of each positive integer given on the command line,
+or reads numbers from standard input if called without arguments.
+
+ --help displays this help message
+ --version displays the version information
+
+Part of the Math::Prime::Util $Math::Prime::Util::VERSION package, wrapping
+the factor() function. See 'man Math::Prime::Util' for more information.
+EOU
+}
diff --git a/bin/primes.pl b/bin/primes.pl
index 091a71a..f6c3560 100755
--- a/bin/primes.pl
+++ b/bin/primes.pl
@@ -93,10 +93,19 @@ GetOptions(\%opts,
'euclid|A018239',
'provable',
'nompugmp', # turn off MPU::GMP for debugging
+ 'version',
'help',
) || die_usage();
-die_usage() if exists $opts{'help'};
Math::Prime::Util::prime_set_config(gmp=>0) if exists $opts{'nompugmp'};
+if (exists $opts{'version'}) {
+ my $version_str =
+ "primes.pl version 1.2 using Math::Prime::Util $Math::Prime::Util::VERSION";
+ $version_str .= " and MPU::GMP $Math::Prime::Util::GMP::VERSION"
+ if Math::Prime::Util::prime_get_config->{'gmp'};
+ $version_str .= "\nWritten by Dana Jacobsen.\n";
+ die "$version_str";
+}
+die_usage() if exists $opts{'help'};
# Get the start and end values. Verify they're positive integers.
die_usage() unless @ARGV == 2;
diff --git a/examples/bench-factor-extra.pl b/examples/bench-factor-extra.pl
index 8a280da..98d3342 100755
--- a/examples/bench-factor-extra.pl
+++ b/examples/bench-factor-extra.pl
@@ -1,17 +1,37 @@
#!/usr/bin/env perl
use strict;
use warnings;
-use Math::Prime::Util;
+use Math::Prime::Util qw/-nobigint/;
use Benchmark qw/:all/;
use List::Util qw/min max/;
+use Config;
my $count = shift || -2;
my $is64bit = (~0 > 4294967295);
my $maxdigits = ($is64bit) ? 20 : 10; # Noting the range is limited for max.
+my $rgen = sub {
+ my $range = shift;
+ return 0 if $range <= 0;
+ my $rbits = 0; { my $t = $range; while ($t) { $rbits++; $t >>= 1; } }
+ while (1) {
+ my $rbitsleft = $rbits;
+ my $U = 0;
+ while ($rbitsleft > 0) {
+ my $usebits = ($rbitsleft > $Config{randbits}) ? $Config{randbits} : $rbitsleft;
+ $U = ($U << $usebits) + int(rand(1 << $usebits));
+ $rbitsleft -= $usebits;
+ }
+ return $U if $U <= $range;
+ }
+};
+
srand(29);
my $rounds = 400;
my $sqrounds = 256*1024;
-my $hrounds = 2000;
+my $rsqrounds = 32*1024;
+my $p1smooth = 1000;
+my $hrounds = 10000;
+my $num_nums = 1000;
test_at_digits($_) for ( 3 .. $maxdigits );
@@ -21,28 +41,29 @@ sub test_at_digits {
die "Digits has to be >= 1" unless $digits >= 1;
die "Digits has to be <= $maxdigits" if $digits > $maxdigits;
- my @nums = genrand($digits, 1000);
- #my @nums = gensemi($digits, 1000, 23);
+ my @nums = genrand($digits, $num_nums);
+ #my @nums = gensemi($digits, $num_nums, 23);
my $min_num = min @nums;
my $max_num = max @nums;
# Determine success rates
my %nfactored;
my $tfac = 0;
+ # Did we find any non-trivial factors?
my $calc_nfacs = sub { ((scalar grep { $_ > 5 } @_) > 1) ? 1 : 0 };
for (@nums) {
$tfac += $calc_nfacs->(Math::Prime::Util::factor($_));
$nfactored{'prho'} += $calc_nfacs->(Math::Prime::Util::prho_factor($_, $rounds));
$nfactored{'pbrent'} += $calc_nfacs->(Math::Prime::Util::pbrent_factor($_, $rounds));
- $nfactored{'pminus1'} += $calc_nfacs->(Math::Prime::Util::pminus1_factor($_, $rounds));
+ $nfactored{'pminus1'} += $calc_nfacs->(Math::Prime::Util::pminus1_factor($_, $p1smooth));
$nfactored{'squfof'} += $calc_nfacs->(Math::Prime::Util::squfof_factor($_, $sqrounds));
- $nfactored{'rsqufof'} += $calc_nfacs->(Math::Prime::Util::rsqufof_factor($_));
+ $nfactored{'rsqufof'} += $calc_nfacs->(Math::Prime::Util::rsqufof_factor($_, $rsqrounds));
#$nfactored{'trial'} += $calc_nfacs->(Math::Prime::Util::trial_factor($_));
#$nfactored{'fermat'} += $calc_nfacs->(Math::Prime::Util::fermat_factor($_, $rounds));
$nfactored{'holf'} += $calc_nfacs->(Math::Prime::Util::holf_factor($_, $hrounds));
}
- print "factoring 1000 random $digits-digit numbers ($min_num - $max_num)\n";
+ print "factoring $num_nums random $digits-digit numbers ($min_num - $max_num)\n";
print "Factorizations: ",
join(", ", map { sprintf "%s %4.1f%%", $_, 100*$nfactored{$_}/$tfac }
grep { $_ ne 'fermat' }
@@ -60,7 +81,8 @@ sub test_at_digits {
"trial" => sub { Math::Prime::Util::trial_factor($_) for @nums },
};
delete $lref->{'fermat'} if $digits >= 9;
- #delete $lref->{'holf'} if $digits >= 9;
+ delete $lref->{'holf'} if $digits >= 17;
+ delete $lref->{'trial'} if $digits >= 15;
cmpthese($count, $lref);
print "\n";
}
@@ -73,7 +95,7 @@ sub genrand {
my $base = ($digits == 1) ? 0 : int(10 ** ($digits-1));
my $max = int(10 ** $digits);
$max = ~0 if $max > ~0;
- my @nums = map { $base+int(rand($max-$base)) } (1 .. $num);
+ my @nums = map { $base + $rgen->($max-$base) } (1 .. $num);
return @nums;
}
@@ -91,9 +113,8 @@ sub gensemi {
my @factors;
my $n;
while (1) {
- $n = $base + int(rand($max-$base));
- $n += 1 if ($n%2) == 0;
- $n += 3 if ($n%3) == 0;
+ $n = $base + $rgen->($max-$base);
+ $n += (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)[$n%30];
@factors = Math::Prime::Util::factor($n);
next if scalar @factors != 2;
next if $factors[0] < $smallest_factor;
diff --git a/examples/bench-factor-semiprime.pl b/examples/bench-factor-semiprime.pl
index 4e289ad..b9973f0 100755
--- a/examples/bench-factor-semiprime.pl
+++ b/examples/bench-factor-semiprime.pl
@@ -4,23 +4,42 @@ use warnings;
$| = 1; # fast pipes
srand(377);
-use Math::Prime::Util qw/factor/;
+use Math::Prime::Util qw/factor -nobigint/;
use Math::Factor::XS qw/prime_factors/;
use Math::Pari qw/factorint/;
use Benchmark qw/:all/;
use Data::Dumper;
+use Config;
my $digits = shift || 15;
-my $count = shift || -2;
+my $count = shift || -3;
+
+my $rgen = sub {
+ my $range = shift;
+ return 0 if $range <= 0;
+ my $rbits = 0; { my $t = $range; while ($t) { $rbits++; $t >>= 1; } }
+ while (1) {
+ my $rbitsleft = $rbits;
+ my $U = 0;
+ while ($rbitsleft > 0) {
+ my $usebits = ($rbitsleft > $Config{randbits}) ? $Config{randbits} : $rbitsleft;
+ $U = ($U << $usebits) + int(rand(1 << $usebits));
+ $rbitsleft -= $usebits;
+ }
+ return $U if $U <= $range;
+ }
+};
my @min_factors_by_digit = (2,2,3,3,5,11,17,47,97);
my $smallest_factor_allowed = $min_factors_by_digit[$digits];
$smallest_factor_allowed = $min_factors_by_digit[-1] unless defined $smallest_factor_allowed;
-my $numprimes = 100;
+my $numprimes = 200;
die "Digits has to be >= 2" unless $digits >= 2;
die "Digits has to be <= 10" if (~0 == 4294967295) && ($digits > 10);
die "Digits has to be <= 19" if $digits > 19;
+my $skip_mfxs = ($digits > 17);
+
# Construct some semiprimes of the appropriate number of digits
# There are much cleverer ways of doing this, using randomly selected
# nth_primes, and so on, but this works well until we get lots of digits.
@@ -32,10 +51,9 @@ foreach my $i ( 1 .. $numprimes ) {
my @factors;
my $n;
while (1) {
- $n = $base + int(rand($add));
+ $n = $base + $rgen->($add);
next if $n > (~0 - 4);
- $n += 1 if ($n%2) == 0;
- $n += 3 if ($n%3) == 0;
+ $n += (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)[$n%30];
@factors = factor($n);
next if scalar @factors != 2;
next if $factors[0] < $smallest_factor_allowed;
@@ -55,12 +73,16 @@ foreach my $sp (@semiprimes) {
die "wrong for $sp\n" unless ($#factors == 1) && ($factors[0] * $factors[1]) == $sp;
}
print "OK\n";
-print "Verifying Math::Factor::XS $Math::Factor::XS::VERSION ...";
-foreach my $sp (@semiprimes) {
- my @factors = prime_factors($sp);
- die "wrong for $sp\n" unless ($#factors == 1) && ($factors[0] * $factors[1]) == $sp;
+if (!$skip_mfxs) {
+ print "Verifying Math::Factor::XS $Math::Factor::XS::VERSION ...";
+ foreach my $sp (@semiprimes) {
+ my @factors = prime_factors($sp);
+ die "wrong for $sp\n" unless ($#factors == 1) && ($factors[0] * $factors[1]) == $sp;
+ }
+ print "OK\n";
+} else {
+ print "Math::Factor::XS is too slow for $digits digits. Skipping.\n";
}
-print "OK\n";
print "Verifying Math::Pari $Math::Pari::VERSION ...";
foreach my $sp (@semiprimes) {
my @factors;
@@ -70,8 +92,11 @@ foreach my $sp (@semiprimes) {
}
print "OK\n";
-cmpthese($count,{
+my %compare = (
'MPU' => sub { factor($_) for @semiprimes; },
'MFXS' => sub { prime_factors($_) for @semiprimes; },
'Pari' => sub { foreach my $n (@semiprimes) { my @factors; my ($pn,$pc) = @{factorint($n)}; push @factors, (int($pn->[$_])) x $pc->[$_] for (0 .. $#{$pn}); } }
-});
+);
+delete $compare{'MFXS'} if $skip_mfxs;
+
+cmpthese($count, \%compare);
diff --git a/examples/bench-factor.pl b/examples/bench-factor.pl
index 04ea83a..e9e5536 100755
--- a/examples/bench-factor.pl
+++ b/examples/bench-factor.pl
@@ -1,19 +1,42 @@
#!/usr/bin/env perl
use strict;
use warnings;
+
+# One of the things this test shows is the impact of the '-nobigint' option.
+# "MPU" results are the timings without the '-nobigint' option
+# "MPUxs" results are the timings _with_ the '-nobigint' option
use Math::Prime::Util qw/factor/;
+
+# Compare to Math::Factor::XS, which uses trial division.
use Math::Factor::XS qw/prime_factors/;
+
use Benchmark qw/:all/;
use List::Util qw/min max reduce/;
+use Config;
my $count = shift || -2;
my $is64bit = (~0 > 4294967295);
my $maxdigits = ($is64bit) ? 20 : 10; # Noting the range is limited for max.
my $semiprimes = 0;
-my $howmany = 10000;
-
+my $howmany = 1000;
+
+my $rgen = sub {
+ my $range = shift;
+ return 0 if $range <= 0;
+ my $rbits = 0; { my $t = $range; while ($t) { $rbits++; $t >>= 1; } }
+ while (1) {
+ my $rbitsleft = $rbits;
+ my $U = 0;
+ while ($rbitsleft > 0) {
+ my $usebits = ($rbitsleft > $Config{randbits}) ? $Config{randbits} : $rbitsleft;
+ $U = ($U << $usebits) + int(rand(1 << $usebits));
+ $rbitsleft -= $usebits;
+ }
+ return $U if $U <= $range;
+ }
+};
srand(29);
-for my $d ( 3 .. $maxdigits ) {
+for my $d ( 17 .. $maxdigits ) {
print "Factor $howmany $d-digit numbers\n";
test_at_digits($d, $howmany);
}
@@ -27,7 +50,7 @@ sub test_at_digits {
my @rnd = genrand($digits, $quantity);
my @smp = gensemi($digits, $quantity);
- # verify
+ # verify (can be _really_ slow for 18+ digits)
foreach my $p (@rnd, @smp) {
my @mpxs = prime_factors($p); push @mpxs, $p if $p < 2;
@@ -77,7 +100,7 @@ sub genrand {
my $base = ($digits == 1) ? 0 : int(10 ** ($digits-1));
my $max = int(10 ** $digits);
$max = ~0 if $max > ~0;
- my @nums = map { $base+int(rand($max-$base)) } (1 .. $num);
+ my @nums = map { $base + $rgen->($max-$base) } (1 .. $num);
return @nums;
}
@@ -98,9 +121,9 @@ sub gensemi {
my @factors;
my $n;
while (1) {
- $n = $base + int(rand($max-$base));
- $n += 1 if ($n%2) == 0;
- $n += 3 if ($n%3) == 0;
+ $n = $base + $rgen->($max-$base);
+ # Skip past all multiples of 2, 3, and 5.
+ $n += (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)[$n%30];
@factors = Math::Prime::Util::factor($n);
next if scalar @factors != 2;
next if $factors[0] < $smallest_factor;
diff --git a/examples/test-bpsw.pl b/examples/test-bpsw.pl
new file mode 100755
index 0000000..383822a
--- /dev/null
+++ b/examples/test-bpsw.pl
@@ -0,0 +1,145 @@
+#!/usr/bin/env perl
+use strict;
+use warnings;
+$| = 1; # fast pipes
+
+use Math::Prime::Util;
+use Math::Primality;
+use Config;
+
+my $nlinear = 10000;
+my $nrandom = shift || 20000;
+my $randmax = ~0;
+
+# I was using Math::BigInt::Random::OO, but on my machine:
+# my $gen = Math::BigInt::Random::OO -> new(length => 23);
+# generates only even numbers.
+my $rgen = sub {
+ my $range = shift;
+ return 0 if $range <= 0;
+ my $rbits = 0; { my $t = $range; while ($t) { $rbits++; $t >>= 1; } }
+ while (1) {
+ my $rbitsleft = $rbits;
+ my $U = $range - $range; # 0 or bigint 0
+ while ($rbitsleft > 0) {
+ my $usebits = ($rbitsleft > $Config{randbits}) ? $Config{randbits} : $rbitsleft;
+ $U = ($U << $usebits) + int(rand(1 << $usebits));
+ $rbitsleft -= $usebits;
+ }
+ return $U if $U <= $range;
+ }
+};
+my $rand_ndigit_gen = sub {
+ my $digits = shift;
+ die "Digits must be > 0" unless $digits > 0;
+ my $howmany = shift || 1;
+ my ($base, $max);
+
+ if ( 10**$digits < ~0) {
+ $base = ($digits == 1) ? 0 : int(10 ** ($digits-1));
+ $max = int(10 ** $digits);
+ $max = ~0 if $max > ~0;
+ } else {
+ $base = Math::BigInt->new(10)->bpow($digits-1);
+ $max = Math::BigInt->new(10)->bpow($digits) - 1;
+ }
+ my @nums = map { $base + $rgen->($max-$base) } (1 .. $howmany);
+ return (wantarray) ? @nums : $nums[0];
+};
+
+
+
+if (1) {
+print "OK for first 1";
+my $dig = 1;
+my $i = 9;
+foreach my $n (2 .. $nlinear) {
+ die "MR(2) failure for $n" unless Math::Prime::Util::is_strong_pseudoprime($n,2) == Math::Primality::is_strong_pseudoprime($n,2);
+ die "SLPSP failure for $n" unless Math::Prime::Util::is_strong_lucas_pseudoprime($n) == Math::Primality::is_strong_lucas_pseudoprime($n);
+ die "Prime failure for $n" unless Math::Prime::Util::is_prime($n) == Math::Primality::is_prime($n);
+ if (--$i == 0) {
+ print "0";
+ $dig++;
+ $i = (10 ** $dig) - (10 ** ($dig-1));
+ }
+}
+print " numbers\n";
+print "Testing random numbers from $nlinear to ", $randmax, "\n";
+
+foreach my $r (1 .. $nrandom) {
+ my $n = $nlinear + 1 + int(rand($randmax - $nlinear));
+ my $rand_base = 2 + $rgen->($n-4);
+ die "MR(2) failure for $n" unless Math::Prime::Util::is_strong_pseudoprime($n,2) == Math::Primality::is_strong_pseudoprime($n,2);
+ die "MR($rand_base) failure for $n" unless Math::Prime::Util::is_strong_pseudoprime($n,$rand_base) == Math::Primality::is_strong_pseudoprime($n,$rand_base);
+ die "SLPSP failure for $n" unless Math::Prime::Util::is_strong_lucas_pseudoprime($n) == Math::Primality::is_strong_lucas_pseudoprime($n);
+ die "Prime failure for $n" unless (Math::Prime::Util::is_prime($n)?1:0) == Math::Primality::is_prime($n);
+ print "." if ($r % 256) == 0;
+}
+print "\n";
+
+use bigint try => 'GMP';
+my $big_base = 2**64 + 1;
+my $range = 2**1024 - 1;
+my $end_base = $big_base + $range;
+print "Testing random numbers from $big_base to $end_base\n";
+
+foreach my $r (1 .. int($nrandom/50)) {
+ my $n = $big_base + $rgen->($range);
+ my $rand_base = 2 + $rgen->($n-4);
+ die "MR(2) failure for $n" unless Math::Prime::Util::is_strong_pseudoprime($n,2) == Math::Primality::is_strong_pseudoprime("$n","2");
+ die "MR($rand_base) failure for $n" unless Math::Prime::Util::is_strong_pseudoprime($n,$rand_base) == Math::Primality::is_strong_pseudoprime($n,$rand_base);
+ die "SLPSP failure for $n" unless Math::Prime::Util::is_strong_lucas_pseudoprime($n) == Math::Primality::is_strong_lucas_pseudoprime("$n");
+ die "Prime failure for $n" unless (Math::Prime::Util::is_prime($n)?1:0) == Math::Primality::is_prime("$n");
+ #print "SUCCESS with $n\n";
+ print "." if ($r % 16) == 0;
+}
+print "\n";
+}
+
+use bigint try => 'GMP';
+my $num_rns = 50;
+my $len_rns = 100;
+my $count = -.1;
+
+my @rns; # make the primality tests at least lift a finger.
+while (@rns < $num_rns) {
+ my $n = $rand_ndigit_gen->($len_rns);
+ next unless $n%2 && $n%3 && $n%5 && $n%7 && $n%11 && $n%13;
+ push @rns, $n;
+}
+
+use Benchmark qw/:all/;
+print "Starting benchmarks, $num_rns $len_rns-digit random numbers...\n";
+
+if (1) {
+ print "\nMiller-Rabin, one base:\n";
+ cmpthese($count, {
+ "MPU:PP" => sub { Math::Prime::Util::PP::miller_rabin($_,2) for @rns; },
+ "MPU:GMP" => sub { Math::Prime::Util::GMP::is_strong_pseudoprime($_,2) for @rns; },
+ "MPU" => sub { Math::Prime::Util::is_strong_pseudoprime($_,2) for @rns; },
+ "MP" => sub { Math::Primality::is_strong_pseudoprime("$_","2") for @rns; },
+ });
+}
+
+if (1) {
+ print "\nStrong Lucas test:\n";
+ cmpthese($count, {
+ "MPU:PP" => sub { Math::Prime::Util::PP::is_strong_lucas_pseudoprime($_) for @rns;},
+ "MPU:GMP" => sub { Math::Prime::Util::GMP::is_strong_lucas_pseudoprime($_) for @rns;},
+ "MPU" => sub { Math::Prime::Util::is_strong_lucas_pseudoprime($_) for @rns;},
+ "MP" => sub { Math::Primality::is_strong_lucas_pseudoprime("$_") for @rns;},
+ });
+}
+
+if (1) {
+ print "\nBPSW test:\n";
+ cmpthese($count, {
+ "MPU:PP" => sub { my $sum = 0;
+ do { $sum += ( Math::Prime::Util::PP::miller_rabin($_, 2) &&
+ Math::Prime::Util::PP::is_strong_lucas_pseudoprime($_) )
+ ? 1 : 0 } for @rns; },
+ "MPU:GMP" => sub { Math::Prime::Util::GMP::is_prob_prime($_) for @rns; },
+ "MPU" => sub { Math::Prime::Util::is_prob_prime($_) for @rns;},
+ "MP" => sub { Math::Primality::is_prime("$_") for @rns;},
+ });
+}
diff --git a/examples/test-euler-pari.pl b/examples/test-euler-pari.pl
new file mode 100755
index 0000000..a2a0893
--- /dev/null
+++ b/examples/test-euler-pari.pl
@@ -0,0 +1,37 @@
+#!/usr/bin/env perl
+use strict;
+use warnings;
+$| = 1; # fast pipes
+
+use Math::Prime::Util qw/-nobigint/;
+use Math::Pari;
+
+Math::Prime::Util::prime_precalc(10_000_000);
+my $nlinear = 1000000;
+my $nrandom = shift || 1000000;
+my $randmax = ~0;
+
+# Looks like MPU is 4x faster than Pari for 1 .. 1M, 5x faster for 1M - 10_000M.
+#
+print "OK for first 1";
+my $dig = 1;
+my $i = 9;
+foreach my $n (2 .. $nlinear) {
+ die "failure for eulerphi($n)" unless Math::Prime::Util::euler_phi($n) == Math::Pari::eulerphi($n);
+ die "failure for moebius($n)" unless Math::Prime::Util::moebius($n) == Math::Pari::moebius($n);
+ if (--$i == 0) {
+ print "0";
+ $dig++;
+ $i = (10 ** $dig) - (10 ** ($dig-1));
+ }
+}
+print " numbers\n";
+print "Testing random numbers from $nlinear to ", $randmax, "\n";
+
+while ($nrandom-- > 0) {
+ my $n = $nlinear + 1 + int(rand($randmax - $nlinear));
+ die "failure for eulerphi($n)" unless Math::Prime::Util::euler_phi($n) == Math::Pari::eulerphi($n);
+ die "failure for moebius($n)" unless Math::Prime::Util::moebius($n) == Math::Pari::moebius($n);
+ print "." if ($nrandom % 256) == 0;
+}
+print "\n";
diff --git a/examples/test-factor-gnufactor.pl b/examples/test-factor-gnufactor.pl
new file mode 100755
index 0000000..dc9c6af
--- /dev/null
+++ b/examples/test-factor-gnufactor.pl
@@ -0,0 +1,147 @@
+#!/usr/bin/env perl
+use strict;
+use warnings;
+use Math::Prime::Util qw/factor/;
+use Math::Pari qw/factorint/;
+use File::Temp qw/tempfile/;
+use Math::BigInt try => 'GMP,Pari';
+use Config;
+use autodie;
+use Text::Diff;
+use Time::HiRes qw(gettimeofday tv_interval);
+my $maxdigits = 50;
+$| = 1; # fast pipes
+srand(87431);
+my $num = 1000;
+
+my $rgen = sub {
+ my $range = shift;
+ return 0 if $range <= 0;
+ my $rbits = 0; { my $t = $range; while ($t) { $rbits++; $t >>= 1; } }
+ while (1) {
+ my $rbitsleft = $rbits;
+ my $U = $range - $range; # 0 or bigint 0
+ while ($rbitsleft > 0) {
+ my $usebits = ($rbitsleft > $Config{randbits}) ? $Config{randbits} : $rbitsleft;
+ $U = ($U << $usebits) + int(rand(1 << $usebits));
+ $rbitsleft -= $usebits;
+ }
+ return $U if $U <= $range;
+ }
+};
+
+{ # Test from 2 to 10000
+ print " 2 - 1000"; test_array( 2 .. 1000);
+ print " 1001 - 5000"; test_array( 1001 .. 5000);
+ print " 5001 - 10000"; test_array( 5001 .. 10000);
+}
+
+foreach my $digits (5 .. $maxdigits) {
+ printf "%5d %2d-digit numbers", $num, $digits;
+ my @narray = gendigits($digits, $num);
+ test_array(@narray);
+ $num = int($num * 0.9) + 1; # reduce as we go
+}
+
+sub test_array {
+ my @narray = @_;
+ my($start, $mpusec, $gnusec, $parisec, $diff);
+
+ print ".";
+ $start = [gettimeofday];
+ my @mpuarray = mpu_factors(@narray);
+ $mpusec = tv_interval($start);
+
+ print ".";
+ $start = [gettimeofday];
+ my @gnuarray = gnu_factors(@narray);
+ $gnusec = tv_interval($start);
+
+ print ".";
+ $start = [gettimeofday];
+ my @pariarray = pari_factors(@narray);
+ $parisec = tv_interval($start);
+
+ print ".";
+ die "MPU got ", scalar @mpuarray, " factors. GNU factor got ",
+ scalar @gnuarray, "\n" unless $#mpuarray == $#gnuarray;
+ die "MPU got ", scalar @mpuarray, " factors. Pari factor got ",
+ scalar @pariarray, "\n" unless $#mpuarray == $#pariarray;
+ foreach my $n (@narray) {
+ my @mpu = @{shift @mpuarray};
+ my @gnu = @{shift @gnuarray};
+ my @pari = @{shift @pariarray};
+ die "mpu array is for the wrong n?" unless $n == shift @mpu;
+ die "gnu array is for the wrong n?" unless $n == shift @gnu;
+ die "pari array is for the wrong n?" unless $n == shift @pari;
+ $diff = diff \@mpu, \@gnu, { STYLE => 'Table' };
+ die "factor($n): MPU/GNU\n$diff\n" if length($diff) > 0;
+ my $diff = diff \@mpu, \@pari, { STYLE => 'Table' };
+ die "factor($n): MPU/Pari\n$diff\n" if length($diff) > 0;
+ }
+ print ".";
+ # We should ignore the small digits, since we're comparing direct
+ # Perl functions with multiple command line invocations. It really
+ # doesn't make sense until we're over 1ms per number.
+ printf "OK MPU:%7.2f ms GNU:%7.2f ms Pari:%7.2f ms\n",
+ (($mpusec*1000) / scalar @narray),
+ (($gnusec*1000) / scalar @narray),
+ (($parisec*1000) / scalar @narray);
+}
+
+sub gendigits {
+ my $digits = shift;
+ die "Digits must be > 0" unless $digits > 0;
+ my $howmany = shift;
+ my ($base, $max);
+
+ if ( 10**$digits < ~0) {
+ $base = ($digits == 1) ? 0 : int(10 ** ($digits-1));
+ $max = int(10 ** $digits);
+ $max = ~0 if $max > ~0;
+ } else {
+ $base = Math::BigInt->new(10)->bpow($digits-1);
+ $max = Math::BigInt->new(10)->bpow($digits) - 1;
+ }
+ my @nums = map { $base + $rgen->($max-$base) } (1 .. $howmany);
+ return @nums;
+}
+
+sub mpu_factors {
+ my @piarray;
+ push @piarray, [$_, factor($_)] for @_;
+ @piarray;
+}
+
+sub gnu_factors {
+ my @ns = @_;
+ my @piarray;
+ my $numpercommand = int( (4000-30)/(length($ns[-1])+1) );
+
+ while (@ns) {
+ my $cs = '/usr/bin/factor';
+ foreach my $n (1 .. $numpercommand) {
+ last unless @ns;
+ $cs .= " " . shift @ns;
+ }
+ my $fout = qx{$cs};
+ my @flines = split(/\n/, $fout);
+ foreach my $fline (@flines) {
+ $fline =~ s/^(\d+): //;
+ push @piarray, [$1, split(/ /, $fline)];
+ }
+ }
+ @piarray;
+}
+
+sub pari_factors {
+ my @piarray;
+ foreach my $n (@_) {
+ my @factors;
+ my ($pn,$pc) = @{factorint($n)};
+ # Map the Math::Pari objects returned into Math::BigInts, because Pari will
+ # throw a hissy fit later when we try to compare them to anything else.
+ push @piarray, [ $n, map { (Math::BigInt->new($pn->[$_])) x $pc->[$_] } (0 .. $#$pn) ];
+ }
+ @piarray;
+}
diff --git a/examples/test-factor-mpxs.pl b/examples/test-factor-mpxs.pl
index da6e4ff..ab5f0be 100755
--- a/examples/test-factor-mpxs.pl
+++ b/examples/test-factor-mpxs.pl
@@ -3,12 +3,31 @@ use strict;
use warnings;
$| = 1; # fast pipes
-use Math::Prime::Util qw/factor/;
+use Math::Prime::Util qw/factor -nobigint/;
use Math::Factor::XS qw/prime_factors/;
+use Config;
my $nlinear = 1000000;
my $nrandom = shift || 1000000;
my $randmax = ~0;
+# MFXS is so slow on 17+ digit numbers, skip them.
+$randmax = int(2**55) if $randmax > 2**55;
+
+my $rgen = sub {
+ my $range = shift;
+ return 0 if $range <= 0;
+ my $rbits = 0; { my $t = $range; while ($t) { $rbits++; $t >>= 1; } }
+ while (1) {
+ my $rbitsleft = $rbits;
+ my $U = 0;
+ while ($rbitsleft > 0) {
+ my $usebits = ($rbitsleft > $Config{randbits}) ? $Config{randbits} : $rbitsleft;
+ $U = ($U << $usebits) + int(rand(1 << $usebits));
+ $rbitsleft -= $usebits;
+ }
+ return $U if $U <= $range;
+ }
+};
print "OK for first 1";
my $dig = 1;
@@ -28,7 +47,7 @@ print " numbers\n";
print "Testing random numbers from $nlinear to ", $randmax, "\n";
while ($nrandom-- > 0) {
- my $n = $nlinear + 1 + int(rand($randmax - $nlinear));
+ my $n = $nlinear + 1 + $rgen->($randmax - $nlinear);
my @mfxs = sort { $a<=>$b } prime_factors($n);
my @mpu = sort { $a<=>$b } factor($n);
die "failure for $n" unless scalar @mfxs == scalar @mpu;
diff --git a/examples/test-factor-yafu.pl b/examples/test-factor-yafu.pl
index bec01f2..c17722c 100755
--- a/examples/test-factor-yafu.pl
+++ b/examples/test-factor-yafu.pl
@@ -3,14 +3,31 @@ use strict;
use warnings;
use Math::Prime::Util qw/factor/;
use File::Temp qw/tempfile/;
+use Math::BigInt try => 'GMP,Pari';
+use Config;
use autodie;
use Text::Diff;
-my $maxdigits = (~0 <= 4294967295) ? 10 : 20;
+my $maxdigits = 50;
$| = 1; # fast pipes
my $num = 10000;
my $yafu_fname = "yafu_batchfile_$$.txt";
$SIG{'INT'} = \&gotsig;
+my $rgen = sub {
+ my $range = shift;
+ return 0 if $range <= 0;
+ my $rbits = 0; { my $t = $range; while ($t) { $rbits++; $t >>= 1; } }
+ while (1) {
+ my $rbitsleft = $rbits;
+ my $U = $range - $range; # 0 or bigint 0
+ while ($rbitsleft > 0) {
+ my $usebits = ($rbitsleft > $Config{randbits}) ? $Config{randbits} : $rbitsleft;
+ $U = ($U << $usebits) + int(rand(1 << $usebits));
+ $rbitsleft -= $usebits;
+ }
+ return $U if $U <= $range;
+ }
+};
{ # Test from 2 to 10000
print " 2 - 1000"; test_array( 2 .. 1000);
@@ -22,7 +39,7 @@ foreach my $digits (5 .. $maxdigits) {
printf "%5d %2d-digit numbers", $num, $digits;
my @narray = gendigits($digits, $num);
test_array(@narray);
- $num = int($num * 0.9); # reduce as we go
+ $num = int($num * 0.9) + 1; # reduce as we go
}
sub test_array {
@@ -51,17 +68,23 @@ sub gendigits {
my $digits = shift;
die "Digits must be > 0" unless $digits > 0;
my $howmany = shift;
+ my ($base, $max);
- my $base = ($digits == 1) ? 0 : int(10 ** ($digits-1));
- my $max = int(10 ** $digits);
- $max = ~0 if $max > ~0;
- my @nums = map { $base+int(rand($max-$base)) } (1 .. $howmany);
+ if ( 10**$digits < ~0) {
+ $base = ($digits == 1) ? 0 : int(10 ** ($digits-1));
+ $max = int(10 ** $digits);
+ $max = ~0 if $max > ~0;
+ } else {
+ $base = Math::BigInt->new(10)->bpow($digits-1);
+ $max = Math::BigInt->new(10)->bpow($digits) - 1;
+ }
+ my @nums = map { $base + $rgen->($max-$base) } (1 .. $howmany);
return @nums;
}
sub mpu_factors {
my @piarray;
- push @piarray, [$_, sort {$a<=>$b} factor($_)] for @_;
+ push @piarray, [$_, factor($_)] for @_;
@piarray;
}
@@ -86,7 +109,7 @@ sub yafu_factors {
push @curfactors, $2;
} elsif (/^C\d+ = (\d+)/) {
# Yafu didn't factor this one completely. Sneakily do it ourselves.
- push @curfactors, factor($1);
+ push @curfactors, factor( Math::BigInt->new("$1") );
} elsif (/ans = (\d+)/) {
push @piarray, [shift @ns, sort {$a<=>$b} @curfactors];
@curfactors = ();
--
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