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

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

commit f9b919393f4a49dce665a833b0f49f4faf754cd1
Author: Dana Jacobsen <d...@acm.org>
Date:   Thu Jun 7 16:53:16 2012 -0600

    Factoring updates
---
 XS.xs                      |  37 +++++++++----
 examples/bench-factor.pl   | 112 +++++++++++++++++++++++++++++++++++++++
 examples/bench-is-prime.pl |   9 ++--
 factor.c                   | 129 +++++++++++++++++++++++++--------------------
 factor.h                   |  35 +++---------
 util.h                     |   2 +-
 6 files changed, 223 insertions(+), 101 deletions(-)

diff --git a/XS.xs b/XS.xs
index f15eb6e..9c02513 100644
--- a/XS.xs
+++ b/XS.xs
@@ -230,20 +230,35 @@ factor(IN UV n)
     if (n < 4) {
       XPUSHs(sv_2mortal(newSVuv( n ))); /* If n is 0-3, we're done. */
     } else {
+      UV tlim = 19;  /* Below this we've checked */
       UV factor_stack[MPU_MAX_FACTORS+1];
       int nstack = 0;
-      /* Quick trial divisions.  We could do tricky gcd magic here. */
+      /* Quick trial divisions.  Crude use of GCD to hopefully go faster. */
       while ( (n% 2) == 0 ) {  n /=  2;  XPUSHs(sv_2mortal(newSVuv(  2 ))); }
-      while ( (n% 3) == 0 ) {  n /=  3;  XPUSHs(sv_2mortal(newSVuv(  3 ))); }
-      while ( (n% 5) == 0 ) {  n /=  5;  XPUSHs(sv_2mortal(newSVuv(  5 ))); }
-      while ( (n% 7) == 0 ) {  n /=  7;  XPUSHs(sv_2mortal(newSVuv(  7 ))); }
-      while ( (n%11) == 0 ) {  n /= 11;  XPUSHs(sv_2mortal(newSVuv( 11 ))); }
-      while ( (n%13) == 0 ) {  n /= 13;  XPUSHs(sv_2mortal(newSVuv( 13 ))); }
-      while ( (n%17) == 0 ) {  n /= 17;  XPUSHs(sv_2mortal(newSVuv( 17 ))); }
+      if ( (n >= UVCONST(3*3)) && (gcd_ui(n, UVCONST(3234846615) != 1)) ) {
+        while ( (n% 3) == 0 ) {  n /=  3;  XPUSHs(sv_2mortal(newSVuv(  3 ))); }
+        while ( (n% 5) == 0 ) {  n /=  5;  XPUSHs(sv_2mortal(newSVuv(  5 ))); }
+        while ( (n% 7) == 0 ) {  n /=  7;  XPUSHs(sv_2mortal(newSVuv(  7 ))); }
+        while ( (n%11) == 0 ) {  n /= 11;  XPUSHs(sv_2mortal(newSVuv( 11 ))); }
+        while ( (n%13) == 0 ) {  n /= 13;  XPUSHs(sv_2mortal(newSVuv( 13 ))); }
+        while ( (n%17) == 0 ) {  n /= 17;  XPUSHs(sv_2mortal(newSVuv( 17 ))); }
+        while ( (n%19) == 0 ) {  n /= 19;  XPUSHs(sv_2mortal(newSVuv( 19 ))); }
+        while ( (n%23) == 0 ) {  n /= 23;  XPUSHs(sv_2mortal(newSVuv( 23 ))); }
+        while ( (n%29) == 0 ) {  n /= 29;  XPUSHs(sv_2mortal(newSVuv( 29 ))); }
+        tlim = 31;
+      }
+      if ( (n >= UVCONST(31*31)) && (gcd_ui(n, UVCONST(95041567) != 1)) ) {
+        while ( (n%31) == 0 ) {  n /= 31;  XPUSHs(sv_2mortal(newSVuv( 31 ))); }
+        while ( (n%37) == 0 ) {  n /= 37;  XPUSHs(sv_2mortal(newSVuv( 37 ))); }
+        while ( (n%41) == 0 ) {  n /= 41;  XPUSHs(sv_2mortal(newSVuv( 41 ))); }
+        while ( (n%43) == 0 ) {  n /= 43;  XPUSHs(sv_2mortal(newSVuv( 43 ))); }
+        while ( (n%47) == 0 ) {  n /= 47;  XPUSHs(sv_2mortal(newSVuv( 47 ))); }
+        tlim = 53;
+      }
       do { /* loop over each remaining factor */
-        while ( (n >= (19*19)) && (!is_definitely_prime(n)) ) {
+        while ( (n >= (tlim*tlim)) && (!is_definitely_prime(n)) ) {
           int split_success = 0;
-          if (n > UVCONST(60000000) ) {  /* tune this */
+          if (n > UVCONST(10000000) ) {  /* tune this */
             /* For sufficiently large n, try more complex methods. */
             /* SQUFOF (succeeds 98-99.9%) */
             split_success = squfof_factor(n, factor_stack+nstack, 256*1024)-1;
@@ -259,8 +274,8 @@ factor(IN UV n)
             n = factor_stack[nstack];
           } else {
             /* trial divisions */
-            UV f = 19;
-            UV m = 19;
+            UV f = tlim;
+            UV m = tlim % 30;
             UV limit = sqrt((double) n);
             while (f <= limit) {
               if ( (n%f) == 0 ) {
diff --git a/examples/bench-factor.pl b/examples/bench-factor.pl
new file mode 100755
index 0000000..f88404a
--- /dev/null
+++ b/examples/bench-factor.pl
@@ -0,0 +1,112 @@
+#!/usr/bin/env perl
+use strict;
+use warnings;
+use Math::Prime::Util qw/factor/;
+use Math::Factor::XS qw/prime_factors/;
+use Benchmark qw/:all/;
+use List::Util qw/min max reduce/;
+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;
+
+srand(29);
+
+for my $d ( 3 .. $maxdigits ) {
+  print "Factor $howmany $d-digit numbers\n";
+  test_at_digits($d, $howmany);
+}
+
+sub test_at_digits {
+  my $digits = shift;
+  die "Digits has to be >= 1" unless $digits >= 1;
+  die "Digits has to be <= $maxdigits" if $digits > $maxdigits;
+  my $quantity = shift;
+
+  my @rnd = genrand($digits, $quantity);
+  my @smp = gensemi($digits, $quantity);
+
+  # verify
+  foreach my $p (@rnd, @smp) {
+    my @mpxs = prime_factors($p);  push @mpxs, $p if $p < 2;
+
+    verify_factor($p, \@mpxs, [factor($p)], "Math::Prime::Util 
$Math::Prime::Util::VERSION");
+}
+
+
+  #my $min_num = min @nums;
+  #my $max_num = max @nums;
+  #my $whatstr = "$digits-digit ", $semiprimes ? "semiprime" : "random";
+  #print "factoring 1000 $digits-digit ",
+  #      $semiprimes ? "semiprimes" : "random numbers",
+  #      " ($min_num - $max_num)\n";
+
+  my $lref = {
+    "MPU  rand"  => sub { factor($_) for @rnd },
+    "MPU  semi"  => sub { factor($_) for @smp },
+    "MFXS rand"  => sub { prime_factors($_) for @rnd },
+    "MFXS semi"  => sub { prime_factors($_) for @smp },
+  };
+  cmpthese($count, $lref);
+}
+
+sub verify_factor {
+  my $n = shift;
+  my $aref_master = shift;
+  my $aref_check = shift;
+  my $name = shift;
+
+  my @master = sort {$a<=>$b} @{$aref_master};
+  my @check  = sort {$a<=>$b} @{$aref_check};
+  die "Factor $n master fail!" unless $n == reduce { $a * $b } @master;
+  die "Factor $n fail: $name" unless $#check == $#master;
+  die "Factor $n fail: $name" unless $n == reduce { $a * $b } @check;
+  for (0 .. $#master) {
+    die "Factor $n fail: $name" unless $master[$_] == $check[$_];
+  }
+  1;
+}
+
+sub genrand {
+  my $digits = shift;
+  my $num = shift;
+
+  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);
+  return @nums;
+}
+
+sub gensemi {
+  my $digits = shift;
+  my $num = shift;
+
+  my @min_factors_by_digit = (2,2,3,5,7,13,23,47,97);
+  my $smallest_factor = $min_factors_by_digit[$digits];
+  $smallest_factor = $min_factors_by_digit[-1] unless defined $smallest_factor;
+
+  my $base = ($digits == 1) ? 0 : int(10 ** ($digits-1));
+  my $max = int(10 ** $digits);
+  $max = (~0-4) if $max > (~0-4);
+  my @semiprimes;
+
+  foreach my $i (1 .. $num) {
+    my @factors;
+    my $n;
+    while (1) {
+      $n = $base + int(rand($max-$base));
+      $n += 1 if ($n%2) == 0;
+      $n += 3 if ($n%3) == 0;
+      @factors = Math::Prime::Util::factor($n);
+      next if scalar @factors != 2;
+      next if $factors[0] < $smallest_factor;
+      next if $factors[1] < $smallest_factor;
+      last if scalar @factors == 2;
+    }
+    die "ummm... $n != $factors[0] * $factors[1]\n" unless $n == $factors[0] * 
$factors[1];
+    push @semiprimes, $n;
+  }
+  return @semiprimes;
+}
diff --git a/examples/bench-is-prime.pl b/examples/bench-is-prime.pl
index cc9694a..69feee5 100755
--- a/examples/bench-is-prime.pl
+++ b/examples/bench-is-prime.pl
@@ -4,14 +4,17 @@ use strict;
 #use Math::Primality;
 use Math::Prime::XS;
 use Math::Prime::Util;
-use Math::Pari;
+#use Math::Pari;
 #use Math::Prime::FastSieve;
 use Benchmark qw/:all/;
 use List::Util qw/min max/;
 my $count = shift || -5;
 
+my $is64bit = (~0 > 4294967295);
+my $maxdigits = ($is64bit) ? 20 : 10;  # Noting the range is limited for max.
+
 srand(29);
-test_at_digits($_) for (5..15);
+test_at_digits($_) for (5 .. $maxdigits);
 
 
 sub test_at_digits {
@@ -36,7 +39,7 @@ sub test_at_digits {
     #'M::P::FS' => sub { $sieve->isprime($_) for @nums },
     'M::P::U' => sub { Math::Prime::Util::is_prime($_) for @nums },
     'MPU prob' => sub { Math::Prime::Util::is_prob_prime($_) for @nums },
-    'Math::Pari' => sub { Math::Pari::isprime($_) for @nums },
+    #'Math::Pari' => sub { Math::Pari::isprime($_) for @nums },
   });
   print "\n";
 }
diff --git a/factor.c b/factor.c
index 3ef5039..59e4314 100644
--- a/factor.c
+++ b/factor.c
@@ -3,6 +3,8 @@
 #include <string.h>
 #include <assert.h>
 #include <limits.h>
+#define __STDC_LIMIT_MACROS
+#include <stdint.h>
 #include <math.h>
 
 #include "factor.h"
@@ -77,73 +79,84 @@ int trial_factor(UV n, UV *factors, UV maxtrial)
   return nfactors;
 }
 
-static UV gcd_ui(UV x, UV y) {
-  UV t;
 
-  if (y < x) { t = x; x = y; y = t; }
+#if UINT64_MAX > UV_MAX
 
-  while (y > 0) {
-    x = x % y;
-    t = x; x = y; y = t;
-  }
-  return x;
-}
-
-/* if n is smaller than this, you can multiply without overflow */
-#define HALF_WORD (UVCONST(1) << (BITS_PER_WORD/2))
-
-static UV mulmod(UV a, UV b, UV m) {
-  UV r = 0;
-  while (b > 0) {
-    if (b & 1) {
-      if (r == 0) {
-        r = a;
-      } else {
-        r = m - r;
-        r = (a >= r)  ?  a-r  :  m-r+a;
-      }
-    }
-    a = (a > (m-a))  ?  (a-m)+a  :  a+a;
-    b >>= 1;
-  }
-  return r;
-}
+  /* We have 64-bit available, but UV is 32-bit.  Do the math in 64-bit.
+   * Even if it is emulated, it should be as fast or faster than us doing it.
+   */
+  #define addmod(n,a,m)  (UV)(((uint64_t)(n)+(uint64_t)(a)) % ((uint64_t)(m)))
+  #define mulmod(a,b,m)  (UV)(((uint64_t)(a)*(uint64_t)(b)) % ((uint64_t)(m)))
+  #define sqrmod(n,m)    (UV)(((uint64_t)(n)*(uint64_t)(n)) % ((uint64_t)(m)))
 
-/* n^power mod m */
-static UV powmod(UV n, UV power, UV m) {
-  UV t = 1;
-  if (m < HALF_WORD) {
-    n %= m;
-    while (power) {
-      if (power & 1)
-        t = (t*n)%m;
-      n = (n*n)%m;
-      power >>= 1;
-    }
-  } else {
+  static UV powmod(UV n, UV power, UV m) {
+    UV t = 1;
     while (power) {
       if (power & 1)
         t = mulmod(t, n, m);
-      n = (n < HALF_WORD) ? (n*n)%m : mulmod(n, n, m);
+      n = sqrmod(n, m);
       power >>= 1;
     }
+    return t;
   }
-  return t;
-}
 
-/* n + a mod m */
-static UV addmod(UV n, UV a, UV m) {
-  return ((m-n) > a)  ?  n+a  :  n+a-m;
-}
+#else
+
+  /* UV is the largest integral type available (that we know of). */
+
+  /* if n is smaller than this, you can multiply without overflow */
+  #define HALF_WORD (UVCONST(1) << (BITS_PER_WORD/2))
+
+  static UV _mulmod(UV a, UV b, UV m) {
+    UV r = 0;
+    while (b > 0) {
+      if (b & 1) {
+        if (r == 0) {
+          r = a;
+        } else {
+          r = m - r;
+          r = (a >= r)  ?  a-r  :  m-r+a;
+        }
+      }
+      a = (a > (m-a))  ?  (a-m)+a  :  a+a;
+      b >>= 1;
+    }
+    return r;
+  }
 
-/* n^2 mod m */
-#define powsqr(n, m)  (n < HALF_WORD) ? (n*n)%m : mulmod(n,n,m)
+  #define addmod(n,a,m) ((((m)-(n)) > (a))  ?  ((n)+(a))  :  ((n)+(a)-(m)))
+  #define mulmod(a,b,m) (((a)|(b)) < HALF_WORD) ? ((a)*(b))%(m) : 
_mulmod(a,b,m)
+  #define sqrmod(n,m)   ((n) < HALF_WORD)       ? ((n)*(n))%(m) : 
_mulmod(n,n,m)
+
+  /* n^power mod m */
+  static UV powmod(UV n, UV power, UV m) {
+    UV t = 1;
+    if (m < HALF_WORD) {
+      n %= m;
+      while (power) {
+        if (power & 1)
+          t = (t*n)%m;
+        n = (n*n)%m;
+        power >>= 1;
+      }
+    } else {
+      while (power) {
+        if (power & 1)
+          t = mulmod(t, n, m);
+        n = sqrmod(n,m);
+        power >>= 1;
+      }
+    }
+    return t;
+  }
+
+#endif
 
 /* n^power + a mod m */
-#define powmodadd(n, p, a, m)  addmod(powmod(n,p,m),a,m)
+#define powaddmod(n, p, a, m)  addmod(powmod(n,p,m),a,m)
 
 /* n^2 + a mod m */
-#define powsqradd(n, a, m)     addmod(powsqr(n,m), a, m)
+#define sqraddmod(n, a, m)     addmod(sqrmod(n,m),  a,m)
 
 
 /* Miller-Rabin probabilistic primality test
@@ -175,7 +188,7 @@ int miller_rabin(UV n, const UV *bases, UV nbases)
     if ( (x == 1) || (x == (n-1)) )  continue;
 
     for (r = 0; r < s; r++) {
-      x = powsqr(x, n);
+      x = sqrmod(x, n);
       if (x == 1) {
         return 0;
       } else if (x == (n-1)) {
@@ -323,7 +336,7 @@ int holf_factor(UV n, UV *factors, UV rounds)
   for (i = 1; i <= rounds; i++) {
     s = sqrt(n*i);                      /* TODO: overflow here */
     if ( (s*s) != (n*i) )  s++;
-    m = powsqr(s, n);
+    m = sqrmod(s, n);
     /* Cheaper would be:
      *     if (m is probably not a perfect sequare)  continue;
      *     f = sqrt(m);
@@ -370,7 +383,7 @@ int pbrent_factor(UV n, UV *factors, UV rounds)
   }
 
   for (i = 1; i <= rounds; i++) {
-    Xi = powsqradd(Xi, a, n);
+    Xi = sqraddmod(Xi, a, n);
     f = gcd_ui( (Xi>Xm) ? Xi-Xm : Xm-Xi, n);
     if ( (f != 1) && (f != n) ) {
       factors[0] = f;
@@ -408,9 +421,9 @@ int prho_factor(UV n, UV *factors, UV rounds)
   }
 
   for (i = 1; i <= rounds; i++) {
-    U = powsqradd(U, a, n);
-    V = powsqradd(V, a, n);
-    V = powsqradd(V, a, n);
+    U = sqraddmod(U, a, n);
+    V = sqraddmod(V, a, n);
+    V = sqraddmod(V, a, n);
     f = gcd_ui( (U > V) ? U-V : V-U, n);
     if (f == n) {
       if (in_loop++)     /* Mark that we've been here */
diff --git a/factor.h b/factor.h
index 57cedb0..ad0a757 100644
--- a/factor.h
+++ b/factor.h
@@ -22,34 +22,13 @@ extern int pminus1_factor(UV n, UV *factors, UV maxrounds);
 extern int miller_rabin(UV n, const UV *bases, UV nbases);
 extern int is_prob_prime(UV n);
 
-#if 0
-/* Try to make a quick probable prime test */
-static int is_perhaps_prime(UV n)
-{
-  static const UV bases2[2] = {31, 73};
-  static const UV bases3[3] = {2,7,61};
-  if ( (n == 2) || (n == 3) || (n == 5) || (n == 7) )
-    return 2;
-  if ( (n<2) || ((n% 2)==0) || ((n% 3)==0) || ((n% 5)==0) || ((n% 7)==0) )
-    return 0;
-  if (n < 121)
-    return 2;
-  if (n < UVCONST(9080191))
-    return 2*miller_rabin(n, bases2, 2);
-  else
-    return miller_rabin(n, bases3, 3) * ((n < UVCONST(4759123141)) ? 2 : 1);
+static UV gcd_ui(UV x, UV y) {
+  UV t;
+  if (y < x) { t = x; x = y; y = t; }
+  while (y > 0) {
+    t = y;  y = x % y;  x = t;  /* y1 <- x0 % y0 ; x1 <- y0 */
+  }
+  return x;
 }
 
-static int is_perhaps_prime7(UV n)
-{
-  static const UV bases2[2] = {31, 73};
-  static const UV bases3[3] = {2,7,61};
-  /* n must be >= 73 */
-  if (n < UVCONST(9080191))
-    return 2*miller_rabin(n, bases2, 2);
-  else
-    return miller_rabin(n, bases3, 3) * ((n < UVCONST(4759123141)) ? 2 : 1);
-}
-#endif
-
 #endif
diff --git a/util.h b/util.h
index b74d8af..0ecd27c 100644
--- a/util.h
+++ b/util.h
@@ -25,7 +25,7 @@ extern void           free_prime_segment(void);
 
 /* Above this value, is_prime will do deterministic Miller-Rabin */
 /* With 64-bit math, we can do much faster mulmods from 2^16-2^32 */
-#if BITS_PER_WORD == 32
+#if (BITS_PER_WORD == 32) && (UINT64_MAX <= UV_MAX)
 #define MPU_PROB_PRIME_BEST  UVCONST(100000000)
 #else
 #define MPU_PROB_PRIME_BEST  UVCONST(100000)

-- 
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