Slow factorization of numbers with many small prime factors

[trizen]

Hi Dana.

I noticed that the prime factorization of large integers that have many small prime factors seems to be quite slow, but I think we can do better.

Example:

use 5.010;
use strict;
use warnings;

use Time::HiRes qw(gettimeofday tv_interval);
use Math::Prime::Util::GMP qw(factor consecutive_integer_lcm);

my $n = consecutive_integer_lcm(38861);

say "Length of n = ", length($n);

my $t0 = [gettimeofday];
my @factors = factor($n);
my $elapsed = tv_interval ( $t0, [gettimeofday]);

say "bigomega(n) = ", scalar(@factors);
say "Factorization took $elapsed seconds.";

Output:

Length of n = 16875
bigomega(n) = 4175
Factorization took 7.991219 seconds.

The same factorization in PARI/GP takes only ~40ms:

? vecsum(factor(lcm([1..38861]))[,2])
%1 = 4175
? ##
  ***   last result: cpu time 41 ms, real time 45 ms.
? 

As a solution, I propose the following adaptive trial-division algorithm, which efficiently finds small prime factors of n, and it also stops early when n no longer has small prime factors:

use 5.020;
use strict;
use warnings;

use Math::GMPz;
use Time::HiRes qw(gettimeofday tv_interval);
use Math::Prime::Util::GMP qw(:all);

use experimental qw(signatures);

sub fast_trial_factor ($n, $L = 1e4, $R = 1e6) {

    $n = Math::GMPz->new("$n");

    my @factors;

    my @P = sieve_primes(2, $L);

    my $g = Math::GMPz::Rmpz_init();
    my $t = Math::GMPz::Rmpz_init();

    while (1) {

        Math::GMPz::Rmpz_set_str($g, vecprod(@P), 10);
        Math::GMPz::Rmpz_gcd($g, $g, $n);

        # Early stop when n seems to no longer have small factors
        if (Math::GMPz::Rmpz_cmp_ui($g, 1) == 0) {
            last;
        }

        # Factorize n over primes in P
        foreach my $p (@P) {
            if (Math::GMPz::Rmpz_divisible_ui_p($g, $p)) {

                Math::GMPz::Rmpz_set_ui($t, $p);
                my $valuation = Math::GMPz::Rmpz_remove($n, $n, $t);
                push @factors, ($p) x $valuation;

                # Stop the loop early when no more primes divide `g` (optional)
                Math::GMPz::Rmpz_divexact_ui($g, $g, $p);
                last if (Math::GMPz::Rmpz_cmp_ui($g, 1) == 0);
            }
        }

        # Early stop when n has been fully factored,
        # or when the trial range has been exhausted
        if ($L > $R or Math::GMPz::Rmpz_cmp_ui($n, 1) == 0) {
            last;
        }

        @P = sieve_primes($L + 1, $L << 1);
        $L <<= 1;
    }

    return (\@factors, $n);
}

my $n = consecutive_integer_lcm(38861);
# $n = vecprod($n, Math::GMPz->new(2)**128 + 1);

say "Length of n = ", length($n);

my $t0 = [gettimeofday];
my ($f, $r) = fast_trial_factor($n);
my $elapsed = tv_interval($t0, [gettimeofday]);

say "remainder = $r";
say "bigomega(n) = ", scalar(@$f);
say "Factorization took $elapsed seconds.";

Output:

Length of n = 16875
remainder = 1
bigomega(n) = 4175
Factorization took 0.047779 seconds.

[danaj]

Good idea. I'm thinking for internal use I'll make a trial factor iterator. We're using pre-computed gcds up to 32k to do quick tests, but it's kind of a mess actually doing the trial division efficiently. Factor trees are used later.

I'll have to look at the different use cases and see what I can do. Nice algorithm, should be more efficient to use once past the small gcds.

[danaj]

Ah, I was thinking of the ECPP code which used gcds. That idea got put into the main factoring code only a few weeks ago (in response to another issue you pointed out).

Current code:

$ time perl -Iblib/lib -Iblib/arch ~/tmp/g27-1.pl 
Length of n = 16875
bigomega(n) = 4175
Factorization took 0.405609 seconds.

Your code:

$ time perl -Iblib/lib -Iblib/arch ~/tmp/g27-2.pl 
Length of n = 16875
remainder = 1
bigomega(n) = 4175
Factorization took 0.037066 seconds.

A big improvement over what it was before, but it ends trial factoring quickly. Simply expanding it to 64k gives:

$ make >/dev/null && time perl -Iblib/lib -Iblib/arch ~/tmp/g27-1.pl 
Length of n = 16875
bigomega(n) = 4175
Factorization took 0.033636 seconds.

But I'm sure I could construct examples that would want better factoring after that. I also think there is some micro-improvements to be made in the trial factoring using the gcd result. I'll have to do more testing.