Integer factoring in C using GMP for numbers

Question

This code performs integer factoring using the GMP library, and can also use OpenMP for parallel computation.

I compile it using: gcc -lgmp -fopenmp and whatever optimization flags I wish to try.

#include <stdio.h>
#include <gmp.h>

void next_prime(MP_INT *x) {
    MP_INT r;
    mpz_init(&r);
    if (!mpz_probab_prime_p(x, 128u)) {
        mpz_mod_ui(&r, x, 2ul);
        if (mpz_cmp_ui(&r, 0ul) == 0)
            mpz_add_ui(x, x, 1ul);
        else
            mpz_add_ui(x, x, 2ul);
        while (!mpz_probab_prime_p(x, 128u)) {
            mpz_add_ui(x, x, 2ul);
        }
    }
    mpz_clear(&r);
}

void prev_prime(MP_INT *x) {
    MP_INT r;
    mpz_init(&r);
    if (!mpz_probab_prime_p(x, 128u)) {
        mpz_mod_ui(&r, x, 2ul);
        if (mpz_cmp_ui(&r, 0ul) == 0)
            mpz_sub_ui(x, x, 1ul);
        else
            mpz_sub_ui(x, x, 2ul);
        while (!mpz_probab_prime_p(x, 128u)) {
            mpz_sub_ui(x, x, 2ul);
        }
    }
    mpz_clear(&r);
}

unsigned long mp_log(MP_INT *x, MP_INT *base) {
    MP_INT r;
    mpz_init(&r);
    unsigned long lg = 0;
    mpz_set(&r, x);
    if (base != NULL) {
        while (mpz_cmp_ui(&r, 0ul) != 0) {
            mpz_div(&r, &r, base);
            lg += 1;
        }
    }
    else {
        while (mpz_cmp_ui(&r, 0ul) != 0) {
            mpz_div_ui(&r, &r, 2ul);
            lg += 1;
        }
    }
    mpz_clear(&r);
    return lg;
}

void factor(MP_INT *x, MP_INT *n, MP_INT *base) {
    /* Initialize variables. */
    MP_INT s;
    MP_INT r;
    mpz_init(&s);
    mpz_init(&r);
    /* Compute stuff. */
    mpz_set(x, n);
    if (mpz_probab_prime_p(n, 128u)) {
        goto end;
    }
    mpz_set(&s, base);
    mpz_mod(&r, n, &s);
    if (mpz_cmp_ui(&r, 0ul) == 0) {
        mpz_set(x, &s);
        goto end;
    }
    unsigned long k = mp_log(n, base);
    unsigned long a = mp_log(n, NULL);
    if (mpz_cmp(x, n) == 0) {
        /* Try different powers of n. */
        MP_INT p;
        mpz_init(&p);
        mpz_pow_ui(&p, &s, k);
        unsigned long i = 1ul;
        while (mpz_cmp(x, n) == 0) {
            #pragma omp parallel for
            for (unsigned long j = i; j < i + 2*a*a; ++j) {
                if (mpz_cmp(x, n) == 0) {
                    MP_INT g;
                    mpz_init(&g);
                    mpz_powm_ui(&g, n, j, &p);
                    mpz_gcd(&g, &g, n);
                    if (mpz_cmp_ui(&g, 1ul) != 0 && mpz_cmp(&g, n) != 0) {
                        mpz_set(x, &g);
                    }
                    mpz_clear(&g);
                }
            }
            i += 2*a*a;
        }
        mpz_clear(&p);
    }
    /* Clear variables. */
    end:
    mpz_clear(&s);
    mpz_clear(&r);
}

int main() {
    char str [65536];
    char bas [65536];
    while (fscanf(stdin, "%s %s", str, bas) != EOF) {
        MP_INT n;
        MP_INT base;
        mpz_init(&n);
        mpz_init(&base);
        mpz_set_str(&n, str, 10u);
        mpz_set_str(&base, bas, 10u);
        mpz_out_str(stdout, 10u, &n);
        int p = mpz_probab_prime_p(&n, 128u);
        printf(" %d ", p);
        mpz_out_str(stdout, 10u, &base);
        printf("\n");
        while (!p) {
            MP_INT x;
            mpz_init(&x);
            factor(&x, &n, &base);
            mpz_out_str(stdout, 10u, &x);
            printf(" ");
            mpz_div(&n, &n, &x);
            p = mpz_probab_prime_p(&n, 128u);
            mpz_clear(&x);
        }
        mpz_out_str(stdout, 10u, &n);
        printf("\n");
        mpz_clear(&n);
        mpz_clear(&base);
    }
    return 0;
}

Answer 1

You can improve the readability of your code, your indentation is great but you have magic numbers and the variable names aren't very clear. You probably wouldn't need comments if your variable names were clear and you had symbolic constants rather than numbers.

Magic Numbers:

You have a lot of raw numbers : 65536, 10u, 128u, 1ul, 2ul ...

For example you can have a symbolic constant CHARACTER_BUFFER_SIZE for 65536. Using the gcc compiler you would write your constant as

#define CHARACTER_BUFFER_SIZE 65536 

When compiling with g++ you would use

static const CHARACTER_BUFFER_SIZE = 65536;

this would change the first 2 lines in main to

    char str[CHARACTER_BUFFER_SIZE];
    char bas[CHARACTER_BUFFER_SIZE];

If you use this number often, it makes it easier to change the value by just changing it where you define the constant. It also makes the code easier to read. The standard for symbolic constants in C and C++ is all caps with '_' between words.

Function Names: Underscore in function names is a lot less common than it used to be, Camel Case is more common

next_prime => nextPrime
prev_prime => prevPrime

Variable Names:

The base variable name is almost understandable, but what kind of base, does it provide a starting point for the function, or is it a numeric base like Octal, Hex or Decimal?

It's really not clear from the code what x is.

In main() what are the variables str (almost obvious), bas, n, p, base?

Future Flexibility:

You main could be a little more flexible, currently you can only redirect input and a output, you could use command switches in your to select input and output files as well. In main it is possible to add a couple of FILE * variable to make it a little more general

    FILE *inputFile = stdin;
    FILE *outputFile = stdout;

    while(fscanf(inputFile, "%s %s", str, bas) != EOF) {
        ...
        mpz_out_str(outputFile, 10u, &base);

Answer 2

This code is 137 lines. Of which, only 5 are empty blank lines. A little bit of breathing room can do a lot for the readability of your code.

Your use of braces is also inconsistent. You omit them for one-line if and else statements, but maintain them for one line while loops. You really shouldn't consider them optional--you should implement them every where.

So keeping those two things in mind, here's what a more readable implementation of your first method looks like:

void next_prime(MP_INT *x) {
    MP_INT r;
    mpz_init(&r);

    if (!mpz_probab_prime_p(x, 128u)) {
        mpz_mod_ui(&r, x, 2ul);

        if (mpz_cmp_ui(&r, 0ul) == 0) {
            mpz_add_ui(x, x, 1ul);
        }
        else {
            mpz_add_ui(x, x, 2ul);
        }

        while (!mpz_probab_prime_p(x, 128u)) {
            mpz_add_ui(x, x, 2ul);
        }
    }

    mpz_clear(&r);
}

There's still more work to be done though. This code is still hard to read, but it's mostly a result of unclear variable and function names.

I also can't see any good reason why this method takes an integer by reference and returns void. Taking a variable by reference in order to set it is a useful tactic when you need to return multiple pieces of information (for example, an error code or something). But if your return type is void... why don't you just take the integer's value and return the "next prime" value that you calculate? When I read the function name "next prime", that's exactly what I expect to see... a function which takes an integer argument, and returns the same type, representing the first prime number appearing after the number I passed in.