[libmath-prime-util-perl] 09/11: More factoring changes
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:44:17 UTC 2015
This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.03
in repository libmath-prime-util-perl.
commit a74009093120d159b79f9f8f0939373bd0234a2d
Author: Dana Jacobsen <dana at acm.org>
Date: Wed Jun 6 19:19:30 2012 -0600
More factoring changes
---
MANIFEST | 2 +
examples/bench-factor-extra.t | 98 +++++++++++++++++++++++++++++++++++++++++++
examples/bench-miller-rabin.t | 36 ++++++++++++++++
factor.c | 41 +++++++++++-------
ptypes.h | 44 +++++++++++++++++++
5 files changed, 206 insertions(+), 15 deletions(-)
diff --git a/MANIFEST b/MANIFEST
index ecb1096..f12eb59 100644
--- a/MANIFEST
+++ b/MANIFEST
@@ -15,8 +15,10 @@ sieve.c
util.h
util.c
examples/bench-nthprime.pl
+examples/bench-factor-extra.t
examples/bench-factor-semiprime.pl
examples/bench-is-prime.pl
+examples/bench-miller-rabin.t
examples/test-factoring.pl
t/01-load.t
t/02-can.t
diff --git a/examples/bench-factor-extra.t b/examples/bench-factor-extra.t
new file mode 100755
index 0000000..555ac97
--- /dev/null
+++ b/examples/bench-factor-extra.t
@@ -0,0 +1,98 @@
+#!/usr/bin/perl -w
+
+use strict;
+use Math::Prime::Util;
+use Benchmark qw/:all/;
+use List::Util qw/min max/;
+my $count = shift || -2;
+
+srand(29);
+my $rounds = 100;
+my $sqrounds = 64*1024;
+test_at_digits($_) for (3..10);
+
+
+sub test_at_digits {
+ my $digits = shift;
+
+ die "Digits has to be >= 1" unless $digits >= 1;
+ die "Digits has to be <= 10" if (~0 == 4294967295) && ($digits > 10);
+ die "Digits has to be <= 19" if $digits > 19;
+
+ my @nums = genrand($digits, 1000);
+ #my @nums = gensemi($digits, 1000, 23);
+ my $min_num = min @nums;
+ my $max_num = max @nums;
+
+ # Determine success rates
+ my %nfactored;
+ for (@nums) {
+ $nfactored{'prho'} += int((scalar grep { $_ > 5 } Math::Prime::Util::prho_factor($_, $rounds))/2);
+ $nfactored{'pbrent'} += int((scalar grep { $_ > 5 } Math::Prime::Util::pbrent_factor($_, $rounds))/2);
+ $nfactored{'pminus1'} += int((scalar grep { $_ > 5 } Math::Prime::Util::pminus1_factor($_, $rounds))/2);
+ $nfactored{'fermat'} += int((scalar grep { $_ > 5 } Math::Prime::Util::fermat_factor($_, $rounds))/2);
+ $nfactored{'squfof'} += int((scalar grep { $_ > 5 } Math::Prime::Util::squfof_factor($_, $sqrounds))/2);
+ #$nfactored{'trial'} += int((scalar grep { $_ > 5 } Math::Prime::Util::trial_factor($_))/2);
+ }
+ my $tfac = $nfactored{'fermat'};
+
+ print "factoring 1000 random $digits-digit numbers ($min_num - $max_num)\n";
+ print "Factorizations: ",
+ join(", ", map { sprintf "%s %4.1f%%", $_, 100*$nfactored{$_}/$tfac }
+ grep { $_ ne 'fermat' }
+ sort {$nfactored{$a} <=> $nfactored{$b}} keys %nfactored),
+ "\n";
+
+ my $lref = {
+ "prho" => sub { Math::Prime::Util::prho_factor($_, $rounds) for @nums },
+ "pbrent" => sub { Math::Prime::Util::pbrent_factor($_, $rounds) for @nums },
+ "pminus1" => sub { Math::Prime::Util::pminus1_factor($_, $rounds) for @nums },
+ "fermat" => sub { Math::Prime::Util::fermat_factor($_, $rounds) for @nums },
+ "squfof" => sub { Math::Prime::Util::squfof_factor($_, $sqrounds) for @nums },
+ "trial" => sub { Math::Prime::Util::trial_factor($_) for @nums },
+ };
+ delete $lref->{'fermat'} if $digits >= 9;
+ cmpthese($count, $lref);
+ print "\n";
+}
+
+
+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 $smallest_factor = shift;
+
+ 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-miller-rabin.t b/examples/bench-miller-rabin.t
new file mode 100755
index 0000000..ddb55d9
--- /dev/null
+++ b/examples/bench-miller-rabin.t
@@ -0,0 +1,36 @@
+#!/usr/bin/perl -w
+
+use strict;
+#use Math::Primality;
+use Math::Prime::XS;
+use Math::Prime::Util;
+#use Math::Prime::FastSieve;
+use Benchmark qw/:all/;
+use List::Util qw/min max/;
+my $count = shift || -5;
+
+srand(29);
+test_at_digits($_) for (5..10);
+
+
+sub test_at_digits {
+ my $digits = shift;
+ die "Digits must be > 0" unless $digits > 0;
+
+ 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 $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 "miller_rabin for 1000 random $digits-digit numbers ($min_num - $max_num)\n";
+
+ cmpthese($count,{
+ 'M::P::U' => sub { Math::Prime::Util::miller_rabin($_,2,3,5,7,11,13,17) for @nums },
+ });
+ print "\n";
+}
diff --git a/factor.c b/factor.c
index 2dfc9c0..31b33c3 100644
--- a/factor.c
+++ b/factor.c
@@ -137,9 +137,15 @@ static UV addmod(UV n, UV a, UV m) {
return ((m-n) > a) ? n+a : n+a-m;
}
+/* n^2 mod m */
+#define powsqr(n, m) (n < HALF_WORD) ? (n*n)%m : mulmod(n,n,m)
+
/* n^power + a mod m */
#define powmodadd(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)
+
/* Miller-Rabin probabilistic primality test
* Returns 1 if probably prime relative to the bases, 0 if composite.
@@ -170,7 +176,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 = (x < HALF_WORD) ? (x*x) % n : mulmod(x, x, n);
+ x = powsqr(x, n);
if (x == 1) {
return 0;
} else if (x == (n-1)) {
@@ -299,6 +305,7 @@ int fermat_factor(UV n, UV *factors, UV rounds)
if ( (r != 1) && (r != n) ) {
factors[0] = r;
factors[1] = n/r;
+ assert( factors[0] * factors[1] == n );
return 2;
}
factors[0] = n;
@@ -327,11 +334,12 @@ int pbrent_factor(UV n, UV *factors, UV rounds)
}
for (i = 1; i < rounds; i++) {
- Xi = powmodadd(Xi, 2, a, n);
- f = gcd_ui(Xi - Xm, n);
+ Xi = powsqradd(Xi, a, n);
+ f = gcd_ui( (Xi>Xm) ? Xi-Xm : Xm-Xi, n);
if ( (f != 1) && (f != n) ) {
factors[0] = f;
factors[1] = n/f;
+ assert( factors[0] * factors[1] == n );
return 2;
}
if ( (i & (i-1)) == 0) /* i is a power of 2 */
@@ -363,14 +371,14 @@ int prho_factor(UV n, UV *factors, UV rounds)
}
for (i = 1; i < rounds; i++) {
- U = powmodadd(U, 2, a, n);
- V = powmodadd(V, 2, a, n);
- V = powmodadd(V, 2, a, n);
-
+ U = powsqradd(U, a, n);
+ V = powsqradd(V, a, n);
+ V = powsqradd(V, a, n);
f = gcd_ui( (U > V) ? U-V : V-U, n);
if ( (f != 1) && (f != n) ) {
factors[0] = f;
factors[1] = n/f;
+ assert( factors[0] * factors[1] == n );
return 2;
}
}
@@ -386,18 +394,20 @@ int prho_factor(UV n, UV *factors, UV rounds)
int pminus1_factor(UV n, UV *factors, UV rounds)
{
UV f, i;
- UV b = 13;
+ UV kf = 13;
assert( (n >= 3) && ((n%2) != 0) );
for (i = 1; i < rounds; i++) {
- b = powmod(b, i, n);
- f = gcd_ui(b-1, n);
- if (n == 1) continue;
- if (f == n) break; /* or we could continue */
- factors[0] = f;
- factors[1] = n/f;
- return 2;
+ kf = powmod(kf, i, n);
+ if (kf == 0) kf = n;
+ f = gcd_ui(kf-1, n);
+ if ( (f != 1) && (f != n) ) {
+ factors[0] = f;
+ factors[1] = n/f;
+ assert( factors[0] * factors[1] == n );
+ return 2;
+ }
}
factors[0] = f;
return 1;
@@ -529,6 +539,7 @@ int squfof_factor(UV n, UV *factors, UV rounds)
factors[0] = p;
factors[1] = q;
+ assert( factors[0] * factors[1] == n );
return 2;
}
diff --git a/ptypes.h b/ptypes.h
new file mode 100644
index 0000000..32887e8
--- /dev/null
+++ b/ptypes.h
@@ -0,0 +1,44 @@
+#ifndef MPU_PTYPES_H
+#define MPU_PTYPES_H
+
+#include "EXTERN.h"
+#include "perl.h"
+
+/* From perl.h, wrapped in PERL_CORE */
+#ifndef U32_CONST
+# if INTSIZE >= 4
+# define U32_CONST(x) ((U32TYPE)x##U)
+# else
+# define U32_CONST(x) ((U32TYPE)x##UL)
+# endif
+#endif
+
+/* From perl.h, wrapped in PERL_CORE */
+#ifndef U64_CONST
+# ifdef HAS_QUAD
+# if INTSIZE >= 8
+# define U64_CONST(x) ((U64TYPE)x##U)
+# elif LONGSIZE >= 8
+# define U64_CONST(x) ((U64TYPE)x##UL)
+# elif QUADKIND == QUAD_IS_LONG_LONG
+# define U64_CONST(x) ((U64TYPE)x##ULL)
+# else /* best guess we can make */
+# define U64_CONST(x) ((U64TYPE)x##UL)
+# endif
+# endif
+#endif
+
+
+#ifdef HAS_QUAD
+ #define BITS_PER_WORD 64
+ #define UVCONST(x) U64_CONST(x)
+#else
+ #define BITS_PER_WORD 32
+ #define UVCONST(x) U32_CONST(x)
+#endif
+
+#define MAXBIT (BITS_PER_WORD-1)
+#define NWORDS(bits) ( ((bits)+BITS_PER_WORD-1) / BITS_PER_WORD )
+#define NBYTES(bits) ( ((bits)+8-1) / 8 )
+
+#endif
--
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