/*
Calculate terms for OEIS sequence A036236: Least inverse of A015910: smallest
integer k > 0 such that 2^k mod k = n, or 0 if no such k exists.
This program prints a table, with rows 'n k', for 0 <= n < TABSIZE = 10001.
The program can be run multiple times, for increasing ranges of k, and the
outputs merged into a single table.
Compile on Linux as, e.g.: clang -O3 -fopenmp -o a036236 a036236.c -lgmp
Invoke on Linux as, e.g.: ./a036236 1 1000000 4 >a036236.txt
where the first two parameters are, respectively, lower and upper bounds
for k, and the third parameter is the number of threads OpenMP should use.
Created by Nick Hobson (n@nickhobson.com), Jan 29 2024.
*/

#include <stdio.h>
#include <stdlib.h>
#include <limits.h>
#include <errno.h>
#include <gmp.h>
#include <omp.h>

#define TABSIZE (1 + 10000)

static unsigned long lower_k;
static unsigned long upper_k;
static unsigned long num_threads;
static unsigned long (*tab)[TABSIZE];

// For each n < TABSIZE, print smallest k > 0 found, across all threads.
// Print -1 for those entries where k = 0; i.e., a(n) is unknown.
static void print_table(void)
{
    puts("0 1");
    puts("1 0");
    for (size_t n = 2; n < TABSIZE; n++) {
        unsigned long min_k = ULONG_MAX;
        for (size_t t = 0; t < num_threads; t++) {
            if (tab[t][n] != 0 && tab[t][n] < min_k) {
                min_k = tab[t][n];
            }
        }
        if (min_k == ULONG_MAX)
            printf("%zu -1\n", n);
        else
            printf("%zu %lu\n", n, min_k);
    }
}

// For each n < TABSIZE, store smallest k in the half-open interval
// [lower_k, upper_k), in this thread, for which 2^k (mod k) = n.
static void a036236(const unsigned long start)
{
    mpz_t a, k, n;
    mpz_init_set_ui(a, 2);
    mpz_inits(k, n, NULL);
    for (unsigned long k_ui = lower_k + start; k_ui < upper_k;
         k_ui += num_threads) {
        mpz_set_ui(k, k_ui);
        mpz_powm(n, a, k, k);  // n = 2^k (mod k)
        if (mpz_cmp_ui(n, TABSIZE) < 0) {
            unsigned long n_ui = mpz_get_ui(n);
            if (tab[start][n_ui] == 0) {
                tab[start][n_ui] = k_ui;
            }
        }
    }
    mpz_clears(a, k, n, NULL);
}

// Call a036236 in parallel for different starting values. The output will be
// merged in print_table.
int main(int argc, char *argv[])
{
    // No validation for input parameters, which should be positive integers.
    lower_k = (argc > 1) ? strtoul(argv[1], NULL, 10) : 1;
    upper_k = (argc > 2) ? strtoul(argv[2], NULL, 10) : 1000000000;
    num_threads = (argc > 3) ? strtoul(argv[3], NULL, 10) : 4;

    tab = calloc(num_threads, sizeof(*tab));
    if (!tab) {
        fprintf(stderr,
                "Error: cannot allocate table of size %lu * %zu bytes.\n",
                num_threads, sizeof(*tab));
        return EXIT_FAILURE;
    }

#pragma omp parallel for schedule(static) num_threads(num_threads)
    for (unsigned long start = 0; start < num_threads; start++) {
        a036236(start);
    }

    print_table();
    free(tab);
    return EXIT_SUCCESS;
}