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.

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. — nhahtdh
(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. — Sergey Kalinichenko
@dasblinkenlight: There are actually some factorization problem on online judge/competitive programming world that requires Pollard Rho to solve. — nhahtdh
@nhahtdh Pollard Rho is probabilistic, there is no guarantee that it finds a factor (of course if you let it run long enough, it would). — Sergey Kalinichenko
@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. — Daniel Fischer

user448810 user448810 · Accepted Answer · 2012-10-06T12:19:14

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.

Fast Prime Factorization Algorithm

2 Answers