Fast Prime Factorization Algorithm

Question

I'm writing a code in C that returns the number of times a positive integer can be expressed as sums of perfect squares of two positive integers.

R(n) is the number of couples (x,y) such that x² + y² = n where x, y, n are all 
non negative integers.

To compute R(n), I need to first find the prime factorization of n.

The problem is that I've tried a lot of algorithm for prime factorization that I can use on C but I need my code to be as fast as possible, so I would appreciate it if anyone can give me what he/she considers as the fastest algorithm to compute the prime factorization of a number as large as 2147483742.

Answer 1

• All prime numbers greater than 2 can be written in the form of (4k ± 1).
• All prime numbers greater than 3 can be written in the form of (6k ± 1).
• All prime numbers greater than 13 can be written in the form of (30k ± { 1 | 7 | 11 | 13 }).

This pattern can be continued/exploited ad-nauseum, but there are diminishing returns. Calculating the square root repeatedly seems inefficient but does significantly reduce the total amount of divisions, when implemented using binary shifts.

public static IEnumerable<T> EnumeratePrimeFactors<T>(this T value) where T : IBinaryInteger<T>, IUnsignedNumber<T> {
    if (BinaryIntegerConstants<T>.Four > value) { yield break; }
    if (BinaryIntegerConstants<T>.Five == value) { yield break; }
    if (BinaryIntegerConstants<T>.Seven == value) { yield break; }
    if (BinaryIntegerConstants<T>.Eleven == value) { yield break; }
    if (BinaryIntegerConstants<T>.Thirteen == value) { yield break; }

    var index = value;

    while (T.Zero == (index & T.One)) { // enumerate factors of 2
        yield return BinaryIntegerConstants<T>.Two;

        index >>= 1;
    }
    while (T.Zero == (index % BinaryIntegerConstants<T>.Three)) { // enumerate factors of 3
        yield return BinaryIntegerConstants<T>.Three;

        index /= BinaryIntegerConstants<T>.Three;
    }
    while (T.Zero == (index % BinaryIntegerConstants<T>.Five)) { // enumerate factors of 5
        yield return BinaryIntegerConstants<T>.Five;

        index /= BinaryIntegerConstants<T>.Five;
    }
    while (T.Zero == (index % BinaryIntegerConstants<T>.Seven)) { // enumerate factors of 7
        yield return BinaryIntegerConstants<T>.Seven;

        index /= BinaryIntegerConstants<T>.Seven;
    }
    while (T.Zero == (index % BinaryIntegerConstants<T>.Eleven)) { // enumerate factors of 11
        yield return BinaryIntegerConstants<T>.Eleven;

        index /= BinaryIntegerConstants<T>.Eleven;
    }
    while (T.Zero == (index % BinaryIntegerConstants<T>.Thirteen)) { // enumerate factors of 13
        yield return BinaryIntegerConstants<T>.Thirteen;

        index /= BinaryIntegerConstants<T>.Thirteen;
    }

    var factor = BinaryIntegerConstants<T>.Seventeen;
    var limit = index.SquareRoot();

    if (factor <= limit) {
        do {
            while (T.Zero == (index % factor)) { // enumerate factors of (30k - 13)
                yield return factor;

                index /= factor;
            }

            factor += BinaryIntegerConstants<T>.Two;

            while (T.Zero == (index % factor)) { // enumerate factors of (30k - 11)
                yield return factor;

                index /= factor;
            }

            factor += BinaryIntegerConstants<T>.Four;

            while (T.Zero == (index % factor)) { // enumerate factors of (30k - 7)
                yield return factor;

                index /= factor;
            }

            factor += BinaryIntegerConstants<T>.Six;

            while (T.Zero == (index % factor)) { // enumerate factors of (30k - 1)
                yield return factor;

                index /= factor;
            }

            factor += BinaryIntegerConstants<T>.Two;

            while (T.Zero == (index % factor)) { // enumerate factors of (30k + 1)
                yield return factor;

                index /= factor;
            }

            factor += BinaryIntegerConstants<T>.Six;

            while (T.Zero == (index % factor)) { // enumerate factors of (30k + 7)
                yield return factor;

                index /= factor;
            }

            factor += BinaryIntegerConstants<T>.Four;

            while (T.Zero == (index % factor)) { // enumerate factors of (30k + 11)
                yield return factor;

                index /= factor;
            }

            factor += BinaryIntegerConstants<T>.Two;

            while (T.Zero == (index % factor)) { // enumerate factors of (30k + 13)
                yield return factor;

                index /= factor;
            }

            factor += BinaryIntegerConstants<T>.Four;
            limit = index.SquareRoot();
        } while (factor <= limit);
    }

    if ((index != T.One) && (index != value)) {
        yield return index;
    }
}

Answer 2

What an odd limit; 2147483742 = 2^31 + 94.

As others have pointed out, for a number this small trial division by primes is most likely fast enough. Only if it isn't, you could try Pollard's rho method:

/* WARNING! UNTESTED CODE! */
long rho(n, c) {
    long t = 2;
    long h = 2;
    long d = 1;

    while (d == 1) {
        t = (t*t + c) % n;
        h = (h*h + c) % n;
        h = (h*h + c) % n;
        d = gcd(t-h, n); }

    if (d == n)
        return rho(n, c+1);
    return d;
}

Called as rho(n,1), this function returns a (possibly-composite) factor of n; put it in a loop and call it repeatedly if you want to find all the factors of n. You'll also need a primality checker; for your limit, a Rabin-Miller test with bases 2, 7 and 61 is proven accurate and reasonably fast. You can read more about programming with prime numbers at my blog.

But in any case, given such a small limit I think you are better off using trial division by primes. Anything else might be asymptotically faster but practically slower.

EDIT: This answer has received several recent upvotes, so I'm adding a simple program that does wheel factorization with a 2,3,5-wheel. Called as wheel(n), this program prints the factors of n in increasing order.

long wheel(long n) {
    long ws[] = {1,2,2,4,2,4,2,4,6,2,6};
    long f = 2; int w = 0;

    while (f * f <= n) {
        if (n % f == 0) {
            printf("%ld\n", f);
            n /= f;
        } else {
            f += ws[w];
            w = (w == 10) ? 3 : (w+1);
        }
    }
    printf("%ld\n", n);

    return 0;
}

I discuss wheel factorization at my blog; the explanation is lengthy, so I won't repeat it here. For integers that fit in a long, it is unlikely that you will be able to significantly better the wheel function given above.

Answer 3

There's a fast way to cut down the number of candidates. This routine tries 2, then 3, then all the odd numbers not divisible by 3.

long mediumFactor(n)
{
    if ((n % 2) == 0) return 2;
    if ((n % 3) == 0) return 3;
    try = 5;
    inc = 2;
    lim = sqrt(n);
    while (try <= lim)
    {
       if ((n % try) == 0) return try;
       try += inc;
       inc = 6 - inc;  // flip from 2 -> 4 -> 2
    }
    return 1;  // n is prime
}

The alternation of inc between 2 and 4 is carefully aligned so that it skips all even numbers and numbers divisible by 3. For this case: 5 (+2) 7 (+4) 11 (+2) 13 (+4) 17

Trials stop at sqrt(n) because at least one factor must be at or below the square root. (If both factors were > sqrt(n) then the product of the factors would be greater than n.)

Number of tries is sqrt(m)/3, where m is the highest possible number in your series. For a limit of 2147483647, that yields a maximum of 15,448 divisions worst case (for a prime near 2147483647) including the 2 and 3 tests.

If the number is composite, total number of divisions is usually much less and will very rarely be more; even taking into account calling the routine repeatedly to get all the factors.