FORTH GO

Some of the Project Euler problems involve factoring numbers which are either large or small depending on your perspective. Problems are generally limited to 64-bit integers (about 18 digits) which are big numbers for most of us, but in the field of factorization those numbers are terribly small compared to the 100+ digit numbers security protocols deal with. Most advanced methods deal with optimizing the factoring of those huge numbers and don’t mind significant amount of overhead, but I want to know what’s fastest for 64-bit integers.

To find out, I ran some tests on some variations on three basic, low-overhead methods: Trial Division, Fermat’s method, and Pollard’s Rho method. All of these take a long time if the number being factored is actually prime, so it’s worthwhile to add in a fourth component which is a Miller-Rabin primality check. Here are my timing results for 400,000 random 64-bit integers. Actually, only the first test uses 400,000 numbers since I limited each test to 1 hour and extrapolated beyond that.

I was really surprised at how well the Rho method worked in practice. It’s a probabilistic method that’s basically like trial division except it chooses numbers at random instead of sequentially. However, the “random” generator uses a polynomial such that lots of the values can be tested at once using some fancy number theory.

Fermat’s Method works best when there are two divisors near √n, which apparently doesn’t happen very often. Here is my Rho code, which is adapted from some pseudocode in a gamedev forum thread.

   long rhoFactor(long n, int c, int max) {
        int check = 10;
        long x1 = 2;
        long x2 = 4 + c;
        long range = 1;
        long product = 1;
        int terms = 0;

        for (int i = 0; i < max; i++) {
            for (int j = 0; j < range; j++) {
                x2 = (squareMod(x2, n) + c) % n;
                long next = multiplyMod(product, Math.abs(x1 - x2), n);
                if (next == 0) {
                    return Math.max(gcd(product, n), gcd(Math.abs(x1 - x2), n));
                }
                else if (++terms == check || j == range - 1) {
                    long g = gcd(next, n);
                    if (g > 1)
                        return g;
                    product = 1;
                    terms = 0;
                }
                else
                    product = next;
            }

            x1 = x2;
            range *= 2;
            for (int j = 0; j < range; j++) {
                x2 = (squareMod(x2, n) + c) % n;
            }
        }
        return 1;
    }

For the parameters, I used small odd numbers for c, the polynomial constant term, and 16 - 20 for max which limits the generated values at around 2^max. If the factorization fails, I increase c by 2 and try again. For max = 16, it failed to find a factor about once for every 10,000 numbers and never failed twice in my tests. And those numbers had already had any small factors (less than about 50,000) removed with trial division.