[libmath-prime-util-perl] 02/11: Make squfof recurse, and fix some issues with it
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:44:15 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 39f6ea14bb76b0c002f79009d328d862b1dbe838
Author: Dana Jacobsen <dana at acm.org>
Date: Tue Jun 5 21:01:48 2012 -0600
Make squfof recurse, and fix some issues with it
---
XS.xs | 54 ++++++++++++++++++--------------------
examples/bench-factor-semiprime.pl | 2 +-
examples/test-factoring.pl | 37 ++++++++++++++++++++++++++
factor.c | 45 +++++++++++++++++++++++--------
factor.h | 2 +-
lib/Math/Prime/Util.pm | 2 +-
6 files changed, 99 insertions(+), 43 deletions(-)
diff --git a/XS.xs b/XS.xs
index 242522c..2be5131 100644
--- a/XS.xs
+++ b/XS.xs
@@ -226,24 +226,40 @@ erat_simple_primes(IN UV low, IN UV high)
void
factor(IN UV n)
+ PREINIT:
PPCODE:
/* Exit if n is 0 or 1 */
if (n < 2) {
XPUSHs(sv_2mortal(newSVuv( n )));
} else {
- /* Quick trial divisions */
+ /* Quick trial divisions. We could do tricky gcd magic here. */
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 ))); }
- if (n > 1) {
- if (n <= UVCONST(0xFFFFFFFF)) { /* tune this */
- /* trial factorization */
- UV f = 17;
- UV m = 17;
- UV limit = sqrt((double) n);
+ while ( (n%17) == 0 ) { n /= 17; XPUSHs(sv_2mortal(newSVuv( 17 ))); }
+ if (n < (19*19)) {
+ if (n != 1) XPUSHs(sv_2mortal(newSVuv( n )));
+ } else {
+ UV startfactor = 19;
+ UV factor_list[MPU_MAX_FACTORS+1];
+ int n_list_factors = 1;
+ int i;
+ if (n > (UVCONST(0xFFFFFFFF)) ) { /* tune this */
+ /* Use SQUFOF to crack some factors, recurse on large factors */
+ n_list_factors = squfof_factor(n, factor_list, 800000, (373*373)+1);
+ } else {
+ factor_list[0] = n;
+ }
+ for (i = 0; i < n_list_factors; i++) {
+ /* trial factorization for each item in the list */
+ UV f, m, limit;
+ n = factor_list[i];
+ f = startfactor;
+ m = startfactor % 30;
+ limit = sqrt((double) n);
while (f <= limit) {
if ( (n%f) == 0 ) {
do {
@@ -256,27 +272,6 @@ factor(IN UV n)
}
if (n != 1)
XPUSHs(sv_2mortal(newSVuv( n )));
- } else {
- UV factors[MPU_MAX_FACTORS+1];
- UV sqfactors, f1, f2;
- int nfactors;
- int i;
- /* Use SQUFOF to crack a factor off */
- sqfactors = squfof_factor(n, factors+nfactors, 800000);
- assert(sqfactors <= 2);
- if (sqfactors == 1) {
- n = factors[nfactors];
- nfactors += trial_factor(n, factors+nfactors, UV_MAX);
- } else {
- UV n1 = factors[nfactors];
- n = factors[nfactors+1];
- nfactors += trial_factor(n1, factors+nfactors, UV_MAX);
- nfactors += trial_factor(n, factors+nfactors, UV_MAX);
- }
- if (nfactors < 1) croak("No factors");
- for (i = 0; i < nfactors; i++) {
- XPUSHs(sv_2mortal(newSVuv( factors[i] )));
- }
}
}
}
@@ -302,7 +297,8 @@ squfof_factor(IN UV n)
int nfactors;
int i;
PPCODE:
- nfactors = squfof_factor(n, factors, 800000);
+ /* have squfof recurse on each factor found */
+ nfactors = squfof_factor(n, factors, 800000, 2);
if (nfactors < 1)
croak("No factors from squfof_factor");
for (i = 0; i < nfactors; i++) {
diff --git a/examples/bench-factor-semiprime.pl b/examples/bench-factor-semiprime.pl
index 59dcbbd..9942cde 100755
--- a/examples/bench-factor-semiprime.pl
+++ b/examples/bench-factor-semiprime.pl
@@ -6,7 +6,7 @@ $| = 1; # fast pipes
use Math::Prime::Util qw/factor/;
use Math::Factor::XS qw/prime_factors/;
use Benchmark qw/:all/;
-my $digits = shift || 10;
+my $digits = shift || 15;
my $count = shift || -2;
my @min_factors_by_digit = (2,2,3,3,5,11,17,47,97);
diff --git a/examples/test-factoring.pl b/examples/test-factoring.pl
new file mode 100755
index 0000000..f8c5d3d
--- /dev/null
+++ b/examples/test-factoring.pl
@@ -0,0 +1,37 @@
+#!/usr/bin/perl -w
+use strict;
+use warnings;
+$| = 1; # fast pipes
+
+use Math::Prime::Util qw/factor/;
+use Math::Factor::XS qw/prime_factors/;
+
+my $nlinear = 1000000;
+my $nrandom = ~0;
+my $randmax = 1_000_000_000_000;
+
+print "OK for first 1";
+my $dig = 1;
+my $i = 9;
+foreach my $n (2 .. $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;
+ for (0 .. $#mfxs) { die "failure for $n" unless $mfxs[$_] == $mpu[$_]; }
+ if (--$i == 0) {
+ print "0";
+ $dig++;
+ $i = (10 ** $dig) - (10 ** ($dig-1));
+ }
+}
+print " numbers\n";
+print "Testing random numbers from $nlinear to ", $nlinear+$randmax, "\n";
+
+while ($nrandom-- > 0) {
+ my $n = $nlinear + int(rand($randmax));
+ my @mfxs = sort { $a<=>$b } prime_factors($n);
+ my @mpu = sort { $a<=>$b } factor($n);
+ die "failure for $n" unless scalar @mfxs == scalar @mpu;
+ for (0 .. $#mfxs) { die "failure for $n" unless $mfxs[$_] == $mpu[$_]; }
+ print "." if ($nrandom % 1024) == 0;
+}
diff --git a/factor.c b/factor.c
index fa2f4e7..18a269a 100644
--- a/factor.c
+++ b/factor.c
@@ -291,16 +291,17 @@ int pminus1_factor(UV n, UV *factors, UV rounds)
/* My modification of Ben Buhrow's modification of Bob Silverman's SQUFOF code.
* I like Jason P's code a lot -- I should put it in. */
-static IV qqueue[100];
+static IV qqueue[100+1];
static IV qpoint;
static void enqu(IV q, IV *iter) {
qqueue[qpoint] = q;
if (++qpoint >= 100) *iter = -1;
}
-int squfof_factor(UV n, UV *factors, UV rounds)
+int squfof_factor(UV n, UV *factors, UV rounds, UV refactor_above)
{
int nfactors = 0;
+ UV rounds2 = rounds/16;
UV temp;
IV iq,ll,l2,p,pnext,q,qlast,r,s,t,i;
IV jter, iter;
@@ -363,6 +364,7 @@ int squfof_factor(UV n, UV *factors, UV rounds)
r = (int)sqrt((double)q); /* r = floor(sqrt(q)) */
if (q != r*r) continue;
if (qpoint == 0) break;
+ qqueue[qpoint] = 0;
reloop = 0;
for (i=0; i<qpoint-1; i+=2) { /* treat queue as list for simplicity*/
if (r == qqueue[i]) { reloop = 1; break; }
@@ -372,10 +374,15 @@ int squfof_factor(UV n, UV *factors, UV rounds)
break;
} /* end of main loop */
+ if (jter >= rounds) {
+ factors[nfactors++] = n;
+ return nfactors;
+ }
+
qlast = r;
p = p + r*((s - p)/r);
q = (n - (p*p)) / qlast; /* q = (n - p*p)/qlast (div is exact)*/
- for (iter=0; iter<(rounds/16); iter++) { /* unrolled second main loop */
+ for (iter=0; iter<rounds2; iter++) { /* unrolled second main loop */
iq = (s + p)/q;
pnext = iq*q - p;
if (p == pnext) break;
@@ -406,24 +413,40 @@ int squfof_factor(UV n, UV *factors, UV rounds)
p = pnext;
}
- if (iter >= (rounds/20)) {
+ if (iter >= rounds2) {
factors[nfactors++] = n;
return nfactors;
}
if ((q & 1) == 0) q/=2; /* q was factor or 2*factor */
- p = n/q;
- if (p < q) {
+ if ( (q == 1) || (q == n) ) {
+ factors[nfactors++] = n;
+ return nfactors;
+ }
+
+ p = n/q;
+ /* smallest factor first */
+ if (q < p) { t = p; p = q; q = t; }
+
+ /* printf(" squfof found %lu = %lu * %lu in %ld/%ld rounds\n", n, p, q, jter, iter); */
+
+#if 0
+ factors[nfactors++] = p;
+ factors[nfactors++] = q;
+#else
+ if (p >= refactor_above)
+ nfactors += squfof_factor(p, factors+nfactors, rounds, refactor_above);
+ else
factors[nfactors++] = p;
+
+ if (q >= refactor_above)
+ nfactors += squfof_factor(q, factors+nfactors, rounds, refactor_above);
+ else
factors[nfactors++] = q;
- } else {
- factors[nfactors++] = q;
- factors[nfactors++] = p;
- }
+#endif
return nfactors;
}
/* TODO: Add Jason Papadopoulos's racing SQUFOF */
-
diff --git a/factor.h b/factor.h
index 39d2001..dcd4ee3 100644
--- a/factor.h
+++ b/factor.h
@@ -10,7 +10,7 @@ extern int trial_factor(UV n, UV *factors, UV maxtrial);
extern int fermat_factor(UV n, UV *factors);
-extern int squfof_factor(UV n, UV *factors, UV rounds);
+extern int squfof_factor(UV n, UV *factors, UV rounds, UV refactor_above);
extern int pbrent_factor(UV n, UV *factors, UV maxrounds);
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 56ad66a..f83f40f 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -420,7 +420,7 @@ q are the same number of digits), then this will be very fast.
Produces factors, not necessarily prime, of the positive number input. It
is possible the function will be unable to find a factor, in which case a
-single factor (the input) is returned.
+single element, the input, is returned.
=head2 prho_factor
--
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