/* Kaprekar map: sort digits into descending and ascending order and subtract; look for cycles. Compile with: gcc -std=gnu99 -g -O2; modify BASE to use for a different base in the range 2 to 62. */ #define BASE 10 #include #include #include #include #include #include #include #include #include #define NORETURN __attribute__((__noreturn__)) #define PRINTF_FORMAT(X, Y) __attribute__((__format__(__printf__, X, Y))) #define UNUSED __attribute__((__unused__)) #define MALLOC_LIKE __attribute__((__malloc__)) /* Name of the program (filled in by util_init). */ static const char *program_name = "(util_init not called)"; /* util_init - initialise the utility functions given argc and argv: fill in program_name. */ static void util_init(int argc UNUSED, char **argv) { if (argv[0] == NULL) { program_name = "(no program name available)"; } else { program_name = strrchr(argv[0], '/'); if (program_name != NULL) program_name++; else program_name = argv[0]; } } /* die - print an error message, preceded by "program_name: error: " and followed by a newline, and exit with error status. */ static NORETURN PRINTF_FORMAT(1, 2) void die(const char *format, ...) { va_list args; /* We do not check for errors here when writing error messages, since there's nothing useful we could do if we detected one. */ fprintf(stderr, "%s: error: ", program_name); va_start(args, format); vfprintf(stderr, format, args); putc('\n', stderr); va_end(args); exit(EXIT_FAILURE); } /* panic - print an error that should never happen and exit. */ static NORETURN PRINTF_FORMAT(1, 2) void panic(const char *format, ...) { va_list args; /* We do not check for errors here when writing error messages, since there's nothing useful we could do if we detected one. */ fprintf(stderr, "%s: panic (should never happen): ", program_name); va_start(args, format); vfprintf(stderr, format, args); putc('\n', stderr); va_end(args); exit(EXIT_FAILURE); } /* xvfprintf - print something (as vfprintf, but with file name argument for error messages), checking for errors. */ static PRINTF_FORMAT(3, 0) void xvfprintf(FILE *stream, const char *filename, const char *format, va_list args) { int ret; ret = vfprintf(stream, format, args); if (ret < 0) die("writing %s: %s", filename, strerror(errno)); } /* xprintf - print something (as printf), checking for errors. */ static PRINTF_FORMAT(1, 2) void xprintf(const char *format, ...) { va_list args; va_start(args, format); xvfprintf(stdout, "standard output", format, args); va_end(args); } /* xputc - print a character (as putc, but with file name argument for error messages), checking for errors. */ static void xputc(int c, FILE *stream, const char *filename) { if (putc(c, stream) != c) die("%s: error writing file: %s", filename, strerror(errno)); } /* xflush - flush a stream (as fflush, but with file name argument for error messages), checking for errors. */ static void xflush(FILE *stream, const char *filename) { fflush(stream); if (ferror(stream)) die("%s: error writing file: %s", filename, strerror(errno)); } /* xmalloc - allocate memory, checking for errors. */ static MALLOC_LIKE void * xmalloc(size_t size) { void *tmp; tmp = malloc(size); if (tmp == NULL) die("out of memory in malloc(%zu)", size); return tmp; } /* xrealloc - reallocate memory, checking for errors. */ static MALLOC_LIKE void * xrealloc(void *buf, size_t size) { void *tmp; tmp = realloc(buf, size); if (tmp == NULL) die("out of memory in realloc(%zu)", size); return tmp; } /* arithmetic_overflow - give an "arithmetic overflow" error. */ static NORETURN void arithmetic_overflow(void) { die("arithmetic overflow"); } /* size_mul - safely multiply size_t values. */ static size_t size_mul(size_t size1, size_t size2) { size_t ret = size1 * size2; if (__builtin_expect((size1 | size2) >= (((size_t)1) << (CHAR_BIT * sizeof(size_t) / 2)), 0)) { if (ret / size2 != size1) arithmetic_overflow(); } return ret; } /* size_muladd - safely compute size1 * size2 + size3. */ static size_t size_muladd(size_t size1, size_t size2, size_t size3) { size_t ret = size_mul(size1, size2) + size3; if (ret < size3) arithmetic_overflow(); return ret; } /* xrealloc2 - safely reallocate (size1 * size2) memory. */ static MALLOC_LIKE void * xrealloc2(void *buf, size_t size1, size_t size2) { size_t size = size_mul(size1, size2); return xrealloc(buf, size); } /* xstrtoimax - read a signed integer within the range of intmax_t, expecting no trailing junk. Second argument is a description of what is being read, for error messages. */ static intmax_t xstrtoimax(const char *s, const char *desc) { intmax_t ret; char *endptr; errno = 0; ret = strtoimax(s, &endptr, 0); if (errno != 0) die("%s value '%s': %s", desc, s, strerror(errno)); if (endptr == s || *endptr != 0) die("%s value '%s' is not numeric or has trailing junk", desc, s); return ret; } /* Print a base-BASE digit. */ static void print_base (int d) { int c; #if BASE <= 10 c = '0' + d; #elif BASE <= 36 c = (d < 10 ? '0' + d : 'A' - 10 + d); #elif BASE <= 62 c = (d < 10 ? '0' + d : d < 36 ? 'A' - 10 + d : 'a' - 36 + d); #else #error "Can't print for base over 62" #endif xputc (c, stdout, "standard output"); } /* The array IN_COUNTS counts the number of each digit in a number NDIGITS long. Compute the results of applying the map once and store it in the array OUT, in little-endian form, and the counts in OUT_COUNTS. Return the new number of digits. */ static int one_iteration (int ndigits, const int in_counts[BASE], int out[ndigits], int out_counts[BASE]) { memset (out, 0, ndigits * sizeof (int)); memset (out_counts, 0, BASE * sizeof (int)); int pos = in_counts[0]; for (int i = 1; i < BASE; i++) { int npos = pos + in_counts[i]; for (int j = pos; j < npos; j++) { out[j] += i; out[ndigits - 1 - j] -= i; } pos = npos; } int carry = 0; for (int i = 0; i < ndigits; i++) { /* This value will always be in the range [-BASE, BASE-1]. */ int j = out[i] + carry; if (j < 0) { carry = -1; j += BASE; } else carry = 0; out[i] = j; out_counts[j]++; } if (carry != 0) panic ("carry at end of iteration"); while (out[ndigits - 1] == 0 && ndigits > 1) { ndigits--; out_counts[0]--; } return ndigits; } /* The array IN_COUNTS counts the number of each digit in a number NDIGITS long that is a member of a cycle. Determine the length of the cycle and its least element and call CYCLE_FUNC. */ static void check_cycle (int ndigits, const int in_counts[BASE], void cycle_func (int cf_ndigits, int cycle_length, const int cycle_counts[BASE], const int cycle_digits[ndigits])) { int omin[ndigits], out[ndigits]; int c0[BASE], c1[BASE], cmin[BASE]; memcpy (c0, in_counts, BASE * sizeof (int)); int i; for (i = 1; ; i++) { int n = one_iteration (ndigits, c0, out, c1); if (n != ndigits) panic ("number of digits went down in cycle"); bool new_min = true; if (i > 1) { for (int j = 0; j < ndigits; j++) { if (out[ndigits - 1 - j] < omin[ndigits - 1 - j]) break; if (out[ndigits - 1 - j] > omin[ndigits - 1 - j]) { new_min = false; break; } } } if (new_min) { memcpy (omin, out, ndigits * sizeof (int)); memcpy (cmin, c1, BASE * sizeof (int)); } if (memcmp (in_counts, c1, BASE * sizeof (int)) == 0) break; memcpy (c0, c1, BASE * sizeof (int)); } cycle_func (ndigits, i, cmin, omin); } /* The array IN_COUNTS counts the number of each digit in a number NDIGITS long. Iterate the map until a cycle, or a number smaller than IN_COUNTS in lexicographic order of sorted digits or in total number of digits, is reached; if a cycle, call CYCLE_FUNC. If SMALLER_OK, do not stop until a cycle is reached. */ static void iterate_one_input (int ndigits, const int in_counts[BASE], bool smaller_ok, void cycle_func (int cf_ndigits, int cycle_length, const int cycle_counts[BASE], const int cycle_digits[ndigits])) { int out[ndigits]; int slow1[BASE], slow2[BASE], fast1[BASE], fast2[BASE]; int sndigits = ndigits, fndigits = ndigits; memcpy (slow1, in_counts, BASE * sizeof (int)); memcpy (fast1, in_counts, BASE * sizeof (int)); while (1) { int c; c = one_iteration (fndigits, fast1, out, fast2); if (smaller_ok) fndigits = c; else { if (c < fndigits) return; for (int i = 0; i < BASE; i++) { if (fast2[i] < in_counts[i]) break; if (fast2[i] > in_counts[i]) return; } } sndigits = one_iteration (sndigits, slow1, out, slow2); c = one_iteration (fndigits, fast2, out, fast1); if (smaller_ok) fndigits = c; else { if (c < fndigits) return; for (int i = 0; i < BASE; i++) { if (fast1[i] < in_counts[i]) break; if (fast1[i] > in_counts[i]) return; } } if (memcmp (fast1, slow2, BASE * sizeof (int)) == 0) { check_cycle (fndigits, fast1, cycle_func); return; } c = one_iteration (fndigits, fast1, out, fast2); if (smaller_ok) fndigits = c; else { if (c < fndigits) return; for (int i = 0; i < BASE; i++) { if (fast2[i] < in_counts[i]) break; if (fast2[i] > in_counts[i]) return; } } sndigits = one_iteration (sndigits, slow2, out, slow1); c = one_iteration (fndigits, fast2, out, fast1); if (smaller_ok) fndigits = c; else { if (c < fndigits) return; for (int i = 0; i < BASE; i++) { if (fast1[i] < in_counts[i]) break; if (fast1[i] > in_counts[i]) return; } } if (memcmp (fast1, slow1, BASE * sizeof (int)) == 0) { check_cycle (fndigits, fast1, cycle_func); return; } } } static bool ignore_13 = false; struct cycle { int ndigits; int len; int counts[BASE]; int least_found; struct cycle *next; }; static size_t number_of_cycles, number_in_cycles; struct cycle_length_count { int len; size_t count; struct cycle_length_count *next; }; static struct cycle_length_count *cycle_length_counts; #define HASHSIZE 1000000 static struct cycle *cycle_hash[HASHSIZE]; /* The array CYCLE_COUNTS gives the number of each digit in the least element of a cycle of length CYCLE_LENGTH; the least element's digits are in CYCLE_DIGITS. Report the cycle if not already seen. */ static void report_cycle_if_new (int ndigits, int cycle_length, const int cycle_counts[BASE], const int cycle_digits[ndigits]) { if (ignore_13 && (cycle_length == 1 || cycle_length == 3)) return; size_t hashval = 0; hashval = hashval * 31 + ndigits; hashval = hashval * 31 + cycle_length; for (int j = 0; j < BASE; j++) hashval = hashval * 31 + cycle_counts[j]; hashval %= HASHSIZE; for (struct cycle *cl = cycle_hash[hashval]; cl != NULL && cl->ndigits == ndigits; cl = cl->next) { if (cl->len == cycle_length && memcmp (cl->counts, cycle_counts, BASE * sizeof (int)) == 0) return; } struct cycle *cnew = xmalloc (sizeof (struct cycle)); cnew->ndigits = ndigits; cnew->len = cycle_length; memcpy (cnew->counts, cycle_counts, BASE * sizeof (int)); cnew->least_found = 0; cnew->next = cycle_hash[hashval]; cycle_hash[hashval] = cnew; number_of_cycles++; number_in_cycles += cycle_length; struct cycle_length_count **p; for (p = &cycle_length_counts; *p != NULL; p = &((*p)->next)) if ((*p)->len == cycle_length) { (*p)->count++; break; } if (*p == NULL) { *p = xmalloc (sizeof (struct cycle_length_count)); (*p)->len = cycle_length; (*p)->count = 1; } xprintf ("%d-cycle: ", cycle_length); for (int j = 0; j < ndigits; j++) print_base (cycle_digits[ndigits - 1 - j]); xprintf ("\n"); xflush (stdout, "standard output"); int c0[BASE], c1[BASE], digs[ndigits]; memcpy (c0, cycle_counts, BASE * sizeof (int)); for (int i = 1; i < cycle_length; i++) { int n = one_iteration (ndigits, c0, digs, c1); if (n != ndigits) panic ("number of digits went down in cycle"); xprintf ("%d-cycle-next: ", cycle_length); for (int j = 0; j < ndigits; j++) print_base (digs[ndigits - 1 - j]); xprintf ("\n"); xflush (stdout, "standard output"); memcpy (c0, c1, BASE * sizeof (int)); } } static int current_in_counts[BASE]; static int current_ndigits; /* The array CYCLE_COUNTS gives the number of each digit in the least element of a cycle of length CYCLE_LENGTH; the least element's digits are in CYCLE_DIGITS. Report the least number reaching this cycle if not already done. */ static void report_least_reaching_cycle_if_new (int ndigits, int cycle_length, const int cycle_counts[BASE], const int cycle_digits[ndigits]) { if (ndigits < current_ndigits) return; size_t hashval = 0; hashval = hashval * 31 + ndigits; hashval = hashval * 31 + cycle_length; for (int j = 0; j < BASE; j++) hashval = hashval * 31 + cycle_counts[j]; hashval %= HASHSIZE; for (struct cycle *cl = cycle_hash[hashval]; cl != NULL && cl->ndigits == ndigits; cl = cl->next) { if (cl->len == cycle_length && memcmp (cl->counts, cycle_counts, BASE * sizeof (int)) == 0) { if (cl->least_found) return; cl->least_found = 1; number_of_cycles--; xprintf ("Least: "); for (int j = 1; j < BASE; j++) if (current_in_counts[j]) { print_base (j); current_in_counts[j]--; break; } for (int j = 0; j < BASE; j++) while (current_in_counts[j]) { print_base (j); current_in_counts[j]--; } xprintf (" %d ", cycle_length); for (int j = 0; j < ndigits; j++) print_base (cycle_digits[ndigits - 1 - j]); xprintf ("\n"); xflush (stdout, "standard output"); return; } } panic ("least number reaching new cycle"); } /* Just report a cycle. */ static void just_report_cycle (int ndigits, int cycle_length, const int cycle_counts[BASE] UNUSED, const int cycle_digits[ndigits]) { xprintf ("%d-cycle: ", cycle_length); for (int j = 0; j < ndigits; j++) print_base (cycle_digits[ndigits - 1 - j]); xprintf ("\n"); xflush (stdout, "standard output"); } /* Check for cycles among numbers of NDIGITS digits, where USED digits have been used for the digits below NEXT; COUNTS has the counts so far, TOTAL is the total of the digits so far and DIFF is the total of the absolute values of the differences between the numbers of digits d and BASE-1-d (for 1 <= d <= BASE-2), which must not end up exceeding 4. */ static void check_for_digits_recursive (int ndigits, int total, int diff, int used, int next, int counts[BASE]) { if ((used == ndigits || next == BASE - 1) && ((total % (BASE - 1)) != 0)) return; if (next == BASE - 1) { counts[next] = ndigits - used; iterate_one_input (ndigits, counts, false, report_cycle_if_new); return; } if (next <= ((BASE - 1) / 2 + 1)) { int half = used; if ((BASE - 1) % 2 == 0 && next == ((BASE - 1) / 2 + 1)) half -= counts[(BASE - 1) / 2]; if (used + half - 4 > ndigits) return; } for (int i = ndigits - used; i >= 0; i--) { int ndiff = diff; if (next > (BASE - 1) / 2) { ndiff += abs (i - counts[BASE - 1 - next]); if (ndiff > 4) continue; } counts[next] = i; check_for_digits_recursive (ndigits, total + i * next, ndiff, used + i, next + 1, counts); } } /* Check for cycles among numbers of NDIGITS digits. */ static void check_for_digits (int ndigits) { int counts[BASE]; xprintf ("Considering %d digits.\n", ndigits); xflush (stdout, "standard output"); check_for_digits_recursive (ndigits, 0, 0, 0, 0, counts); if (!ignore_13) { xprintf ("Number-of-cycles %d %zu\n", ndigits, number_of_cycles); xprintf ("Number-in-cycles %d %zu\n", ndigits, number_in_cycles); } for (struct cycle_length_count *p = cycle_length_counts; p != NULL; p = p->next) { xprintf ("Cycle-length-count %d %d %zu\n", ndigits, p->len, p->count); p->count = 0; } xflush (stdout, "standard output"); } /* Find the least number reaching each cycle of NDIGITS digits, where USED digits have been used for the digits below NEXT; COUNTS has the counts so far. */ static void find_least_sources_recursive (int ndigits, int used, int next, int counts[BASE]) { if (next == BASE - 1) { counts[next] = ndigits - used; memcpy (current_in_counts, counts, BASE * sizeof (int)); iterate_one_input (ndigits, counts, true, report_least_reaching_cycle_if_new); return; } for (int i = ndigits - used; i >= 0; i--) { counts[next] = i; find_least_sources_recursive (ndigits, used + i, next + 1, counts); if (number_of_cycles == 0) return; } } /* Find the least number reaching each cycle of NDIGITS digits. */ static void find_least_sources (int ndigits) { int counts[BASE]; current_ndigits = ndigits; if (number_of_cycles == 0) return; if (ndigits == 1) { memset (counts, 0, BASE * sizeof (int)); counts[0] = 1; memcpy (current_in_counts, counts, BASE * sizeof (int)); iterate_one_input (ndigits, counts, true, report_least_reaching_cycle_if_new); if (number_of_cycles == 0) return; } for (int i = 1; i < BASE; i++) /* Least nonzero digit is i; there are j 0s. */ for (int j = ndigits - 1; j >= 0; j--) { memset (counts, 0, BASE * sizeof (int)); counts[0] = j; if (i == BASE - 1) { counts[i] = ndigits - j; memcpy (current_in_counts, counts, BASE * sizeof (int)); iterate_one_input (ndigits, counts, true, report_least_reaching_cycle_if_new); if (number_of_cycles == 0) return; } else for (int k = ndigits - j; k >= 1; k--) { counts[i] = k; find_least_sources_recursive (ndigits, j + k, i + 1, counts); if (number_of_cycles == 0) return; } } } /* Determine the result of the sequence starting with an NDIGITS number whose first digit is 1, whose last is D and whose other digits are 0. */ static void check_10d (int ndigits, int d) { int counts[BASE] = { 0 }; counts[0] = ndigits - 2; counts[1] = 1; counts[d]++; iterate_one_input (ndigits, counts, true, just_report_cycle); } /* Apply the map once to a value represented as an int, returning the result as an int. */ static int one_iteration_int (int value) { int ndigits = 0; int counts[BASE] = { 0 }; while (value != 0) { int d = value % BASE; value /= BASE; counts[d]++; ndigits++; } if (ndigits == 0) { counts[0]++; ndigits++; } int out[ndigits]; int out_counts[BASE]; ndigits = one_iteration (ndigits, counts, out, out_counts); int ret = 0; for (int i = ndigits - 1; i >= 0; i--) ret = ret * BASE + out[i]; return ret; } /* Find the trajectory of a value represented as an int, reporting the (all in decimal) immediate image, the first element of a cycle reached, the length of the preperiod and the length of the preperiod+cycle. */ static void trajectory_int (int value) { static int *trajectory_array = NULL; static size_t trajectory_alloc = 0; int value_succ = one_iteration_int (value); if (trajectory_alloc < 1) { trajectory_alloc = 16; trajectory_array = xrealloc2 (trajectory_array, trajectory_alloc, sizeof (int)); } trajectory_array[0] = value; size_t count = 1; size_t cycle_length = 0; while (1) { int s = one_iteration_int (trajectory_array[count - 1]); for (size_t j = 1; j <= count; j++) if (s == trajectory_array[count - j]) { cycle_length = j; break; } if (cycle_length) break; if (trajectory_alloc <= count) { trajectory_alloc = size_muladd (trajectory_alloc, 2, 16); trajectory_array = xrealloc2 (trajectory_array, trajectory_alloc, sizeof (int)); } trajectory_array[count++] = s; } xprintf ("%d %d %d %zu %zu\n", value, value_succ, trajectory_array[count - cycle_length], count - cycle_length, count); xflush (stdout, "standard output"); } /* Read an int argument ARG (description DESC) within the range MIN to MAX inclusive. */ static int read_int_arg (const char *arg, const char *desc, int min, int max) { intmax_t val = xstrtoimax (arg, desc); if (val < min || val > max) die ("%s value '%s' outside range [%d, %d]", desc, arg, min, max); return val; } int main (int argc, char **argv) { util_init (argc, argv); int dig_min = 1, dig_max = INT_MAX / BASE; int ndigits; bool compute_10d_trajectories = false; bool find_least = false; bool compute_range_trajectories = false; int dval = 0; int ret; while ((ret = getopt (argc, argv, "iltm")) != -1) switch (ret) { case 'i': ignore_13 = true; break; case 'l': find_least = true; break; case 't': compute_10d_trajectories = true; break; case 'm': compute_range_trajectories = true; break; case '?': /* Error message already printed. */ exit (EXIT_FAILURE); default: panic ("bad return %d from getopt", ret); } if (ignore_13 + find_least + compute_10d_trajectories + compute_range_trajectories > 1) die ("incompatible options passed"); if (compute_range_trajectories) { int val_min = 0; int val_max = BASE; while (val_max < INT_MAX / BASE) val_max *= BASE; val_max--; if (argc > optind) { val_min = read_int_arg (argv[optind], "least value", val_min, val_max); optind++; } if (argc > optind) { val_max = read_int_arg (argv[optind], "greatest value", val_min, val_max); optind++; } for (int i = val_min; i <= val_max; i++) trajectory_int (i); exit (EXIT_SUCCESS); } if (compute_10d_trajectories) { dig_min = 2; if (argc > optind) { dval = read_int_arg (argv[optind], "last digit", 0, BASE - 1); optind++; } } if (argc > optind) { dig_min = read_int_arg (argv[optind], "minimum number of digits", dig_min, dig_max); optind++; } if (argc > optind) { dig_max = read_int_arg (argv[optind], "maximum number of digits", dig_min, dig_max); optind++; } if (argc > optind) die ("stray command-line argument"); for (ndigits = dig_min; ndigits <= dig_max; ndigits++) { if (compute_10d_trajectories) check_10d (ndigits, dval); else { number_of_cycles = 0; number_in_cycles = 0; check_for_digits (ndigits); if (find_least) find_least_sources (ndigits); } } exit (EXIT_SUCCESS); }