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

_______________________________________________
Pkg-perl-cvs-commits mailing list
Pkg-perl-cvs-commits@lists.alioth.debian.org
http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-perl-cvs-commits

Reply via email to