Prime factorization by trial division.

primes.c

/**
 * @file primes.c
 * @author Jose Fernando Lopez Fernandez <[email protected]>
 * @brief Prime factorization using trial division algorithm.
 * @version 0.1
 * @date 2022-05-18
 * 
 * @copyright Copyright (c) 2022
 * 
 * Compile and link using the following command.
 * 
 *    gcc -std=c17 -Wall -Wextra -Wpedantic -D_GNU_SOURCE -O3 -march=native -fexpensive-optimizations -o primes primes.c -lgmp
 * 
 */

#include <stdbool.h>
#include <stddef.h>
#include <stdint.h>
#include <stdlib.h>
#include <stdio.h>

#include <gmp.h>

int main(int argc, char *argv[])
{
    if (argc == 1) {
        fprintf(stderr, "Argument(s) expected.\n");
        return EXIT_FAILURE;
    }

    /** Return two, since it is a prime number. */
    printf("2\n");

    /**
     * Initialize the upper limit of the range in which we
     * will look for prime numbers.
     */
    mpz_t limit;
    // TODO: Ensure that the limit is a positive number.
    mpz_init_set_str(limit, argv[1], 10);
    // TODO: Remove this or add a verbosity flag.
    //gmp_printf("Limit: %Zd\n", limit);

    /** Initialize the integer iterator. */
    mpz_t n;
    mpz_init_set_ui(n, 3);

    /**
     * Continue iterating over every integer in the given
     * range until we hit the limit.
     */
    while (mpz_cmp(n, limit) <= 0) {
        /** Initialize the counter variable. */
        mpz_t i;
        mpz_init_set_ui(i, 2);

        /** Initialize the loop termination variable. */
        mpz_t l;
        mpz_init(l);

        /** Iterate over every value from 2 to n - 1. */
        // TODO: Set the loop termination variable to one more than the square root of the current value of n.
        mpz_set(l, n);
        mpz_sub_ui(l, l, 1);

        /** Initialize the division result variable. */
        mpz_t r;
        mpz_init(r);

        /** Initially, we believe all numbers are prime. */
        bool is_prime = true;

        /**
         * Perform trial division on each number from 3 to
         * the square root of n, plus one.
         */
        while (mpz_cmp(i, l) < 0) {
            /** Perform trial division of n by i. */
            mpz_cdiv_r(r, n, i);

            /**
             * Check whether n is perfectly divisible by i.
             */
            if (mpz_cmp_ui(r, 0) == 0) {
                /**
                 * A prime number is, by definition,
                 * divisible only by itself and one, so if
                 * we hit this point, we know that this
                 * number is definitely not prime.
                 */
                is_prime = false;

                /**
                 * Once we've detected that the number is
                 * not prime, we can simply break out of
                 * the loop to avoid performing a lot of
                 * unnecessary computations.
                 */
                break;
            }

            /** Increment the loop counter. */
            // TODO: Increment the loop counter by two to only check odd numbers.
            mpz_add_ui(i, i, 1);
        }

        /**
         * Once we have finished the computations in the
         * loop, we print the number to standard output if
         * it's prime.
         */
        if (is_prime) {
            gmp_printf("%Zd\n", n);
        }

        /** Reset the division result variable. */
        mpz_clear(r);

        /** Reset the counter variable. */
        mpz_clear(i);

        /** Reset the loop counter limit variable. */
        mpz_clear(l);

        /** Increment the current number. */
        mpz_add_ui(n, n, 1);
    }

    return EXIT_SUCCESS;
}