[libmath-prime-util-perl] 02/06: Factoring updates
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:44:20 UTC 2015
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 <dana at 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
More information about the Pkg-perl-cvs-commits
mailing list