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 <[email protected]> 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 _______________________________________________ Pkg-perl-cvs-commits mailing list [email protected] http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-perl-cvs-commits
