Trial division (C)
From LiteratePrograms
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 (inclusive) are found, the number is prime.
Contents |
Implementations
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.
With squaring
Rather than testing if , we can square both sides and test if . 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.
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; }
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>
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 |