/*
Calculate terms for OEIS sequence A155894: Least positive integer such that
a(n)! starts with n, both written in base 16.
This program prints a table, with columns 'n a(n)', for 1 <= n < 16^upper_n.
Compile on Linux as, e.g.: clang -O3 -o a155894 a155894.c -lgmp
Invoke on Linux as, e.g.: ./a155894 4
where the parameter is log_16 of the upper bound for n.
Created by Nick Hobson (n@nickhobson.com), Mar 04, 2024.
*/

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

#define LOG2_16 4

static void a155894(unsigned long *tab, const unsigned long upper_n,
                    const unsigned long upper_log16_n)
{
    mpz_t fac, q;
    mpz_init_set_ui(fac, 1);
    mpz_init(q);
    unsigned long c = 0;
    // Loop until elements 1..upper_n - 1 of the table are filled. The existence
    // of a table entry for each n was proven by John E. Maxfield in 1970.
    for (unsigned long n = 1; c < upper_n - 1; n++) {
        mpz_mul_ui(fac, fac, n);
        // Isolate the high upper_log16_n nibbles of fac.
        size_t sz = mpz_sizeinbase(fac, 16);
        mpz_tdiv_q_2exp(
            q, fac, (sz > upper_log16_n) ? LOG2_16 * (sz - upper_log16_n) : 0);
        unsigned long high_nibbles = mpz_get_ui(q);
        assert(high_nibbles < upper_n);
        while (high_nibbles != 0) {
            if (tab[high_nibbles] == 0) {
                tab[high_nibbles] = n;
                c++;
            }
            high_nibbles /= 16;
        }
    }
    mpz_clears(fac, q, NULL);
}

int main(int argc, char *argv[])
{
    // Minimal validation for input parameter, which should be an integer.
    unsigned long upper_log16_n = (argc > 1) ? strtoul(argv[1], NULL, 10) : 4;
    unsigned long upper_n = 1UL << (LOG2_16 * upper_log16_n);
    if (upper_log16_n >= sizeof(upper_n) * CHAR_BIT / LOG2_16) {
        fprintf(stderr,
                "Parameter upper_log16_n = %lu must be less than %zu.\n",
                upper_log16_n, sizeof(upper_log16_n) * CHAR_BIT / LOG2_16);
        return EXIT_FAILURE;
    }
    unsigned long *tab = calloc(upper_n, sizeof(*tab));
    if (!tab) {
        fprintf(stderr, "Cannot allocate %lu elements of size %zu.\n", upper_n,
                sizeof(*tab));
        return EXIT_FAILURE;
    }
    a155894(tab, upper_n, upper_log16_n);
    for (unsigned long n = 1; n < upper_n; n++) {
        printf("%lu %lu\n", n, tab[n]);
    }
    free(tab);
    return EXIT_SUCCESS;
}