3
\$\begingroup\$

A book1 that I'm reading states Goldbach's conjecture as follows:

  1. Every even integer > 2 could be represented as a sum of two prime numbers.
  2. Every integer > 17 could be represented as a sum of three unique prime numbers.
  3. Every integer could be represented as a sum consisted of maximally six prime numbers.
  4. Every odd integer > 5 could be represented as a sum of three prime numbers.
  5. Every even integer could be represented as a difference between two prime numbers.

and then asks for a code that checks each of the five statements. Here is the code:

Goldbach.c

// upper bound of the interval we search for primes
#define MAX 100

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "HelperArray.h"
#include "HelperPrime.h"
#include "Goldbach.h"

int main(){
    testGoldbackHypothesis();
    return 0;
}

HelperArray.h

#ifndef HELPERARRAY_H
#define HELPERARRAY_H

/*
Function: pi()

It returns the approximate
number of primes up to the 
paramter x. 
(To be used to estimate the size of the array to store primes.)
*/
long int pi(int x){
    return x / (log((double) x) - 1);
}

//-----------------------------------------------------------

/*
Function: arraySize()

It returns the size of the 
array that will hold primes.

x/logx always > prime number density pi(x)/x.
*/
long int arraySize(int x){
    return x / log((double) x);
}
//-----------------------------------------------------------

/*
initArray();

*/
void initArray(int* primes, unsigned int size, int initValue){
    unsigned int i;
    for (i = 0; i < size; ++i){
        primes[i] = initValue;  
    }
}
//-----------------------------------------------------------

/*
Function: printArray();

*/
void printArray(int* primes, unsigned int size){
    unsigned int i;
    int* ptrToArray = primes;
    int fieldWidth = 1 + log10((double)MAX);

    printf("{");
    for (i = 0; i < size; ++i){
        printf("%*d", fieldWidth, ptrToArray[i]);

        if (i < size - 1){
            // exclude unassigned values at the end
            if (ptrToArray[i+1] == 0){
                break;
            }
            printf(", ");
        }

        if (i % 20 == 0 && i != 0){
           printf("\n");
        }
    }
    printf(" }\n");
}
//-----------------------------------------------------------

/*
Function: binarySearch()

It returns true if targer is 
found in the array named primes. 
Otherwise returns false.
*/
char binarySearch(unsigned int target, int* primes, unsigned int size){
    int* ptrToArray = primes;
    int first = 0;
    int last = size;

    while (first <= last){
        unsigned int middle = first + (last - first) / 2;

        if (ptrToArray[middle] == target){
           return 1;
        }

        if (ptrToArray[middle] < target){
            first = middle + 1;
        }else{
            last =  middle - 1;
        }
    }
    return 0;
}

#endif

HelperPrime.h

#ifndef HELPERPRIME_H
#define HELPERPRIME_H

/*
Function isPrime();

It returns true if the argument
is prime, otherwise returns false.
*/
char isPrime(int n){
unsigned int denom = 2;

    if (n < 2){
        return 0;
    }

    if (n == 2){
        return 1;
    }

    while (denom <= sqrt((double) n)){
        if (n % denom == 0){
            return 0;
        }
        ++denom;
    }
    return 1;
}
//-----------------------------------------------------------

/*
Function: findPrimesTill()

Finds all primes up to given number, n,
and returns them collected in array.
*/
void findPrimesTill(int* primes, unsigned int size, unsigned int upperBound{
    unsigned int index = 0;
    //int* ptrToArray = primes;

    unsigned int i = 0;
    for (i = 2; i < upperBound; ++i){
        if (isPrime(i)){
            primes[index++] = i;

            if (index >= size){
                printf("realloc on i = %d.\n", i);
                break;
            }
        }
    }
}
//-----------------------------------------------------------

/*
Function: isSumOfTwoPrimes()

Checks if argument is a sum of two prime.
*/
char isSumOfTwoPrimes(unsigned int target, int* primes, unsigned int size){
    unsigned int i;
    unsigned int remainder;
    int* ptrToArray = primes;

    for (i = 0; i < size; ++i){
        if (ptrToArray[i] < target){
            remainder = target - ptrToArray[i];
        }else{
            break;
        }

        if (binarySearch(remainder, primes, size)){
            printf("%d = %d + %d", target, ptrToArray[i], remainder);
            return 1;
        }
    } 
    return 0;
}
//-----------------------------------------------------------

/*
isSumOfUniqueThreePrimes();

*/
char isSumOfUniqueThreePrimes(unsigned int target, int* primes, unsigned int size){
    unsigned int i;
    unsigned int j;
    unsigned int remainder;
    int* ptrToArray = primes;

    for (i = 0; i < size; ++i){
        for (j = 0; j < size; ++j){
            if (ptrToArray[i] + ptrToArray[j] < target){
                remainder = target - ptrToArray[i] - ptrToArray[j];
            }else{
                break;
            }
            // check uniqueness 
            if (ptrToArray[i] != ptrToArray[j] && ptrToArray[j] != remainder && ptrToArray[i] != remainder){
                if (binarySearch(remainder, primes, size)){
                    printf("%d = %d + %d + %d", target, ptrToArray[i], ptrToArray[j], remainder);
                    return 1;
                }
            }
        }
    }
    return 0;
}
//-----------------------------------------------------------

/*
isSumOfThreePrimes();

*/
char isSumOfThreePrimes(unsigned int target, int* primes, unsigned int size{
    unsigned int i;
    unsigned int j;
    unsigned int remainder;
    int* ptrToArray = primes;

    for (i = 0; i < size; ++i){
        for (j = 0; j < size; ++j){
            if (ptrToArray[i] + ptrToArray[j] < target){
                remainder = target - ptrToArray[i] - ptrToArray[j];
            }else{
                break;
            }

            if (binarySearch(remainder, primes, size)){
                printf("%d = %d + %d + %d", target, ptrToArray[i], ptrToArray[j], remainder);
                return 1;
            }
        }
    }
    return 0;
}
//-----------------------------------------------------------

/*
Functiop: isSumOfTheMostSixPrimes();

It could probably be a recursive function.

Complexity: O(n^6)
*/
char isSumOfTheMostSixPrimes(unsigned int target, int* primes, unsigned int size){
    int* ptrToArray = primes;
    unsigned int bound = 6;
    unsigned int i, j, k , l ,m ,n;
    unsigned int currentSum = 0;

    for (i = 0; i < size; ++i){
        unsigned int currentSum = ptrToArray[i];
        if (currentSum == target) return 1;
        else if (currentSum > target) break;
        for (j = 0; j < size; ++j){
            currentSum = ptrToArray[i] + ptrToArray[j];
            if (currentSum == target) return 1;
            else if (currentSum > target) break;
            for (k = 0; k < size; ++k){
                currentSum = ptrToArray[i] + ptrToArray[j] + ptrToArray[k];
                if (currentSum == target) return 1;
                else if (currentSum > target) break;
                for (l = 0; l < size; ++l){
                    currentSum = ptrToArray[i] + ptrToArray[j] + ptrToArray[k] + ptrToArray[l]; 
                    if (currentSum == target) return 1;
                    else if (currentSum > target) break;
                    for (m = 0; m < size; ++m){
                       currentSum = ptrToArray[i] + ptrToArray[j] + ptrToArray[k] + ptrToArray[l] + ptrToArray[m];
                        if (currentSum == target) return 1;
                        else if (currentSum > target) break;
                        for (n = 0; n < size; ++n){
                            currentSum = ptrToArray[i] + ptrToArray[j] + ptrToArray[k] + ptrToArray[l] + ptrToArray[m] +  + ptrToArray[n];
                            if (currentSum == target) return 1;
                            else if (currentSum > target) break;
                        }
                    }
                }
            }
        }
    }
    return 0;
}
//-----------------------------------------------------------

/*
Function: isDifferenceOfPrimes();

*/
char isDifferenceOfPrimes(unsigned int target, int* primes, unsigned int size){
    int* ptrToArray = primes;
    unsigned int i, j;

    for (i = 0; i < size - 1; ++i){
        for (j = i + 1; j < size; ++j){
            if (target == ptrToArray[j] - ptrToArray[i]){
                printf("%d = %d - %d", target, ptrToArray[j], ptrToArray[i]);
                return 1;
            }
        }
    }
    return 0;
}
#endif

GoldBach.h

#ifndef GOLDBACH_H
#define GOLDBACH_H

// probably all uppedBounds in the for loops could be doubled
/*
Function: First();

Test first hypothesis.
*/
void First(int* primes, unsigned int size, unsigned int upperBound){
    unsigned int even;

    for (even = 4; even <= upperBound; even += 2){
        if (isSumOfTwoPrimes(even, primes, size)){
            printf("\nFirst Goldback's hypothesis not disproved!\n");
        }else{
            printf("\n?Exception: %d\n", even);
        }
    }
} 
//-----------------------------------------------------------

/*
Function: Second();

Test first hypothesis.
*/
void Second(int* primes, unsigned int size, unsigned int upperBound){
    unsigned int natural;

    for (natural = 17; natural <= upperBound; ++natural){
        if (isSumOfUniqueThreePrimes(natural, primes, size)){
            printf("\nSecond Goldback's hypothesis not disproved!\n");
        }else{
            printf("\n?Exception:: %d\n", natural);
        }
    }
}
//-----------------------------------------------------------

/*
Function: Third()

*/
void Third(int* primes, unsigned int size, unsigned int upperBound){
    int* ptrToArray = primes;
    unsigned int integer;

    for (integer = 0; integer < upperBound; ++integer){
        if (isSumOfTheMostSixPrimes(integer, primes, size)){
            printf("\nThird Goldback's hypothesis not disproved!\n");
        }else{
            printf("\n?Exception:: %d\n", integer);
        }
    }
}
//-----------------------------------------------------------

/*
Function: Fourth()

*/
void Fourth(int* primes, unsigned int size, unsigned int upperBound){
    unsigned int odd;

    for (odd = 7; odd <= upperBound; odd += 2){
        if (isSumOfThreePrimes(odd, primes, size)){
            printf("\nFourth Goldback's hypothesis not disproved!\n");
        }else{
            printf("\n?Exception:: %d\n", odd);
        }
    }
}
//-----------------------------------------------------------

/*
Function: Fifth();

*/
void Fifth(int* primes, unsigned int size, unsigned int upperBound){
    unsigned int even;

    for (even = 2; even <= upperBound; even += 2){
        if(isDifferenceOfPrimes(even, primes, size)){
            printf("\nFifth  Goldback's hypothesis not disproved!\n");
        }else{
            printf("\n?Exception:: %d\n", even);
        }
    }
}
//-----------------------------------------------------------

/*
Function: testFirstGoldbackHypothesis(void)

*/
void testGoldbackHypothesis(void){
    // calculate size of array
    int error = MAX / 10; // uses adding of error rather than memory reallocation 
    int size = arraySize(MAX) + error;

    // allocate memory for the array storing the primes
    int* primes = 0;
    primes = (int*)malloc(sizeof(int) * size);

    // check allocation
    if (!primes){
        printf("Failed to allocate memory for array!\n");
    }
    initArray(primes, size, 0);

    findPrimesTill(primes, size, MAX);

    printArray(primes, size);

    // First(primes, size, MAX);
    // Second(primes, size, MAX);
    // Third(primes, size, MAX);
    // Fourth(primes, size, MAX);
    Fifth(primes, size, MAX);

    printf("\nup to the number: %d.\n", MAX);

    // free allocated memory
    free(primes);
}

#endif 

Questions:

  1. Would it be better if memory is reallocated for each prime outside of the current array size? (Currently, there are few garbage values at the end of the array.)

  2. Is the current approach of checking right, are there more effective algorithms?

  3. Is the code written according to the C coding standard?


1. Progamming = ++Algorithms

\$\endgroup\$
  • 1
    \$\begingroup\$ Regarding #1, you could a) implement a linked list b) use a trick similar to C++'s std::vector, which reallocates to twice its current size when it runs out of room. \$\endgroup\$ – Mateen Ulhaq Sep 2 '16 at 11:42
  • 2
    \$\begingroup\$ Regarding #3, you can check by setting your compiler to its strictest settings/turning on warnings. \$\endgroup\$ – Mateen Ulhaq Sep 2 '16 at 11:42
3
\$\begingroup\$

Would it be better if memory is reallocated for each prime outside of the current array size?

No. See following.

Is the current approach of checking right, are there more effective algorithms?

  1. char isPrime(int n){ is barely run-time efficient. Rather than seek the next prime with ++denom, maintain a prime list and use the next one. For code using only int, could use a bit accessed array uint8_t IsPrime[(INT_MAX-1)/8 + 1] (about 512M bytes) or something smaller based on MAX, populated with Sieve of Eratosthenes. The question becomes what is "efficient". Is that speed, memory usage, code space usage, source code terseness, small stack usage? OP did not specify - assume speed.

  2. A compiler may not recognize that sqrt() has no side effects and that n is constant, thus repetitively calling sqrt(). Call sqrt() once. Better to round() result too.

    //while (denom <= sqrt((double) n)){
    unsigned limit = (unsigned) round(sqrt(n));
    while (denom <= limit) {
    
        ...if (n % denom == 0){
    }
    
  3. Good use of unsigned int middle = first + (last - first) / 2; to avoid overflow issues - even though (first + last) / 2; appears faster, the latter can fail.

  4. Minor. Returning char rather than int/unsigned is rarely faster/less code as that type is usually the processor's "preferred" type. Return int or bool. Profile code if this optimization is in doubt.

    // char isSumOfTheMostSixPrimes(...
    int isSumOfTheMostSixPrimes(...
    

Is the code written according to the C coding standard?

  1. .h files are best for define and declarations. It is non-standard practice to put code in a .h file.

  2. HelperArray.h does not include needed .h files. .h files should not depend on the code that includes them to have included certain files. This file should include them.

    // add for `log()`, etc.
    #include <math.h>
    
    long int pi(int x){
      return x / (log((double) x) - 1);
    }
    
  3. MAX is used, but not defined in HelperArray.h. The .h file needs to 1) not depend on other non-included declarations or defines or 2) should error intelligently.

    #ifndef MAX
      #error Define `MAX` before including HelperArray.h
    #endif
    void printArray(int* primes, unsigned int size){
        unsigned int i;
        int* ptrToArray = primes;
        int fieldWidth = 1 + log10((double)MAX);
        ...
    
  4. Change of type without checking range - candidate bug. Similar unqualified type changes used. The loose use of int/unsigned permeates code.

    char binarySearch(unsigned int target, int* primes, unsigned int size){
        ...
        // What if `last` > INT_MAX?
        int last = size;  
    
  5. Invalid code. Likely missing )

    void findPrimesTill(int* primes, unsigned int size, unsigned int upperBound{
    char isSumOfThreePrimes(unsigned int target, int* primes, unsigned int size{
    
  6. Unused variables are OK, but why have them?

    char isSumOfTheMostSixPrimes(unsigned int target, int* primes, unsigned int size){
        // Unused
        unsigned int bound = 6;
        unsigned int currentSum = 0;
    
    void Third(int* primes, unsigned int size, unsigned int upperBound){
        // Unused
        int* ptrToArray = primes;
    
  7. Cast is OK, but not needed per the standard

    // primes = (int*)malloc(sizeof(int) * size);
    // Better 
    primes = malloc(sizeof(int) * size);
    // Even better: Avoid coding the wrong type.
    primes = malloc(sizeof *primes * size);
    
  8. Unclear why functions return type long. The C standard uses size_t as the Goldilocks type (not too narrow, not too wide) for array indexing and size.

    // long int pi(int x){
    // long int arraySize(int x){
    size_t pi(int x){
    size_t arraySize(int x){
    
\$\endgroup\$
  • \$\begingroup\$ Thank you for the great feedback! 3 and 13 are connected, if last is unsigned int in some cases results in overflow and strange behaviours. 4 my compiler doesn't support bool, but I should probably use int as returning type. I hope I didn't miss something useful in points 5 - 9. 10: implementations will be moved in .c files. 11: I was wondering if I should do with the missing libraries what it is done with MAX in 12, first check and then include/error message? 14: just wrong editing when copying. \$\endgroup\$ – Ziezi Sep 2 '16 at 15:24
  • 1
    \$\begingroup\$ @Ziezi Points 5-9 not posted - reserve for later additions if any. first, middle, last can all be unsigned with minor changes to code. with no over/underflow. Something like size_t low = 0; size_t highp1 = n; while (low < highp1) { size_t mid = low + (highp1 - low) / 2; if (values[mid] > value) low = mid + 1; else if (values[mid] < value) highp1 = mid; else return mid + offset; }. 11: Sounds good. 4: What compiler are you using? \$\endgroup\$ – chux Sep 2 '16 at 17:27
  • \$\begingroup\$ Strike the + offset part. \$\endgroup\$ – chux Sep 2 '16 at 17:45
  • \$\begingroup\$ So, for the Binary search algorithm, I'm replacing the upper bound of the range just with the middle, rather than the middle - 1. 4: I am using an "old" VC++2010 because for some reason VC++2015 couldn't install C++ project templates. \$\endgroup\$ – Ziezi Sep 2 '16 at 18:33
  • 1
    \$\begingroup\$ @Ziezi There exist on the Internet some nice stdbool.h, stdint.h, etc. files for various VC environments. That helps one code to the C99/C11 standard closer even when using older MS compilers. \$\endgroup\$ – chux Sep 2 '16 at 18:50
1
\$\begingroup\$

Not every number can be prime

    unsigned int i = 0;
    for (i = 2; i < upperBound; ++i){
        if (isPrime(i)){
            primes[index++] = i;

Your findPrimesTill checks every integer, even though we can trivially discard all even numbers greater than two.

    primes[index++] = 2;

    unsigned int i;
    for (i = 3; i < upperBound; i += 2) {
        if (isPrime(i)) {
            primes[index++] = i;

Less trivially we can discard all odd numbers divisible by three other than three.

    primes[index++] = 2;
    primes[index++] = 3;

    unsigned int increment = 4;
    for (unsigned int i = 5; i < upperBound; i += increment) {
        increment = 6 - increment;
        if (isPrime(i)) {
            primes[index++] = i;

The increment = 6 - increment alternates between 2 and 4. So 5, 7, 11 (skip 9, which is divisible by 3), 13, 17 (skip 15), etc.

I also moved the declaration of i into the for loop. While older C compilers didn't support this, modern versions would.

Similarly, you use trial division with every integer in isPrime. You could use the same tricks with it.

char isPrime(int n){
    if (n < 2){
        return 0;
    }

    if (n == 2 || n == 3) {
        return 1;
    }

    if (n % 2 || n % 3) {
        return 0;
    }

    unsigned int increment = 4;
    unsigned int upper_bound = (unsigned int)sqrt((double) n);
    for (unsigned int factor_candidate = 5; factor_candidate <= upper_bound; factor_candidate += increment) {
        if (n % factor_candidate == 0) {
            return 0;
        }

        increment = 6 - increment;
    }

    return 1;
}

I moved taking the square root out of the loop so that we don't have to rely on the compiler to optimize it out of the loop. Note that for small values of MAX, e.g. 100, it's quicker to multiply on each iteration of the loop than to take the square root once.

Use what you know

The truth is that we can actually do better than this isPrime function though. You are generating a list of all prime numbers in order. Use that to test if the number is prime.

#include <stdbool.h>

bool isPrime(int candidate, int *primes, unsigned int size) {
    for (int i = 0; i < size; ++i) {
        if (candidate % primes[i] != 0) {
            return false;
        }

        if (candidate / primes[i] <= primes[i]) {
            return true;
        }
    }

    return true;
}

Note that this is not a general isPrime function. It is very specific to this particular problem.

While the original C did not have a Boolean type, it has been supported via library since C99.

The candidate / primes[i] <= primes[i] checks if we've reached the square root yet. It's done in that odd way because many systems will calculate the quotient and remainder at the same time. That makes the division essentially free (since we had to check the remainder anyway). And we would do a comparison regardless. You can also put the reverse of that check in the for loop check, but it can be harder to read there. For this specific problem, you could even replace the i < size with candidate / primes[i] >= primes[i]. You'll always hit the latter before running out of primes.

Note that it should be >= in the for loop, as the logic is reversed. There it is saying when to continue rather than when to end. It should be = in both places, as one place it is happening before the remainder check and the other place it is happening after.

Use a sieve

But the most efficient way to find every prime number up to some upper bound is probably a sieve. The best known is the Sieve of Eratosthenes. Then you don't need an isPrime method at all. The sieve sieves out the composite (non-prime) numbers, leaving only the primes.

void generatePrimesTill(unsigned int upperBound, int **primes, int *count) {
    *count = 0;
    if (upperBound < 2) {
        return;
    }

    bool* composites = calloc(upperBound, sizeof *composites);
    if (!composites) {
        // panic
        exit(-1);
    }

    composites[0] = true;
    composites[1] = true;

    for (int i = 2; i < upperBound; i++) {
        if (!composites[i]) {
            (*count)++;
            // this is prime, so sieve out all multiples starting with the square
            for (int j = i * i; j < upperBound; j += i) {
                composites[j] = true;
            }
        }
    }

    *primes = malloc(*count * sizeof **primes);
    if (!*primes) {
        // panic
        exit(-1);
    }

    for (int i = 2, j = 0; j < *count && i < upperBound; ++i) {
        if (!composites[i]) {
            (*primes)[j] = i;
            ++j;
        }
    }

    free(composites);
}

I haven't tested this version, but it is the basic pattern for a sieve solution.

Note that this is called like

    generatePrimesTill(MAX, &primes, &size);

instead of

    findPrimesTill(primes, size, MAX);

This allows it to return values for primes and size. Originally you had to estimate size and allocate primes first. This will have the side effect of always allocating a correctly sized primes.

It might be better to return EXIT_FAILURE or similar when panicking.

\$\endgroup\$
  • \$\begingroup\$ Thank you for the feedback! 1: That is a neat way to optimise the isPrime function, the alternating incrementation is something that I will remember. The rest two algorithms are equally good, except the tricky pointers in the last one that need to be considered when calling. \$\endgroup\$ – Ziezi Sep 2 '16 at 18:58
  • \$\begingroup\$ Btw, have you tested the optimised version of isPrime in point 1, if am not mistaken there should be something wrong as I'm testing it against other two to check run time but the number of found primes is different only in the version described in point 1. In addition it counts the multiple of 5 as prime \$\endgroup\$ – Ziezi Sep 2 '16 at 19:49
  • \$\begingroup\$ @Ziezi The problem looks to be that I didn't change the initial number from 2 to 5. Updated. Still not tested. \$\endgroup\$ – mdfst13 Sep 2 '16 at 21:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.