Trial division (C)

From LiteratePrograms
Jump to: navigation, search
Other implementations: C | Java

Trial division, perhaps the simplest algorithm for factoring and primality testing, simply attempts to divide a number by all possible factors. Divisions that have no remainder (go in evenly) reveal factors. We only have to test up to the square root of n, since if x is a factor, so is n/x. If no factors between 2 and \sqrt{n} (inclusive) are found, the number is prime.


Contents

[edit] Implementations

[edit] With floating-point square root

The simplest possible implementation tries each integer between 2 and the square root of n, computed using the library function sqrt(). We use the modulus operator % to find the remainder after division by each trial factor:

<<trial division with sqrt>>=
unsigned int trial_division_sqrt(unsigned int n) {
    unsigned int sqrt_n = (unsigned int)sqrt((double)n);
    unsigned int x;
    for(x=2; x <= sqrt_n; x++) {
        if ((n % x) == 0) {
            return x;
        }
    }
    return IS_PRIME;
}

<<header files>>=
#include <math.h>  /* For sqrt() */

The reserved value IS_PRIME is not a possible factor:

<<constants>>=
#define IS_PRIME 0

This works but has a few issues. The square root operation, although only done during initialization, is very expensive; more subtlely, we have no way of guaranteeing that we can convert n to a double without losing precision.

[edit] With squaring

Rather than testing if x \leq \sqrt{n}, we can square both sides and test if x^2 \leq n. We can additionally avoid performing a multiply in every loop by noting that (x + 1)2 = x2 + (2(x + 1) − 1). The resulting trial division is ideal for very small numbers:

<<trial division with squaring>>=
unsigned int trial_division_squaring(unsigned int n) {
    unsigned int x, x_squared;
    for(x=2, x_squared=4;
        x_squared > 2*x - 1 && x_squared <= n;
        x++, x_squared += 2*x - 1)
    {
        if ((n % x) == 0) {
            return x;
        }
    }
    return IS_PRIME;
}

There is some subtlety here; we have to ensure that x2 does not overflow before it exceeds n, which might occur if n is close to its maximum value. The approach used above is to make sure that x_squared is larger than the value just added to it; since x never exceeds sqrt(UINT_MAX) this is safe.

[edit] With less trial divisors

The approaches so far work, but are much slower than they need to be: if 2 does not divide n, no multiple of 2 will either, but we try all even numbers in the range. We can save time by skipping all of them but 2:

<<trial division with odd divisors only>>=
unsigned int trial_division_odd(unsigned int n) {
    unsigned int sqrt_n = (unsigned int)sqrt((double)n);
    unsigned int x;
    if ((n % 2) == 0) return 2;
    for(x=3; x <= sqrt_n; x+=2) {
        if ((n % x) == 0) {
            return x;
        }
    }
    return IS_PRIME;
}

More generally, we only have to try prime divisors. If UINT_MAX = 232−1, then the list of all 6542 possible prime factors, each fitting in 16 bits, occupies only 13084 bytes. We could either include a hard-coded list of the primes, or use a simple one-time Sieve of Eratosthenes to generate it, like this:

<<Sieve of Eratosthenes for all 16 bit primes>>=
typedef unsigned short u16;
typedef unsigned int u32;

u16* primes;

void generate_prime_list(int max) {
    unsigned char* is_not_prime = (unsigned char*)calloc(max+1, 1);
    int i, j, num_primes = 1;
    for(i = 2; i <= max; i++) {
        if (is_not_prime[i]) continue;
        num_primes++;
        for(j = i + i; j <= max; j += i) {
            is_not_prime[j] = 1;
        }
    }
    primes = (u16*)malloc(sizeof(u16) * (num_primes + 1));
    j = 0;
    for(i = 2; i <= max; i++) {
        if (!is_not_prime[i]) {
            primes[j] = i;
            j++;
        }
    }
    primes[j] = 0;
    free(is_not_prime);
}

We use a zero terminator to indicate the end of the table. We could make the sieve more efficient in time and space at the expense of clarity by omitting all even numbers from the sieve array. We cast the result of malloc for C++ compatibility. Memory allocation functions require a header:

<<header files>>=
#include <stdlib.h>   /* for malloc, calloc, free */

We can now write a version of trial division for 32-bit numbers that only checks prime divisors:

<<trial division with prime divisors only>>=
u32 trial_division_primes(u32 n) {
    u32 sqrt_n = (u32)sqrt((double)n);
    u16* prime = primes;
    while (*prime != 0 && *prime <= sqrt_n) {
        if ((n % *prime) == 0) {
            return *prime;
        }
        prime++;
    }
    return IS_PRIME;
}

[edit] Test main

This test main allows us to exercise all of our implementations and compare them on numbers supplied on the command line:

<<trial_division.c>>=
header files

constants

trial division with sqrt

trial division with squaring

trial division with odd divisors only

Sieve of Eratosthenes for all 16 bit primes

trial division with prime divisors only

int main(int argc, char* argv[]) {
    int i;
    unsigned int n = atoi(argv[2]);

    if (strcmp(argv[1], "trial_division_sqrt") == 0) {
        for(i=0; i<10000 - 1; i++) trial_division_sqrt(n);
        printf("%u\n", trial_division_sqrt(n));
    } else if (strcmp(argv[1], "trial_division_squaring") == 0) {
        for(i=0; i<10000 - 1; i++) trial_division_squaring(n);
        printf("%u\n", trial_division_squaring(n));
    } else if (strcmp(argv[1], "trial_division_odd") == 0) {
        for(i=0; i<10000 - 1; i++) trial_division_odd(n);
        printf("%u\n", trial_division_odd(n));
    } else if (strcmp(argv[1], "trial_division_primes") == 0) {
        generate_prime_list(65536);
        for(i=0; i<10000 - 1; i++) trial_division_primes(n);
        printf("%u\n", trial_division_primes(n));
    } else {
        printf("Invalid algorithm selection.\n");
    }
    return 0;
}

Since trial division on such small numbers is very quick, we place each call in a long loop for better performance measuring. This also makes the cost of producing the primes table in the last case small on average. We need headers for I/O and string comparison:

<<header files>>=
#include <string.h>
#include <stdio.h>

[edit] Performance

We try each of the algorithms on a number of integer values. Note that trial division's time depends on the magnitude of the smallest factor (the result), as well as the number being factored. Runs were performed on a Gentoo Linux desktop with a Pentium 4 2.8 GhZ processor and 1 GB of RAM.

Input Result sqrt squaring odd primes
2 PRIME 20 ns 3 ns 20 ns 22 ns
23 PRIME 57 ns 41 ns 32 ns 51 ns
281 PRIME 220 ns 208 ns 107 ns 99 ns
2741 PRIME 648 ns 661 ns 354 ns 233 ns
27737 PRIME 2.0 μs 2.1 μs 1.1 μs 0.52 μs
241333 PRIME 5.9 μs 6.2 μs 3.1 μs 1.2 μs
2327911 PRIME 18.2 μs 19.1 μs 9.6 μs 3.1 μs
24988597 PRIME 59 μs 63 μs 31 μs 8.5 μs
290045981 PRIME 200 μs 215 μs 106 μs 25 μs
2144892011 PRIME 550 μs 583 μs 289 μs 61 μs
2144892013 11 144 ns 150 ns 83 ns 87 ns
2144892019 29287 350 μs 368 μs 183 μs 41 μs

As expected, the function that uses only prime divisors quickly dominates for larger inputs, despite the added one-time initialization overhead. Some of the other methods such as the squaring and odds excelled for very small inputs/outputs up to about 500, but squaring is noticably worse for large prime inputs. As expected, both input and output size influence speed, and in particular trial division is extremely quick for numbers with small factors regardless of the input size.

Download code
hijacker
hijacker
hijacker
hijacker