Optimization of prime factorization in C

I am trying to teach myself some programming in C by doing the problems of project Euler. I am trying to find the largest prime factor of the number 600851475143 (Euler problem 3). I have written a code which successfully does this for numbers which are not too large:

#include <stdio.h>
int primearray(int,int);

main()
{
int alpha, k, j, i=1, prime[200000], count = 0;
prime[0]=2;

/* Define number which counts number of primes in interval [0,alpha] */
int totprime=1;

printf("Enter the number of which you want the prime factorization\n");
scanf("%d",&alpha);

/* This procedure makes a list of all prime numbers up to alpha */
for (j=3 ; j<alpha ; j++)
{
for (k=0 ; k<i ; k++)
{
if ((j/prime[k])*prime[k]==j)
break;
else
count = count + 1;

}
if (count == i)
{
prime[i]=j;
i = i + 1;
count = 0;
totprime = totprime + 1;
}
else
count = 0;
}

/* Print prime numbers */
printf("The primes are:\n");
for (j=0 ; j<totprime ; j++) printf("%d\n",prime[j]);

printf("Total number of primes: %d\n",totprime);

int primefac[totprime];

for (j=0 ; j<totprime ; j++)
{
if ( (alpha/prime[j])*prime[j] == alpha)
{
primefac[count] = prime[j];
count = count + 1;
}
else
continue;
}
printf("Primefactors are\n");
for (j=0 ; j<count ; j++) printf("%d\n",primefac[j]);
}


But it fails to work for large numbers, probably because it has been programmed horribly inefficient. I was wondering if you can help me think of a way to improve this program such that it is able to handle large numbers in a reasonable amount of time. All tips are welcome!

-
For the first j loop you can count up by 2 since you only need to check odd numbers for primality. On the second j loop, count backwards through the primes until you find one that divides alpha since you only need the largest. –  Godric Seer Apr 26 at 16:34
what OS and what compiler are you using? –  Ahmed Masud Apr 26 at 23:29
there are SOOO many issues with your code though :) –  Ahmed Masud Apr 27 at 14:02
@AhmedMasud: Please elaborate. –  Paul Apr 27 at 17:43

Whether your code is efficient or not, it will not work for any numbers over 32767 as written. This is because the int data type is a signed type of 16 bit length. One bit is used for the sign and 15 are used for the value of the integer, making the largest storable integer 215-1=32767. If you wish to support numbers larger than this, you will need to use a different integer type such as long. Assuming you do not need to store negative integers, you should also use an unsigned type such as unsigned long or even unsigned long long which can store integer values up to 232-1=14294967295 or 264-1=18446744073709551615 respectively.

Edit:
As @Ahmed pointed out in the comments, the C standards do not specify the sizes of integer types so portability can be a problem. To make sure you get the behavior you expect, I would suggest including the <inttypes.h> header. <inttypes.h> subsumes <stdint.h> and includes additional macros, especially relating to I/O. As described here, <stdint.h> defines the types:

int8_t    // 8 bit integer
int16_t   // 16 bit integer
int32_t   // 32 bit integer
uint8_t    // unsigned 8 bit integer
uint16_t   // unsigned 16 bit integer
uint32_t   // unsigned 32 bit integer


Most implementations also provide

int64_t   // 64 bit integer
uint64_t   // unsigned 64 bit integer


These types are portable and reliable across architectures and compilers. I think using an external library is overkill in this case, but if needed, you can see @Ahmed's link in the comments.

-
Size of int is highly dependent on local CPUs as well as the C compiler used by OP. It is better to use a big number library for larget numbers such as gmplib.org –  Ahmed Masud Apr 26 at 23:28
Good point, I added an edit addressing this. I think an external library is overkill in this case though. –  Doug Lipinski Apr 27 at 0:06

Since Doug has pointed out the need for a bigger integer datatype to do the problem you want to tackle, let's talk about the logic of your program and improvements.

Your approach is (A) input alpha as the number to factor, (B) make a list of all prime numbers less than alpha, and (C) see which of these divide alpha and add the ones that do to an array primefac[ ].

This factoring technique is a version trial division. That is a good place to start, at least when you know that alpha itself is not a prime, something you probably want to assume from the context of the Euler Problem. However there are relatively simple and fast ways to be sure that alpha has factors (is composite) before we knock ourselves out trying to factor it, so for general purposes you might want to look into the Miller-Rabin primality test. But that can be another Question.

The basic thing to fix in your program (after longer datatypes) is checking for repeated prime factors. Once you determine that a prime divides alpha, you just continue with the next prime. This means you never determine how many times a repeated prime factor occurs in alpha, and that may be fine for your application, but it also means you never take advantage of making alpha smaller by any of the prime factors that divide it. Smaller numbers to check generally mean fewer checks to make, as we explain in the next paragraph.

You have set the goal to check alpha for divisibility of all primes less than alpha, but this is overkill. You really only want to check for divisibility by primes that are less than or equal to the square root of alpha. Any factor of alpha that is left after checking trial division by primes less than or equal to the square root of alpha must itself be a prime. As the size of alpha drops, so to does the size of its square root, and we can stop when the trial divisor gets to be bigger than this falling target (square root of alpha).

Moreover you are using trial division to build the list of primes to check for divisibility into alpha. This is overkill. You do treat 2 as a special case (known prime), but then you start checking all the consecutive integers starting with 3 for primality. As Godric commented above, you can easily increment by 2 instead of 1 when creating prime candidates, since there will never be another even prime.

Indeed the testing of these candidates for primality wastes too much time. You can generate some odd numbers that include all the odd primes quickly, and use these as your trial divisors. Even though these will include some odd composite numbers, you will save time because you skip testing the trial divisors for primality and go directly into testing alpha for divisibility by these trial divisors.

So I recommend these steps:

(1) Find out if 2 divides alpha, and if so set alpha = alpha/2 until alpha is no longer divisible by 2.

(2) Using the possibly revised value of alpha from the last step, find out if it is divisible by 3, and if so set alpha = alpha/3 until alpha is no longer divisible by 3.

(3) Starting with 5, and continuing until the trial divisors are greater than the square root of the current alpha, generate trial divisors by alternately adding 2 and 4. I.e. generate 5,7=5+2,11=7+4,13=11+2,17=13+4, etc. If alpha is divisible by the current trial divisor T, set alpha = alpha/T until alpha is no longer divisible by T.

(4) The stopping criterion, that T is greater than the square root of current alpha, does not require a square root operation. Instead it makes sense to check whether alpha/T is less than T and stop if that is the case. This is fairly efficient because you would compute alpha/T anyway in the course of checking (exact) divisibility by T as your code already does it, asking if alpha == (alpha/T)*T.

Summary: bigger datatypes, no list of primes, check only up to the square root of alpha. If the final value of alpha is greater than 1, then that is the largest prime factor of the originally input alpha. Otherwise the largest prime factor is the last successful trial divisor T.

-
Thanks for a great and clear reply! This really helped me a lot! –  Funzies Apr 29 at 16:16

Since you are trying to teach yourself C, let's actually show you what the issues are in your code and show you how to attack this problem properly utilizing the complete power of C.

#include <stdio.h>


The following declaration is useless

int primearray(int,int);


The declaration below should be ANSI compliant

main()


int main()

{


You should NOT use such a big array as a stack variable and actually either declare it as a global variable or allocate it dynamically ( much better solution )

  int alpha, k, j, i=1, prime[200000], count = 0;
prime[0]=2;

/* Define number which counts number of primes in interval [0,alpha] */
int totprime=1;


you should divide the code up into separate functions so that you can deal with each piece simply and easily

For instance a function that takes a number and finds all primes using whatever method, that way you can start out with a simple function and as you learn better methods you can actually make that routine faster.

  printf("Enter the number of which you want the prime factorization\n");
scanf("%d",&alpha);

/* This procedure makes a list of all prime numbers up to alpha */


The block below has many issues.

1. First it's suboptimal because it tests far too many numbers for primality.

2. The test for primality is incorrect and will not result in the value always equal to j when they're supposed to be. Also it's very CPU intenstive you can use modulo operator to get results more effectively

3. The tests for storing the primes is wierd.

4. This entire block should be in its own function.

  for (j=3 ; j<alpha ; j++)
{

for (k=0 ; k<i ; k++)
{
if ((j/prime[k])*prime[k]==j)
break;
else
count = count + 1;

}

if (count == i)
{
prime[i]=j;
i = i + 1;
count = 0;
totprime = totprime + 1;
}
else
count = 0;
}


Don't comment obvious code the code shows that you're going to print the primes adding a comment about is frivilous and detracts

/* Print prime numbers */
printf("The primes are:\n");

for (j=0 ; j<totprime ; j++) printf("%d\n",prime[j]);


Using for loop the way you did above is poor form and you should not use the printf statement on the same line.

Generally speaking: Your tabing and block style is inconsistent when you are coding you should use a consistent way to mark your blocks code is read far more often than it's written using a consistent style helps follow the logic

  printf("Total number of primes: %d\n",totprime);


Dynamic arrays are not supported in C. The code below is not standard. if your platform/compiler supports it you will get inconsistent results when you port the code.

  int primefac[totprime];


The block below has many problems.

  for (j=0 ; j<totprime ; j++)
{


This test will fail because you are assuming things about how integer division works (and whether or not it results in a fraction) Not to mentions that it's an extremely expensive operation

    if ( (alpha/prime[j])*prime[j] == alpha)
{
primefac[count] = prime[j];
count = count + 1;
}


The else below is completely not required by any logic

    else
continue;
}

printf("Primefactors are\n");


Same issue as above

for (j=0 ; j<count ; j++) printf("%d\n",primefac[j]);


This should line up with the main function }

/* STILL EDITING */ The correct approach coming up.

Okay now let's figure out how approach this problem correctly.

Stating the problem simply

You want to find the largest prime factor of 600851475143.

That is to say you want to find maximum value of N

$$\left. 600851475143 \equiv 0 \pmod{N} \: \middle| \: {N \not \equiv 0} \pmod {p} \; \forall p \in \mathbb{Z}^+ \right.$$

In C the modulo operator (%) (actually it's the remainder operator) provides the necessary capability.

So how do we do this, first let's create a function that does that. We know that if $a \times b \equiv N \; \forall {a, b} \in \mathbb{Z}^+$ then either $a \leq \sqrt{N}$ or $b \leq \sqrt{N}$

So a simple (crude) algorithm would be to see if there is a prime number $p$ such that: $\left. p \leq \sqrt{600851475143} \; \big| \; p \, \vert \, 600851475143 \right.$

A correct approach

I will not provide you with a solution but walk you through some basic concepts.

First things first, it's a good idea to split things up into smaller chunks and treat each as a function. If you define your approach step by step then a good rule of thumb is to create one function per step. If it's a complex step, with simpler substeps, then that function should be divided into multiple sub functions. You may not get it right the first time, that's okay, you will develop better judgement as you go along as to how things should get divided up.

Okay let's look at the problem step by step using our simple algorithm using a flowchart

Now essentially each box in the flowchart above is going to the statements of your main function:

/**
* name: primes_upto
* parameters:
*    M: finds all primes until M and stores them in zero-terminated list
*    list: the pointer to an array of ints
*          if list is NULL then it will allocate enough memory
*          to store the primes + 1
*          if list is not NULL it will truncate insert a 0 at
*          the index where value of prime is > M
*
* return value: returns number of primes on success or 0 on failure
*/
int primes_upto(int k, int **list);

int main() {

const int original_N = 600851475143;
int N, a, S;
int *L = NULL;    /* list of primes we need to fill */
int *ap = NULL; /* a memory pointer to integers that will point to the next prime */

N = original_N;

do {
S = sqrt(N); /* we have to provide this */

primes_upto(S, &L); /* we have to provide this function */

/* we will make sure that the primes_upto function
injects a zero at the end of the list. Since zero is not a prime
we can use it as a marker to indicate the end of the list
and that way we don't have to worry about how long the list is. */

for (ap = L; ap && *ap != 0; ap++) {
int a = *ap;
if (N % a == 0) {
N = N / a;
break;
}
}

} while (ap && *ap != 0);

printf ("The largest number that divides %d is %u\n", original_N, N);

}

-

you should consider making use of the modulus (%) operator instead of (j/prime[k]*prime[k] == j) say (j%prime[k] == 0). I know this is supported in C++, correct me if it is not present in C.

Furthermore you only need to check primes up to the square root of your target number. So say you want to factor a number N, you only need primes up to sqrt(N) since its guaranteed if the number is not composite that one factor will exist on that interval.

Here is an example:

Say you want to find all prime numbers less than 100. We just need to generate primes 2,3,5,7 through a process of your choice. and then if you write out all the numbers from 8 - 100, and just cross out the multiples of 2,3,5 and 7 (without considering ANY additional primes) the only numbers that won't be crossed out will be the remaining prime numbers.

List every prime number less than 775146

and check each of these against 600851475143 as you generate them.

In order to handle large numerical datatypes i recommend either importing externally OR:

take a look at the vector class (well thats c++) and if thats off limits to then you are going to want to create a numerical type function which essentially converts a string ('101101343433') into a linked list of digits and then create functinos for adding, subtracting, dividing, and multiplying over this linked list device.

LASTLY IF YOU REALLY REALLY NEED HIGH SPEED.

When you generate your list of primes multiply every k consecutive primes together (choose k as something such that the average size of the numbers yielded has the same number of digits as the number you want to factor). And then find the GCF of this target number and your groups of k primes multiplied together. The punchline here is that instead of having to do k separate modulus checks (is Num % prime == 0) or in your case k seperate divisions and multiplications you do k smaller multiplications (since they are smaller numbers) thereby shaving off some time. The euclidean algorithm should be used to find GCF.

Which is:

given A, B to find GCF(a,b)

while(true): if( A >= B)
A = A %B;

if A == 0
return B


if( B >= A) B = B%A;

if B == 0 return A

Thus here is an example for the checking the factors of the number 143.

List

2,3,5,7,11 (these are all primes less than sqrt(143) generated using the same eratosthenes protocol that you have in your code (the sieve of eratosthenes is the formal name for the trial division algorithm you have implemented) )

group:

2x3x5 = 30, 5x7x11 = 385

GCF(143,30) = 1 GCF(143,385) = 11

11 = factor

143/11 = 13

13 is another factor.

13 > 11

13 is the largest factor.

-

There's a more advanced version of the Pythagorean sieve that tracks patterns of composite numbers, and uses those patterns to reduce the number of necessary tests.

For a simple example, consider the numbers 0 through 5. Only 1 and 5 are relatively prime to 6=2*3. You can use this information to eliminate some tests for the numbers 6 through 11. 2+6 = 8, and 2 is not relatively prime to 6, so neither can 8 be, so it must be composite. The same goes for 3+6 = 9 and 4+6 = 10. The others, 1+6 = 7, and 5+6 = 11, might be prime. You need to test these, but in this case they turn out to be prime. The point is you don't need to test 8, 9, or 10.

Similarly for the next batch of 6, the numbers 12 through 17, they operate similarly to the previous interval because they're equivalent to 6 through 11 in mod 6 arithmetic. You need to test 1+12 and 5+12, but not 0+12, 2+12, 3+12, 4+12.

You may continue for the intervals [18,23] and [24,29].

At this point you have covered the interval [0,29]. The next number 30 = 2*3*5. You might be noticing a pattern here. We consider ever expanding groups, all nested within the next larger interval, which will be the product of the first N+1 primes.

So you have identified which numbers mod 30 must be composite: {0,2,3,4,6,8,9,10,12,14,15,16,18,20,21,22,24,26,27,28}. 25 is the first possible prime that turns out to be composite, but that's OK because you would have tested it according to the pattern you found. Because 25 is composite you add it to the list of numbers you will eliminate when you consider the next larger interval.

Similarly to the way you eliminated composites in the interval [0,29], you can do the same for [30,59], [60,89], [90,119], and so forth up to [180,209]. You test possible primes along the way and add to the list the ones that turn out to be composite.

Now your next larger interval 210 = 2*3*5*7, the product of the first 4 primes. You repeat the process up to 2310 = 2*3*5*7*11, then 30030, 510510, etc.

To save space you can record flags, indicating which numbers are composite or not, as individual bits within integers.

If you're interested in other simple algorithms for prime testing, you can check Pollard's Rho algorithm. It's a little abstruse but easy to implement. Wikipedia has a nice article describing it.

-