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.

Missie

Just do a prime sieve up to (square root of the upper limit, the upper limit is 2147483742 here), and use the prime table generated for all the prime factorization. You only need to loop at most a few tens of thousands of loop for the worst case - so it's fast enough. – Tactics

(Un)fortunately, there is no "fast factorization" algorithm: you just need to build a table of the first primes up to sqrt(2147483742), and try them all one by one. That's about 6K items. – Kibosh

@dasblinkenlight: There are actually some factorization problem on online judge/competitive programming world that requires Pollard Rho to solve. – Tactics

@Tactics Pollard Rho is probabilistic, there is no guarantee that it finds a factor (of course if you let it run long enough, it would). – Kibosh

@dasblinkenlight I did that already, I got 4971 primes with 46337 being the maximum. But my thinking is that the memory this table will consume might reduce the execution time. Correct me if I'm wrong. – Missie

@F'OlaYinka You are correct, this is one of these not-so-rare cases when you can pay with memory to significantly speed up computation. The table is (relatively) small - in fact, you can use 16-bit numbers instead of full 32-bit ints to save some memory. – Kibosh

@dasblinkenlight You're saying, execution will be faster but more memory will be consumed? – Missie

@dasblinkenlight: I haven't done it first hand, but it has been done by other people. online-judge.uva.es/board/viewtopic.php?f=42&t=42310 – Tactics

@F'OlaYinka: It's speed-memory tradeoff. You usually can use from 32-128MB (depending on problem) for online judge, but you won't be using all of the memory in most cases. You will usually safe from the threshold when you use correct algorithm or overshoot it by miles. – Tactics

@Tactics that's fair enough, i guess i will just stick with this method – Missie

@F'OlaYinka There are 4792 primes not exceeding 46340. If you just typo-swapped the two middle digits, you're one off (forgot to count 2 since the sieve contained only odd numbers?). If not, you have a serious mistake in your sieve. – Vinavinaceous

Answer 1

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.

Vilayet

Answer 2

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.

Wellfounded

Answer 3

• 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;
    }
}

Magdalenemagdalenian