/* 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; }