/* Kevin Ryde, July 2024 Usage: ./a.out [-v] [n] [n-n] This is some C code using the primesieve library, and best in a 64 bit build, to calculate terms of A065851(n) = maximum number of distinct n-digit primes which can be formed by permuting some list of n decimal digits A373179(n) = smallest n-digit integer whose digit permutations attain that maximum The method is a search through all primes of n digits, with a boost from some table lookups to help digit twiddling. Command Line ------------ The calculation is made for individual n or ranges of n on the command line. Output is a free-form report ./a.out 6 => n=6 max count 148 at digits 1123479 which means A065851(6) = 148 and A373179(6) = 123479. Optional command line parameter "-v" prints more information at the end of each calculation. Save and Resume --------------- The current n, counts, and next prime, are saved to a file saved-state.data at 5 minute intervals and on kill (SIGINT, SIGTERM, SIGHUP). The search of that n can be resumed by ./a.out resume This is intended for a long run attempting say n=14 or n=15 and which might have to be restarted. The file format is a dump of an internal structure and array so is architecture dependent (type sizes, structure padding, endianness). Progress information is written to a text file progress.txt. Compile Time ------------ The compile-time limits on n attempted are n <= MAXN = 16 index_table[] array size n <= 19 so primes fit primesieve uint64_t n <= DIGS_ELEM_MAX = 31 digs 6 bit field MAXN can be increased as desired, but the expectation would be that say n=17 is too big to iterate primes explicitly. DIGS_ELEM_MAX is the way a "digs" is represented and is not easily changed, but is already beyond the largest prime which primesieve can make. Speed ----- On a 3GHz CPU at the time of writing, n=14 took about 38 hours. n=15 would be an estimated at least 400 hours and likely more. Data type "count_t" is 64 bits anticipating that an n=15 count might exceed 32 bits (which suffices for n <= 14). Not sure whether the smaller type might be a little faster by fitting more into L1/L2 cache. Algorithm --------- Primes p which are n digits long are 10^(n-1) <= p < 10^n A006879(n) many p Each p is turned into a "digs", digs[d] = how many times digit d occurs in p, for 0 <= d <= 9 digs = type uint64_t with counts in 6 bits each A digs becomes an index for the state->counts[] array, index = digs_to_index(digs) counts[index] = how many times "digs" occurred 0 <= index < counts_len where counts_len = binomial(n+9,9) After all primes, max_count = A065851(n) = maximum count in counts[] A373179(n) = smallest digs_to_smallest(digs) for those digs attaining max_count It's possible that max_count is attained by more than one digs. This happens for n = 1..4. digs_to_smallest() is the smallest number of the given digs, and then A373179 is the smallest over all the digs which attain max_count. Digs to Index ------------- The index for "digs" is by reckoning digs as a string of digits in ascending order, and sorting those strings in lexicographically increasing order, digits index ----------- -------------- 0,0,...,0,0 0 0,0,...,0,1 1 ... ... 9,9,...,9,9 counts_len - 1 The index is formed by index = Sum index_table[d][i] d=0..8 where i = how many digits > d occur in digs and index_table[d][i] = binomial (i-1 + 9-d, 9-d) i is the number of positions remaining for further digits > d. The binomial is how many digit combinations are skipped due to the next position being >= d+1, rather than = d. If the next position were d there would be i-1 further positions to fill with ascending digits >= d. The binomial is i-1 places which have digits and additional 9-d places anywhere among them where step to next digit. These steps result in digit values d to 9 inclusive, zero or more of each, and total i-1 of them. The sum goes only to d=8, since at d=9 there are no digits > d and no combinations to skip. The sole choice is all 9s there on. At d=8, index_table[8][i] = binomial(i,1) = i itself, and the code simply adds that without a table lookup. In digs_to_index(), the n variable starts as the total number of digits and decreases by digs[d] each time so is the "i". Digs from Index --------------- digs_from_index() is the inverse operation, from index to digs. for d=0..9 find smallest i>=0 for which index >= index_table[d][i] digs[d] = n - i index -= index_table[d][i] n = i (remaining) [end with index = 0] When a small number of digit d occur, there's many combinations skipped (those where d occurred more times) so index_table[d][i] is big. When "index" is big enough to have skipped those, it did skip those. This rule for taking a greedy amount of skip, and hence fewest d, recovers the original digs uniquely. At the largest table entry with index >= index_table[d][i], that i must be taken. Any delay means not enough to get index down to 0, since index_table[d][i] = 1 + Sum index_table[e][i-1] e = d..8 This equality is an identity in binomials. The sum on the right is the largest possible of the d..8, and yet the total is still 1 short of the LHS index_table[d][i]. Speed ----- The key operations are the primesieve iterator and then p_to_digs_with_cache() and digs_to_index(). A 3GHz CPU at the time of writing was about 20 to 30 million primes/second near n=13. But the rate decreases with n, and the number of primes (A006879) grows about 9x for each n. */ /* Set WANT_ASSERT to 1 for various assert() self-checks. Set WANT_DEBUG to 1 for verbose diagnostic prints. */ #define WANT_ASSERT 1 #define WANT_DEBUG 0 #if ! WANT_ASSERT #define NDEBUG #endif #include #include #include #include #include #include #include #include #include #include #include #include #include #include /* ---------------------------------------------------------------------- */ /* Generic Helpers */ #if WANT_DEBUG #define DEBUG(expr) do { expr; } while (0) #else #define DEBUG(expr) do { } while (0) #endif #define numberof(array) (sizeof(array)/sizeof((array)[0])) #define SIZE_T_MAX (~ (size_t) 0) #ifdef __GNUC__ #define LIKELY(cond) __builtin_expect((cond) != 0, 1) #define UNLIKELY(cond) __builtin_expect((cond) != 0, 0) #define ATTRIBUTE_PRINTF __attribute__ ((format (printf, 1, 2))) #else #define LIKELY(cond) (cond) #define UNLIKELY(cond) (cond) #define ATTRIBUTE_PRINTF #endif static noreturn ATTRIBUTE_PRINTF void error (const char *format, ...) { va_list ap; va_start (ap, format); vfprintf (stderr, format, ap); va_end (ap); exit(1); } /* Return x+y or die if that would overflow size_t. */ static size_t size_add_or_die (size_t x, size_t y) { if (UNLIKELY (x > SIZE_T_MAX - y)) error("Add %zu + %zu overflows size_t (of %zu bytes)\n", x, y, sizeof(size_t)); return (x+y); } /* Return product x*y or die if that would overflow size_t. */ static size_t size_multiply_or_die (size_t x, size_t y) { if (UNLIKELY (y > SIZE_T_MAX / x)) error("Multiply %zu * %zu overflows size_t (of %zu bytes)\n", x, y, sizeof(size_t)); return (x*y); } /* with ptr=NULL allowed for initial malloc() too */ static void * realloc_or_die(void *ptr, size_t size) { ptr = (ptr==NULL ? malloc(size) : realloc(ptr, size)); if (UNLIKELY (ptr == NULL)) error("Cannot realloc() to %zu bytes\n", size); return (ptr); } static FILE * fopen_or_die (const char *filename, const char *mode) { FILE *fp = fopen(filename,mode); if (fp==NULL) error("Cannot %s %s\n", mode[0] == 'r' ? "open" : "create", filename); return (fp); } static void ferror_die (FILE *fp, const char *filename, const char *action) { if (ferror(fp)) error("Error %s %s: %s\n", action, filename, strerror(errno)); } static void fseeko_or_die (FILE *fp, off_t offset, int whence, const char *filename) { if (UNLIKELY(fseeko (fp, offset, whence) != 0)) error("Cannot fseeko in %s offset %"PRId64": %s\n", filename, (int64_t) offset, strerror(errno)); } static off_t ftello_or_die (FILE *fp, const char *filename) { off_t pos = ftello(fp); if (pos == -1) error("Cannot ftello in %s: %s\n", filename, strerror(errno)); return pos; } /* size of the file open at "fp" */ static off_t fsize_or_die (FILE *fp, const char *filename) { off_t old_pos = ftello_or_die (fp, filename); fseeko_or_die (fp, 0, SEEK_END, filename); off_t len = ftello_or_die (fp, filename); fseeko_or_die (fp, old_pos, SEEK_SET, filename); return (len); } static void fread_or_die (void *ptr, size_t size, size_t num, FILE *fp, const char *filename) { size_t got = fread (ptr, size,num, fp); ferror_die (fp, filename, "reading"); if (got != num) error("File %s too small: want %zu records of %zu bytes, got %zu\n %s", filename, num, size, got, strerror(errno)); } static void fwrite_or_die (const void *ptr, size_t size, size_t num, FILE *fp, const char *filename) { if (fwrite(ptr,size,num,fp) != num) error("Error writing %s, size %zu num %zu bytes: %s", filename, num, size, strerror(errno)); } static void fclose_or_die (FILE *fp, const char *filename) { if (fclose(fp) != 0) error("Error closing %s\n", filename); } static void rename_or_die (const char *old_filename, const char *new_filename) { if (rename(old_filename, new_filename) != 0) error("Error renaming %s to %s: %s\n", old_filename, new_filename, strerror(errno)); } static void sigaction_or_die (int signum, const struct sigaction *act, struct sigaction *oldact) { if (sigaction (signum, act, oldact) != 0) error("Cannot sigaction signum=%d: %s\n", signum, strerror(errno)); } static void setitimer_or_die (int which, const struct itimerval *new_value, struct itimerval *old_value) { if (setitimer(which, new_value, old_value) != 0) error("Cannot setitimer %d: %s\n", which, strerror(errno)); } /* CPU time consumed by this process so far, in seconds */ static double cputime_increment (void) { static double prev = 0; struct timespec ts; if (clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &ts) != 0) error("Cannot clock_gettime() CPUTIME: %s\n", strerror(errno)); double t = ts.tv_sec + ts.tv_nsec/1e9; double increment = t - prev; prev = t; return (increment); } /* -------------------------------------------------------------------------- */ /* uint64_t arithmetic */ /* Bit mask comprising lowest n bits of uint64_t. Must have 0 <= n <= 64. */ #define UINT64_LOW_MASK(n) \ ((n)==64 ? UINT64_MAX : (~ (UINT64_MAX << (n)))) static uint64_t uint64_mul (uint64_t x, uint64_t y) { if (UNLIKELY(x > UINT64_MAX / y)) error("uint64_mul() overflow\n"); return (x*y); } static const uint64_t uint64_pow10_table[] = { 1, 10, 100, 1000, 10000, 100000, 1000000, 10000000, 100000000, 1000000000, 10000000000, 100000000000, 1000000000000, 10000000000000, 100000000000000, 1000000000000000, 10000000000000000, 100000000000000000, 1000000000000000000 }; static uint64_t uint64_pow10 (int n) { if (! (0 <= n && n < numberof(uint64_pow10_table))) error("uint64_pow10() n=%d out of range\n", n); return (uint64_pow10_table[n]); } static uint64_t uint64_div_exact (uint64_t n, uint64_t d) { if (n%d != 0) error("oops, %"PRIu64" not divisible by %"PRIu64"\n",n,d); return (n/d); } static uint64_t uint64_binomial (int n, int k) { if (! (0 <= k && k <= n)) error("uint64_binomial() is for 0 <= k <= n\n"); uint64_t ret = 1; int i; for (i=1; i<=k; i++) { ret = uint64_mul(ret, n-k + i); ret = uint64_div_exact(ret, i); } return (ret); } static size_t size_binomial (int n, int k) { uint64_t b = uint64_binomial(n,k); if (b > SIZE_T_MAX) error("size_binomial() too big\n"); return (b); } /* -------------------------------------------------------------------------- */ static int option_verbose = 0; #define MAXN 16 /* for holding count of how many primes have a given digs */ typedef uint64_t count_t; #define PRIcount PRIu64 struct state { int n; uint64_t p; /* prime not yet examined */ double total_time; size_t state_size; size_t counts_len; count_t counts[]; }; struct state *state; size_t n_to_counts_len (int n) { return (size_binomial(n+9,9)); } void state_alloc (size_t state_size) { state = realloc_or_die (state, state_size); memset (state, '\0', state_size); } void initial_state (int n) { size_t counts_len = n_to_counts_len(n); size_t state_size = size_add_or_die (sizeof(struct state), size_multiply_or_die (counts_len, sizeof(state->counts[0]))); state_alloc (state_size); state->n = n; state->state_size = state_size; state->counts_len = counts_len; } void free_state (void) { free (state); state = NULL; } /* total of all counts, and which is total primes considered */ static uint64_t counts_total (void) { uint64_t ret=0; size_t i; for (i=0; i < state->counts_len; i++) ret += state->counts[i]; return (ret); } /* number of different digs seen, meaning how many count[i] != 0 */ static size_t counts_num_digs_seen (void) { size_t i, ret=0; for (i=0; i < state->counts_len; i++) ret += (state->counts[i] != 0); return (ret); } /* how many times "target" occurs as a count */ static size_t counts_num_equal_to (count_t target) { size_t i, ret=0; for (i=0; i < state->counts_len; i++) ret += (state->counts[i] == target); return (ret); } /* maximum count, which is A065851 */ static count_t counts_max (void) { size_t i, ret=0; for (i=0; i < state->counts_len; i++) if (state->counts[i] > ret) ret = state->counts[i]; return (ret); } /* -------------------------------------------------------------------------- */ void check_A065851_max_count (int n, count_t max_count) { static const count_t want[] = { 0, /* no n=0 */ 1, 2, 4, 11, 39, 148, 731, 4333, 26519, 152526, /* n=10 */ 1251724, /* n=11 */ 7627713, /* n=12 */ 49545041, /* n=13 */ 342729012 /* n=14 */ }; if (! (n>=1)) error("check_A065851_max_count() is for n>=1\n"); if (n >= numberof(want)) { if (option_verbose) printf(" (beyond A065851 data)\n"); return; } if (max_count == want[n]) { if (option_verbose) printf(" (good vs A065851 data)\n"); return; } error("oops, n=%d A065851 count want %"PRIu64" got %"PRIu64"\n", n, want[n], max_count); } void check_A373179_digits (int n, uint64_t digits) { static const uint64_t want[] = { 0, /* no n=0 */ 2, 13, 149, 1237, 13789, 123479, 1235789, 12345679, 102345679, 1123456789, 10123456789, /* n=11 */ 1011233456789, /* n=12 */ 1012334567789, /* n=13 */ 10123345677899 /* n=14 */ }; if (! (n>=1)) error("check_A373179_digits() is for n>=1\n"); if (n >= numberof(want)) { if (option_verbose) printf(" (beyond A373179 data)\n"); return; } if (digits == want[n]) { if (option_verbose) printf(" (good vs A373179 data)\n"); return; } error("oops, n=%d A373179 digits want %"PRIu64" got %"PRIu64"\n", n, want[n], digits); } static const uint64_t A006879_num_primes_table[] = { 0, 4, 21, 143, 1061, 8363, 68906, 586081, 5096876, 45086079, 404204977, /* n=10 */ 3663002302, /* n=11 */ 33489857205, /* n=12 */ 308457624821, /* n=13 */ 2858876213963, /* n=14 2.8 trillion*/ 26639628671867, /* n=15 */ 249393770611256 /* n=16 */ }; void check_A006879_num_primes (int n, uint64_t num_primes) { if (0 <= n && n < numberof(A006879_num_primes_table)) { uint64_t want = A006879_num_primes_table[n]; if (num_primes != want) error("n=%d num_primes want %"PRIu64" got %"PRIu64"\n", n, want, num_primes); } } /* -------------------------------------------------------------------------- */ static const char saved_state_filename[] = "saved-state.data"; static const char progress_filename[] = "progress.txt"; static const char temp_filename[] = "saved-state.txt.tempfile"; static uint64_t save_last_count; void save (int n, uint64_t p) { double last_cputime = cputime_increment(); state->total_time += last_cputime; state->p = p; /* Write first to temp_filename and then rename, so atomic replacement. */ { FILE *fp = fopen_or_die (temp_filename, "wb"); fwrite_or_die (state, state->state_size, 1, fp, temp_filename); ferror_die (fp, temp_filename, "writing"); fclose_or_die (fp, temp_filename); rename_or_die (temp_filename, saved_state_filename); } /* progress */ { FILE *fp = fopen_or_die (temp_filename, "w"); fprintf(fp, "n = %d\n", n); fprintf(fp, "next prime p = %"PRIu64"\n", p); uint64_t num_primes = counts_total(); uint64_t total_num_primes = (1 <= n && n < numberof(A006879_num_primes_table) ? A006879_num_primes_table[n] : 0); fprintf(fp, "number of primes %"PRIu64" of %"PRIu64" is %.2lf%%\n", num_primes, total_num_primes, (total_num_primes ? num_primes * 100.0 / total_num_primes : -1.0)); fprintf(fp, "\n"); fprintf(fp, "total %.1f seconds CPU time\n", state->total_time); fprintf(fp, "last chunk %.1f seconds CPU time\n", last_cputime); double million_per_second = (last_cputime >= 1 ? save_last_count / last_cputime / 1e6 : -1.0); fprintf(fp, " %"PRIu64" many primes is %.2f million/second\n", save_last_count, million_per_second); save_last_count = 0; ferror_die (fp, temp_filename, "writing"); fclose_or_die (fp, temp_filename); rename_or_die (temp_filename, progress_filename); } if (option_verbose >= 2) printf("save to %s n=%d p=%"PRIu64"\n", saved_state_filename, n, p); } static volatile sig_atomic_t flag_interrupt = 0; static volatile sig_atomic_t flag_terminate = 0; void sig_handler_save_state (int signum) { if (option_verbose >= 2) printf("signal %d received for save state\n", signum); flag_interrupt = 1; } void sig_handler_terminate (int signum) { if (option_verbose) printf("signal %d received for terminate\n", signum); flag_interrupt = 1; flag_terminate = 1; } /* Initialise handlers for saving on SIGINT, SIGHUP, SIGTERM, and on a timer. */ void init_save () { struct sigaction s; s.sa_handler = sig_handler_terminate; sigfillset(&s.sa_mask); s.sa_flags = SA_RESTART; sigaction_or_die(SIGINT, &s, NULL); sigaction_or_die(SIGHUP, &s, NULL); sigaction_or_die(SIGTERM, &s, NULL); s.sa_handler = sig_handler_save_state; sigaction_or_die(SIGALRM, &s, NULL); struct itimerval t; t.it_interval.tv_sec = 5*60; /* 5 minutes */ t.it_interval.tv_usec = 0; t.it_value = t.it_interval; setitimer_or_die(ITIMER_REAL, &t, NULL); } /* -------------------------------------------------------------------------- */ /* must have 10*DIGS_BITS_PER_ELEM <= 64 so fit in uint64_t */ #define DIGS_BITS_PER_ELEM 6 #define DIGS_ELEM_MASK UINT64_LOW_MASK( DIGS_BITS_PER_ELEM) #define DIGS_TOTAL_MASK UINT64_LOW_MASK(10 * DIGS_BITS_PER_ELEM) #define DIGS_ABOVE_MASK (~ DIGS_TOTAL_MASK) #define DIGS_ELEM_MAX DIGS_ELEM_MASK /* Return the total number of digits in digs (sum over d=0..2). */ int digs_validate (uint64_t digs) { if (digs & DIGS_ABOVE_MASK) error("oops, digs outside total mask\n"); return (1); } /* Return the total number of digits in digs (sum over d=0..9). */ int digs_num_digits (uint64_t digs) { assert (digs_validate(digs)); int d, ret=0; for (d=0; d<=9; d++, digs >>= DIGS_BITS_PER_ELEM) ret += digs & DIGS_ELEM_MASK; return (ret); } /* In digs, return the count of how many times d occurs. */ static inline int digs_get (uint64_t digs, int d) { assert (digs_validate(digs)); return ((digs >> (d * DIGS_BITS_PER_ELEM)) & DIGS_ELEM_MAX); } void digs_print_counts (uint64_t digs) { int d; for (d=0;d<=9;d++) printf(" %d", digs_get(digs,d)); if (digs & DIGS_ABOVE_MASK) printf(" (oops, extra above too)"); printf("\n"); } void digs_print_digits (uint64_t digs) { int d,i; for (d=0;d<=9;d++) { int num = digs_get(digs,d); for (i=0; i < num; i++) printf("%d",d); } if (digs & DIGS_ABOVE_MASK) printf(" (oops, extra above too)"); printf("\n"); } /* p is a number. Return a "digs" which counts how many times each decimal digit d=0..9 occurs in p. (Ignoring leading 0s.) */ static inline uint64_t p_to_digs (uint64_t p) { uint64_t digs = 0; while (p) { uint64_t q = p / 10; uint64_t r = p % 10; p = q; digs += (uint64_t)1 << (r * DIGS_BITS_PER_ELEM); } DEBUG(printf(" digs "); digs_print_counts(digs)); return (digs); } /* four_to_digs[p] is a digs for the least significant 4 digits of p. p is taken as 4 digits with initial 0s. */ static uint64_t four_to_digs[10000]; static void init_four_to_digs (void) { uint64_t p; for (p=0; p < numberof(four_to_digs); p++) four_to_digs[p] = p_to_digs(p + 10000) - (1 << DIGS_BITS_PER_ELEM); /* sans "1" at 10^4 = 10000 */ } /* p_to_digs_with_cache() takes a number p. Return a "digs" which counts how many times each decimal digit d=0..9 occurs in p. A caching scheme uses a table lookup for low 4 digits, and if the higher digits are the same as last time then re-use the calculation previously made for them. */ static inline uint64_t p_to_digs_with_cache (uint64_t p) { if (UNLIKELY(p < 1000)) return (p_to_digs(p)); assert (p >= 1000); /* 4 or more digits in p */ uint64_t q = p / 10000; uint64_t r = p % 10000; uint64_t digs = four_to_digs[r]; static uint64_t prev_q = 0; static uint64_t prev_digs = 0; if (UNLIKELY(q != prev_q)) { prev_q = q; prev_digs = p_to_digs(q); } return (digs + prev_digs); } static size_t index_table[10][MAXN+1]; static void init_index_table (int n) { if (n < 1) error("Must have n>=1\n"); if (n > MAXN) error("n=%d bigger than MAXN, increase that and re-compile\n", n); int i,d; for (i=1; i<=n; i++) for (d=0; d<=9; d++) index_table[d][i] = size_binomial (i-1 + 9-d, 9-d); DEBUG(printf("index_table[]\n"); for (i=0; i<=n; i++) { printf(" rem=%d: ", i); for (d=0; d<=9; d++) printf(" %2zu%s", index_table[d][i], d==9?"":","); printf("\n"); }); } static inline size_t digs_to_index (int n, uint64_t digs) { /* if the loop drops n to n=0 then the loop could stop, but with GNU C it seems a touch faster unconditional */ assert (digs_num_digits(digs) == n); size_t ret = 0; int d; for (d=0; d<=7; d++) { /* drop to n = how many digits remaining after d */ n -= digs & DIGS_ELEM_MASK; digs >>= DIGS_BITS_PER_ELEM; ret += index_table[d][n]; } /* second final digit d=8 is index_table[d][n] = n, so just add */ assert (d==8); n -= digs & DIGS_ELEM_MASK; digs >>= DIGS_BITS_PER_ELEM; assert (index_table[d][n] == n); ret += index_table[d][n]; /* final d=9 omitted since drop to remaining n=0 only adds 0 */ d++; assert (d == 9); assert (n == digs); /* 9s remaining */ assert (index_table[d][0] == 0); return (ret); } static uint64_t digs_from_index (int n, size_t index) { DEBUG(printf("digs_from_index() index=%zu\n", index)); uint64_t digs = 0; int d; for (d=0; d<=9; d++) { uint64_t num = 0; /* how many of digit d, want as few as possible */ for (;;) { DEBUG(printf(" d=%d num=%"PRIu64" n=%d index=%zu compare %zu\n", d, num, n, index, index_table[d][n])); if (index >= index_table[d][n]) break; num++; n--; if (n==0) break; } assert (index >= index_table[d][n]); index -= index_table[d][n]; assert (num <= DIGS_ELEM_MAX); digs += num << (d * DIGS_BITS_PER_ELEM); DEBUG(printf(" set d=%d num=%"PRIu64" to n=%d skip %zu to index=%zu\n", d, num, n, index_table[d][n], index)); } DEBUG(printf(" final index=%zu digs %"PRIX64"\n", index, digs)); return (digs); } /* digs[d] is count of how many of digit d. Return the smallest number which has these digs digits, with repetition. The return starts with the smallest non-0 digit, then the rest in ascending order. */ static inline size_t digs_to_smallest (int n, uint64_t digs) { /* find smallest d>=1 digit in the set */ int d=1; while (UNLIKELY (digs_get(digs,d) == 0)) d++; /* digit d then 0s */ uint64_t ret = d * uint64_pow10_table[digs_get(digs,0)]; /* other repetitions of d, after its one initial */ int rep; for (rep = digs_get(digs,d) - 1; rep > 0; rep--) ret = 10*ret + d; /* digits after d */ for (d++; d<=9; d++) for (rep=digs_get(digs,d); rep > 0; rep--) ret = 10*ret + d; return (ret); } void show_configuration (void) { static int done = 0; if (done) return; done = 1; if (option_verbose) { printf("primesieve version %s\n", primesieve_version()); printf("MAXN = %d\n", MAXN); printf("four_to_digs[] table size %zu bytes\n", sizeof(four_to_digs)); printf("assert()s %s\n", WANT_ASSERT ? "enabled" : "disabled"); } else { assert (printf("assert()s enabled\n")); } } /* Among all counts[index] == target, Return the smallest digs_to_smallest() of those indices. */ static uint64_t counts_smallest_at (int n, size_t target) { size_t i; uint64_t ret=UINT64_MAX; for (i=0; i < state->counts_len; i++) { if (state->counts[i] == target) { uint64_t digs = digs_from_index(n,i); assert (digs_num_digits(digs) == n); uint64_t s = digs_to_smallest(n,digs); if (s < ret) ret = s; } } return (ret); } void output (int n) { state->total_time += cputime_increment(); uint64_t num_primes = counts_total(); count_t max_count = counts_max(); uint64_t max_count_at = counts_smallest_at (n,max_count); printf("n=%d max count %"PRIcount" at digits %"PRIu64"\n", n, max_count, max_count_at); if (option_verbose) { size_t num_ways = counts_num_equal_to(max_count); printf(" %zu equal max count, digit combination(s):\n" , num_ways); size_t i; for (i=0; i < state->counts_len; i++) { if (state->counts[i] == max_count) { uint64_t digs = digs_from_index(n,i); printf(" "); digs_print_digits(digs); } } printf(" %"PRIu64" primes with %zu digit combs\n", num_primes, counts_num_digs_seen()); printf(" counts table %zu entries, %.1lf mbytes\n", state->counts_len, state->counts_len * sizeof(state->counts[0]) / 1e6); printf(" time %.1f seconds is %.1lf million primes/second\n", state->total_time, state->total_time==0 ? -1 : num_primes / state->total_time / 1e6); } check_A065851_max_count (n, max_count); check_A373179_digits (n, max_count_at); check_A006879_num_primes (n, counts_total()); if (option_verbose) printf("\n"); } int search_from_state (void) { int n = state->n; init_index_table(n); cputime_increment(); /* dummy initial */ /* state->p is a prime which has not yet been examined, jump the iterator to that p */ uint64_t prime_range_hi = uint64_pow10(n) - 1; primesieve_iterator it; primesieve_init (&it); primesieve_jump_to(&it, state->p, prime_range_hi); for (;;) { uint64_t p = primesieve_next_prime (&it); if (UNLIKELY (p > prime_range_hi)) break; DEBUG(printf("prime p=%"PRIu64" (%"PRIu64" many previous)\n", p, counts_total())); if (UNLIKELY(flag_interrupt)) { flag_interrupt = 0; save (n, p); if (flag_terminate) return(0); } uint64_t digs = p_to_digs_with_cache (p); save_last_count++; size_t index = digs_to_index(n,digs); assert (0 <= index && index < state->counts_len); count_t count = ++ state->counts[index]; if (UNLIKELY(count==0)) error("count overflow\n"); } output (n); free_state (); primesieve_free_iterator (&it); return (1); } int show (int n) { show_configuration(); if (! (n>=1)) error("calculation is for 1 <= n <= %d>=1, got %d\n", n, MAXN); initial_state(n); state->p = uint64_pow10 (n-1); return (search_from_state ()); } void resume (void) { FILE *fp = fopen_or_die (saved_state_filename, "rb"); size_t state_size = fsize_or_die (fp, saved_state_filename); state_alloc (state_size); fread_or_die (state, state_size, 1, fp, saved_state_filename); ferror_die (fp, saved_state_filename, "reading"); fclose_or_die (fp, saved_state_filename); if (state->state_size != state_size || state->counts_len != n_to_counts_len(state->n)) error("oops, wrong lengths in %s\n", saved_state_filename); if (option_verbose >= 2) printf("resume n=%d at p=%"PRIu64"\n", state->n, state->p); search_from_state (); } int main (int argc, char **argv) { setbuf(stdout,NULL); init_save(); init_four_to_digs(); DEBUG (option_verbose = 1); int option_resume = 0; int i, endpos, seen_n=0; int n, nhi; for (i = 1; i < argc; i++) { const char *arg = argv[i]; if (strcmp(arg,"-v") == 0) { option_verbose++; } else if (strcmp(arg,"resume") == 0) { option_resume = 1; } else if (sscanf(arg, "%d%n", &n, &endpos) == 1 && endpos == strlen(arg)) { seen_n = 1; show (n); } else if (sscanf(arg, "%d-%d%n", &n, &nhi, &endpos) == 2 && endpos == strlen(arg)) { seen_n = 1; for ( ; n<=nhi; n++) if (! show(n)) break; } else { error("unrecognised command line option: %s\n", arg); } } if (!seen_n) { if (option_resume) resume(); else for (n=1; n<=6; n++) show (n); } return (0); }