/*
Calculate terms for OEIS sequence A155891: Least positive integer such that
a(n)! starts with n, both written in binary.
This program prints a table, with columns 'n a(n)', for 1 <= n < 2^upper_n.
Compile on Linux as, e.g.: clang -O3 -o a155891 a155891.c -lgmp
Invoke on Linux as, e.g.: ./a155891 16
where the parameter is log_2 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>

static void a155891(unsigned long *tab, const unsigned long upper_n,
                    const unsigned long upper_log2_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_log2_n bits of fac.
        size_t sz = mpz_sizeinbase(fac, 2);
        mpz_tdiv_q_2exp(q, fac, (sz > upper_log2_n) ? sz - upper_log2_n : 0);
        unsigned long high_bits = mpz_get_ui(q);
        assert(high_bits < upper_n);
        while (high_bits != 0) {
            if (tab[high_bits] == 0) {
                tab[high_bits] = n;
                c++;
            }
            high_bits /= 2;
        }
    }
    mpz_clears(fac, q, NULL);
}

int main(int argc, char *argv[])
{
    // Minimal validation for input parameter, which should be an integer.
    unsigned long upper_log2_n = (argc > 1) ? strtoul(argv[1], NULL, 10) : 16;
    unsigned long upper_n = 1UL << upper_log2_n;
    if (upper_log2_n >= sizeof(upper_n) * CHAR_BIT) {
        fprintf(stderr, "Parameter upper_log2_n = %lu must be less than %zu.\n",
                upper_log2_n, sizeof(upper_log2_n) * CHAR_BIT);
        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;
    }
    a155891(tab, upper_n, upper_log2_n);
    for (unsigned long n = 1; n < upper_n; n++) {
        printf("%lu %lu\n", n, tab[n]);
    }
    free(tab);
    return EXIT_SUCCESS;
}