[libmath-prime-util-perl] 08/59: Speedup for powmod/mulmod (helps factor, isprime, etc.)

[Partha P. Mukherjee]

This is an automated email from the git hooks/post-receive script.

ppm-guest pushed a commit to annotated tag v0.10
in repository libmath-prime-util-perl.

commit b4d8c78761a567efc1167e2822fafea5b246593d
Author: Dana Jacobsen <dana at acm.org>
Date:   Thu Jun 28 13:45:30 2012 -0600

    Speedup for powmod/mulmod (helps factor, isprime, etc.)
---
 examples/bench-is-prime.pl |  5 +++--
 examples/bench-pp-sieve.pl |  5 ++++-
 lib/Math/Prime/Util/PP.pm  | 42 +++++++++++++++++++-----------------------
 3 files changed, 26 insertions(+), 26 deletions(-)

diff --git a/examples/bench-is-prime.pl b/examples/bench-is-prime.pl
index a96a214..83fc71d 100755
--- a/examples/bench-is-prime.pl
+++ b/examples/bench-is-prime.pl
@@ -9,6 +9,7 @@ use Math::Prime::Util;
 use Benchmark qw/:all/;
 use List::Util qw/min max/;
 my $count = shift || -5;
+my $numbers = 1000;
 
 my $is64bit = (~0 > 4294967295);
 my $maxdigits = ($is64bit) ? 20 : 10;  # Noting the range is limited for max.
@@ -24,14 +25,14 @@ sub test_at_digits {
   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 .. 1000);
+  my @nums = map { $base+int(rand($max-$base)) } (1 .. $numbers);
   my $min_num = min @nums;
   my $max_num = max @nums;
 
   #my $sieve = Math::Prime::FastSieve::Sieve->new(10 ** $magnitude + 1);
   #Math::Prime::Util::prime_precalc(10 ** $magnitude + 1);
 
-  print "is_prime for 1000 random $digits-digit numbers ($min_num - $max_num)\n";
+  print "is_prime for $numbers random $digits-digit numbers ($min_num - $max_num)\n";
 
   cmpthese($count,{
     #'Math::Primality' => sub { Math::Primality::is_prime($_) for @nums },
diff --git a/examples/bench-pp-sieve.pl b/examples/bench-pp-sieve.pl
index 5464eb3..681c568 100644
--- a/examples/bench-pp-sieve.pl
+++ b/examples/bench-pp-sieve.pl
@@ -17,12 +17,15 @@ my $sum;
 # This is like counting, but we want an array returned.
 # The subs will compute a sum on the results.
 
+# In practice you would probably want to return a ref to your array, or return
+# a ref to your sieve structure and let the caller decode it as needed.
+
 # Times for 100k.
 # Vs. MPU sieve, as we move from 8k to 10M:
 #   Atkin MPTA, Rosetta 3 & 1, Shootout, Scriptol, DO Array, DJ Array, and
 #   InMany all slow down.  Atkin 2 speeds up (from 65x slower to 54x slower).
 # The DJ string methods have almost no relative slowdown, so stretch out their
-# advantage over the others (In Many, DJ Array, DJ Vec, and DO Array).
+# advantage over the other fast ones (In Many, DJ Array, DJ Vec, and DO Array).
 my $pc_subs = {
   "Rosetta 4" => sub {$sum=0; $sum+=$_ for rosetta4($countarg);$sum;},  #    9/s
   "Atkin MPTA"=> sub {$sum=0; $sum+=$_ for atkin($countarg);$sum;},     #   11/s
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index 3c0e81a..e95934c 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -30,6 +30,12 @@ BEGIN {
 }
 sub _maxbits { $_uv_size }
 
+# If $n < $_half_word, then $n*$n will be exact.
+my $_half_word = (~0 == 18446744073709551615) ? 4294967296 :    # 64-bit
+                 (~0 ==           4294967295) ?      65536 :    # 32-bit
+                 (~0 ==                   -1) ?   1000**10 :    # bignum
+                                                         0 ;    # No idea
+
 my $_precalc_size = 0;
 sub prime_precalc {
   my($n) = @_;
@@ -128,10 +134,10 @@ sub is_prime {
 #        .. do something with the prime $n
 #      }
 #
-# We're using method 4, though sadly it is memory intensive.
-#
-# Even so, it is a lot less memory than making an array entry for each number,
-# and the performance is almost 10x faster than a naive array sieve.
+# We're using method 4, though sadly it is memory intensive relative to the
+# other methods.  I will point out that it is 30-60x less memory than sieves
+# using an array, and the performance of this function is over 10x that
+# of naive sieves like found on RosettaCode.
 
 sub _sieve_erat_string {
   my($end) = @_;
@@ -557,6 +563,7 @@ sub nth_prime {
 
 sub _mulmod {
   my($a, $b, $m) = @_;
+  return (($a * $b) % $m) if ($a|$b) < $_half_word;
   my $r = 0;
   while ($b > 0) {
     if ($b & 1) {
@@ -577,7 +584,7 @@ sub _powmod {
   my($n, $power, $m) = @_;
   my $t = 1;
 
-  if ( (~0 == 18446744073709551615) && ($m < 4294967296) ) {
+  if  ($m < $_half_word) {
     $n %= $m;
     while ($power) {
       $t = ($t * $n) % $m if ($power & 1) != 0;
@@ -633,23 +640,12 @@ sub miller_rabin {
     my $x = _powmod($a, $d, $n);
     next if ($x == 1) || ($x == ($n-1));
 
-    if (~0 == 18446744073709551615) {
-      foreach my $r (1 .. $s) {
-        $x = ($x < 4294967296) ? ($x*$x) % $n : _mulmod($x, $x, $n);
-        return 0 if $x == 1;
-        if ($x == ($n-1)) {
-          $a = 0;
-          last;
-        }
-      }
-    } else {
-      foreach my $r (1 .. $s) {
-        $x = ($x < 65536) ? ($x*$x) % $n : _mulmod($x, $x, $n);
-        return 0 if $x == 1;
-        if ($x == ($n-1)) {
-          $a = 0;
-          last;
-        }
+    foreach my $r (1 .. $s) {
+      $x = ($x < $_half_word) ? ($x*$x) % $n : _mulmod($x, $x, $n);
+      return 0 if $x == 1;
+      if ($x == ($n-1)) {
+        $a = 0;
+        last;
       }
     }
     return 0 if $a != 0;
@@ -863,7 +859,7 @@ sub holf_factor {
   for my $i (1 .. $rounds) {
     my $s = int(sqrt($n * $i));
     $s++ if ($s * $s) != ($n * $i);
-    my $m = _mulmod($s, $s, $n);
+    my $m = ($s < $_half_word) ? ($s*$s) % $n : _mulmod($s, $s, $n);
     # Check for perfect square
     my $mcheck = $m & 127;
     next if (($mcheck*0x8bc40d7d) & ($mcheck*0xa1e2f5d1) & 0x14020a);

-- 
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-perl/packages/libmath-prime-util-perl.git