/* Kevin Ryde, November 2024 Usage: ./a.out [-v] This is some C code for GNU C and a 64 bit build, using its __int128 feature, to find terms of A261870(n) = minimum number of steps required to go from n to 1 under x -> 3x-1 if x odd, x -> 3x-1 or x/2 if x even. A378375(n) = number of different ways for that minimum number of steps A378376(n) = smallest starting x which requires n steps The method is an explore back through predecessors of 1, building a best[x] array of minimum steps and number of ways. The explore is pruned at x already seen in the best[] array at fewer steps, and at x too big to drop down to affect the best[] array within MAX_STEPS. Configuration ------------- The #defines below are compile-time constants for MAX_STEPS explore <= MAX_STEPS MAX_WATCH best[x] holds results for x <= MAX_WATCH WANT_WAYS whether to record how many ways The supplied MAX_STEPS 71 and MAX_WATCH 1000 runs quickly and is enough for 1000 terms of A261870 and A378375. Try MAX_STEPS 98 and MAX_WATCH 10000 for 10000 terms in about 6 minutes on a 3 GHz CPU at the time of writing. These are compile-time constants with the aim of helping code generation as much as possible when attempting large MAX_STEPS. MAX_STEPS determines time taken. On average there's 1+1/3 predecessors of each x so work proportional to roughly (4/3)^MAX_STEPS, though pruning reduces it. MAX_WATCH is the largest n for A261870 steps and A378375 ways, and largest term for A378376 first of n steps. Full results depend on MAX_STEPS big enough to find all those. Larger MAX_WATCH can prune more: see "MAX_WATCH Size" below. WANT_WAYS selects whether to record the A378375 "ways". Doing so re-explores any equal-best x so as to count again. The extra work is maybe 2* or more. Set WANT_WAYS 0 for more speed if seeking just min steps. Output ------ Output is 3 b-files at the end of all exploring. Each stop at their first unknown term (within the configured options). b261870.txt minimum number of steps b378375.txt number of ways for that minimum (if WANT_WAYS) b378376.txt smallest x requiring n steps Optional command line parameter -v prints more verbose output. Doubled -v -v shows best[] progress, within x <= MAX_WATCH. This progress might be too much output if MAX_WATCH is large, but of interest for small MAX_STEPS and/or MAX_WATCH. Save and Resume --------------- To help long runs which might have to be restarted, the state of the search is saved at 15 minute intervals and on SIGINT, SIGTERM, SIGHUP termination. saved-state.data binary data file progress.txt summary text file The run can be resumed from saved-state.data by ./a.out resume saved-state.data is a dump of the "struct state". Its contents depend on the compile-time configuration options, and the machine data type sizes, alignments, endianness, etc. progress.txt includes the stack[] array of pending work. Pushes put s values +2 apart. Places more than +2 are where searching has been down to take an entry and put a new run of +2 starting there, and so shows work advancing. Possibly Never Reaching 1 ------------------------- It's unknown, presently, whether every starting x eventually reach 1. But no such x is known in the feasible range of calculation for A261870 steps and A378375 ways here. The Collatz 3x-1 problem is x -> x/2 when x even, without a choice when x even. Explicit trajectory testing shows every x there reaches 1 or 5 or 17 for all x <= 2*10^12. (The choice x -> 3x-1 when x even allows 5 and 17 to go to 1.) If some x never reaches 1, then the output here for A261870 steps and A378375 ways would always stop at that x, no matter how large MAX_STEPS. But A378376 smallest of s steps is unaffected: it only ever considers those x which do reach 1. Predecessors ------------ A given x has one or two predecessors y which step to x, y = 2*x step y -> y/2 = x y = (x+1)/3 step y -> 3*x-1 = x if (x+1)/3 is an integer, which is x==2 mod 3 search() starts at x=1 and works backward through predecessors of 1, predecessors of predecessors of 1, and so on, to find all x which reach 1 after "s" steps, for s <= MAX_STEPS.. 1 s=0 | predecessors 2 s=1 tree | \ 4 1(again) s=2 | 8 ______________ s=3 | \ 16 3 s=4 | | 32 6 s=5 | \ | 64 11 12 s=6 | |\ | 128 22 4(again) 24 s=7 search() first uses fill_doublings() to apply all "2*x" doubling predecessors of x into the best[] results. x at s steps 2*x at s+1 steps 4*x at s+2 steps ... x*2^e at s+e steps with x*2^e <= MAX_WATCH search() then pushes onto the stack[] array (x+1)/3 predecessor of each of those doublings, when an integer. push (x*2^e + 1)/3 for e >= 0, at s+1+e <= MAX_STEPS, when an integer If x == 0 mod 3 then (x*2^e + 1)/3 is never an integer and nothing to push. The "3" chain in the diagram above is this type. All its doublings are considered, but they don't branch off to additional predecessors. If x == 2 mod 3 then (x*2^e + 1)/3 is an integer iff e is even. So push ( x + 1)/3 at s+1 steps ( 4*x + 1)/3 at s+3 steps (16*x + 1)/3 at s+5 steps ... (x*4^f + 1)/3 at s+1 + 2*f steps for f>=0 and s+1 + 2*f <= MAX_STEPS and can_reach_watch() I x == 1 mod 3 then (x*2^e + 1)/3 is an integer iff e is odd. This is handled by stepping x -> 2*x at s+1 steps, then the same as x==2 mod 3. x at s steps taken from the stack[] pushes at most new entries s+1, s+3, etc. So the stack never holds more than MAX_STEPS entries (or rather half that in fact). num_x count which is shown in progress.txt and in final results counts each x considered by search() for filling and pushing in the manner above. This means each x which goes onto the stack[]. Filling Extent -------------- fill_doublings() updates best[] for all x*2^e <= MAX_WATCH, even if that makes some best[].min_steps > MAX_STEPS. The idea is not to bother with an extra s <= MAX_STEPS. That condition is unlikely to shorten the loop and anything bigger than MAX_STEPS is ignored in final results. Prune by best[].min_steps ------------------------- If some x <= MAX_WATCH is reached at s steps then best[x].min_steps < s means x was seen before at fewer steps, so no improvement, prune the search Searching through predecessors of this x again would only find similarly bigger s than seen before. The simplest case is x=1 reached a second time at s=2 by 1 <- 2 <- 1 predecessor 2*x then (x+1)/3 But 1 is the start and was seen before at 0 steps. Searching its predecessors again would be exactly the same as already in progress, but at bigger s. The tree diagram shown above has these kind of repeats marked "1(again)" and "4(again)". stack[] still holds these repeats (even the initial x=1,s=2). They're checked for "no improvement" when popped. best[].min_steps might have decreased at that time, so more opportunity to prune than looking when pushed. An entry on the stack for more time costs nothing. When WANT_WAYS = 0, the prune is strengthened to "<=" best[x].min_steps <= s means x was seen before at same or fewer steps, so no improvement, prune the search Searching predecessors of this x again would only find the same (or worse) than already known. There can be multiple equivalent paths from n to 1. WANT_WAYS wants to know how many, but min_steps and smallest x of s steps don't. Prune by can_reach_watch() -------------------------- When s is close to MAX_STEPS, it has only a few more potential (x+1)/3 predecessor steps. If x is big then those steps cannot reach <= MAX_WATCH so cannot affect the results and this x can be pruned. can_reach_watch(x,s) returns true if x at s steps is small enough that more steps, through to MAX_STEPS, might potentially reach <= MAX_WATCH. can_reach_watch_table[] is set by thinking in the other direction: 3x-1 successors up from MAX_WATCH, more = MAX_STEPS - s can_reach_watch_table[MAX_STEPS ] = MAX_WATCH can_reach_watch_table[MAX_STEPS - 1] = 3*MAX_WATCH - 1 can_reach_watch_table[MAX_STEPS - 2] = 3*(3*MAX_WATCH - 1) - 1 These are how big MAX_WATCH becomes under "more" many 3x-1 steps. Table entries are capped at VALUE_MAX when 3x-1 overflows value_t. That happens for small s (large "more") and x > VALUE_MAX never happens so the effect is no pruning at small s. can_reach_watch() greatly reduces work done at large x: large and irrelevant that is, for finding results x <= MAX_WATCH. MAX_WATCH Size -------------- MAX_WATCH is always set to the number of terms of A261870 steps and A378375 ways desired, and to the largest term of A378376 smallest of s steps desired. The latter is hard to know in advance and some attempts or likely guess might be necessary. A yet bigger MAX_WATCH can be used for more opportunity to prune as described above. How much bigger is a balance between the work of building and maintaining best[] entries (fill_doublings()), against how much searching is pruned. Very big MAX_WATCH does a lot of work for not much more pruning. CPU memory caching is probably a factor too (on today's typical multi-level RAM), due to slower memory access once outside caches. Might have to experiment on a given machine. Division by 3 ------------- Division by 3, with remainder, for (x+1)/3 predecessors, is handled by an explicit multiply by modular inverse. GCC knows modular inverses for testing remainders, but its version 13 didn't seem to come out as well it could when combining that with moving to x+1 or 2x+1 to make remainder 0, ready for (x+1)/3 or (2x+1)/3. The modular inverse method is M = 2^VALUE_BITS = within the value_t type q = (x * VALUE_DIV3_MUL) mod M x==0 mod 3 iff q <= VALUE_DIV3_CMP_0 x==1 mod 3 iff q >= VALUE_DIV3_CMP_1 x==2 mod 3 otherwise where constants VALUE_DIV3_MUL = 3^-1 mod M = 0xAA...AAB which is (3 * VALUE_DIV3_MUL) == 1 mod M VALUE_DIV3_CMP_0 = floor(VALUE_MAX / 3) VALUE_DIV3_CMP_1 = VALUE_DIV3_MUL These cases and constants are for even VALUE_BITS. Equivalents are easy for odd bits too, but a machine with odd number of bits in its word size is almost unheard of today. If x==1 mod 3 then x going to 2*x, to make 2 mod 3, is as easy as q to 2*q (mod M), as if (2*x)*VALUE_DIV3_MUL (mod M). If x==2 mod 3 then y = (x+1)/3 is an integer = q + VALUE_DIV3_MUL (mod M) as if from (x+1)*VALUE_DIV3_MUL Value Type Size --------------- The maximum value reached in a search is all 2*x predecessors, x <= 2^MAX_STEPS value_t type must be big enough for that. check_configuration() checks this at the start and the rest of the code runs freely knowing overflow will not occur. value_t is set to GNU C extension "unsigned __int128". This is a convenient way to act on values near 100 bits when seeking terms to n=10000 in MAX_STEPS 98 or thereabouts. Change value_t to some other type if desired, as long as 2^MAX_STEPS fits the chosen MAX_STEPS. Ways Counting ------------- When an equal smallest s == best[x].min_steps is found, this is an additional way and so increment best[x].ways. The search continues back from x to increment corresponding ways for further predecessors too (those for which this is still an equal smallest way). There's no attempt to record the predecessors which will be seen by this re-explore. Only those <= MAX_WATCH are of interest, but the worry would be how much memory that requires for each x <= MAX_WATCH and whether a net saving over re-explore. Some best[].ways during the search can be much larger than the final result. Overflow is detected when incrementing. Increase the field or type in "struct best" if necessary. A378376 formulas set out some relations for number of ways based on where n will or can step to, a(n) = a(3n-1) if n odd a(n) = / a(3n-1) \ | or a(x/2) | if n even \ or a(3n-1) + a(x/2) / These are not useful for calculating terms, as such. a(3n-1) is always required for n odd, and might be required for n even. Only an explore for the minimum steps for each n will say which. The ways and min_steps can easily be due to trajectories with terms > MAX_WATCH. Those terms are explored in full here, but the results recorded are limited to MAX_WATCH. Things Not Done --------------- x values become sparse as they increase so pruning might hold some of the bigger previously seen x in a tree or hash. But the expectation would be that costs of lookups and maintaining may easily exceed repeat explores saved. Popping the largest s from the stack[] each time tends to start best[] with some relatively large (poor) min_steps. But taking smaller s from the stack increases the space needed there. The worst case would be taking the very smallest pending s every time, which becomes a breadth-first search back and the number of pending entries would grow somewhere like (4/3)^s. */ /* ---------------------------------------------------------------------- */ /* Configuration */ #define MAX_STEPS 71 #define MAX_WATCH 1000 #define WANT_WAYS 1 /* ---------------------------------------------------------------------- */ /* Set WANT_ASSERT to 1 for various assert() checks (much slower). Set WANT_DEBUG to 1 for verbose diagnostic prints. */ #define WANT_ASSERT 0 #define WANT_DEBUG 0 #if ! WANT_ASSERT #define NDEBUG #endif #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 MAX(x,y) ((x)>=(y) ? (x) : (y)) #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); } 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 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 short\n", filename); } static void fexpect_eof_or_die (FILE *fp, const char *filename) { int c = fgetc(fp); ferror_die(fp, filename, "testing EOF"); if (c != EOF) error("File too long: %s\n", filename); } 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)); } static void clock_gettime_or_die (clockid_t id, struct timespec *t, const char *name) { if (clock_gettime(id, t) != 0) error("Cannot clock_gettime() %s: %s\n", name, strerror(errno)); } static double clock_gettime_double_sec_or_die (clockid_t id, const char *name) { struct timespec t; clock_gettime_or_die (id, &t, name); return ((double) t.tv_sec + t.tv_nsec/1e9); } /* elapsed_cputime() returns the amount of CPU time in seconds used since the last call to elapsed_cputime(), or since program start for the first call. */ double elapsed_cputime (void) { static double last_time = 0.0; double t = clock_gettime_double_sec_or_die (CLOCK_PROCESS_CPUTIME_ID, "CPUTIME"); double elapsed = t - last_time; last_time = t; return (elapsed); } /* -------------------------------------------------------------------------- */ static int option_verbose = 0; void check_A261870_num_steps (int n, int s) { /* OEIS_data("A261870",'bfile) b261870.txt */ static int want[] = { -1, /* no n=0 */ 0, 1, 4, 2, 15, 5, 13, 3, 16, 11, 6, 6, 19, 14, 9, 4, 35, 17, 17, 12, 25, 7, 20, 7, 33, 15, 15, 15, 10, 10, 23, 5, 36, 36, 18, 18, 31, 18, 13, 13, 26, 26, 8, 8, 26, 21, 21, 8, 34, 34, 21, 16, 16, 16, 29, 16, 29, 11, 11, 11, 24, 24, 24, 6, 37, 37, 24, 19, 19, 19, 19, 19, 32, 32, 19, 19, 14, 14, 14, 14, 27, 27, 27, 27, 27, 9, 40, 9, 27, 27, 22, 22, 22, 22, 22, 9, 35, 35, 35, 22, 35, 22, 17, 17, 17, 17, 17, 17, 30, 30, 30, 17, 30, 30, 12, 12, 43, 12, 30, 12, 25, 25, 25, 25, 25, 25, 25, 7, 38, 38, 38, 38, 38, 25, 38, 20, 20, 20, 20, 20, 20, 20, 20, 20, 33, 33, 33, 33, 33, 20, 33, 20, 46, 15, 15, 15, 15, 15, 33, 15, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 10, 10, 41, 41, 41, 10, 41, 28, 28, 28, 23, 23, 23, 23, 54, 23, 23, 23, 23, 23, 23, 10, 36, 36, 36, 36, 36, 36, 36, 23, 36, 36, 18, 23, 18, 18, 18, 18, 49, 18, 18, 18, 36, 18, 31, 18, 31, 31, 31, 31, 31, 31, 31, 18, 44, 31, 31, 13, 44, 13, 44, 13, 44, 44, 13, 13, 31, 31, 31, 13, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 8, 39, 39, 39, 39, 39, 39, 39, 39, 39, 26, 26, 26, 39, 39, 21, 21, 52, 21, 21, 21, 52, 21, 21, 21, 39, 21, 21, 21, 21, 21, 34, 21, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 21, 21, 34, 34, 34, 16, 47, 47, 16, 16, 47, 16, 47, 16, 47, 16, 16, 16, 34, 34, 34, 16, 60, 29, 29, 29, 29, 29, 29, 29, 60, 29, 29, 29, 29, 29, 29, 29, 47, 29, 29, 29, 42, 11, 42, 11, 42, 42, 42, 42, 42, 42, 42, 11, 42, 42, 29, 29, 42, 29, 42, 29, 24, 24, 24, 24, 24, 24, 24, 24, 55, 55, 24, 24, 24, 24, 24, 24, 42, 24, 24, 24, 24, 24, 37, 11, 68, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 24, 24, 37, 37, 37, 37, 50, 19, 50, 24, 50, 19, 19, 19, 19, 19, 50, 19, 50, 50, 19, 19, 37, 19, 37, 19, 37, 37, 19, 19, 32, 32, 32, 19, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 19, 50, 32, 32, 32, 32, 32, 45, 14, 45, 45, 14, 14, 45, 45, 45, 14, 45, 45, 45, 45, 45, 14, 45, 14, 32, 32, 32, 32, 32, 32, 45, 14, 58, 27, 27, 27, 58, 27, 27, 27, 27, 27, 27, 27, 58, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 45, 27, 27, 27, 27, 27, 40, 9, 71, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 27, 40, 27, 40, 27, 40, 40, 40, 22, 22, 22, 53, 22, 53, 53, 22, 22, 53, 22, 22, 22, 53, 53, 22, 22, 22, 22, 22, 22, 40, 40, 22, 22, 40, 22, 40, 22, 35, 22, 22, 22, 35, 35, 35, 22, 66, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 22, 35, 22, 35, 35, 35, 35, 35, 35, 48, 17, 48, 48, 48, 17, 48, 17, 48, 17, 48, 48, 17, 17, 48, 48, 48, 17, 48, 48, 17, 17, 48, 17, 35, 17, 35, 35, 35, 35, 35, 35, 48, 17, 61, 61, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 61, 30, 30, 30, 30, 30, 30, 30, 48, 30, 30, 30, 30, 30, 30, 30, 48, 43, 30, 30, 30, 30, 30, 30, 43, 43, 12, 12, 43, 43, 43, 12, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 12, 43, 43, 43, 43, 30, 30, 43, 30, 43, 43, 30, 30, 43, 43, 43, 25, 56, 25, 25, 25, 56, 25, 56, 25, 56, 25, 25, 25, 25, 25, 25, 25, 56, 56, 56, 25, 56, 25, 25, 25, 25, 25, 25, 25, 43, 25, 25, 25, 43, 43, 25, 25, 43, 25, 25, 25, 38, 25, 25, 25, 38, 38, 38, 12, 69, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 25, 25, 38, 25, 38, 38, 38, 38, 38, 38, 38, 38, 51, 20, 20, 20, 51, 51, 51, 20, 51, 51, 20, 20, 51, 20, 51, 20, 51, 20, 20, 20, 51, 51, 51, 20, 51, 51, 20, 20, 51, 20, 51, 20, 38, 38, 20, 20, 38, 38, 38, 20, 51, 38, 38, 38, 33, 20, 33, 20, 33, 33, 33, 33, 33, 33, 33, 20, 64, 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, 64, 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, 20, 64, 51, 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, 46, 46, 15, 15, 46, 46, 46, 46, 46, 15, 46, 15, 46, 46, 46, 46, 46, 46, 46, 15, 46, 46, 46, 46, 46, 46, 46, 46, 46, 46, 15, 15, 46, 46, 46, 15, 33, 33, 33, 33, 33, 33, 46, 33, 46, 33, 33, 33, 46, 46, 46, 15, 59, 59, 28, 28, 59, 28, 28, 28, 59, 28, 28, 28, 28, 28, 28, 28, 59, 28, 28, 28, 28, 28, 59, 28, 59, 59, 28, 28, 28, 28, 28, 28, 59, 28, 28, 28, 46, 28, 46, 28 }; const char *func = "check_A261870_num_steps"; const char *A_number = "A261870"; int min_n = 1; if (! (n >= min_n)) error("oops, %s() is for n>=%d\n", func, min_n); if (n >= numberof(want)) { if (option_verbose) printf("%s(%d) beyond sample data\n", A_number,n); return; } if (s == want[n]) { if (option_verbose) printf("%s(%d) good vs sample data\n", A_number,n); return; } error("oops %s(%d) want x=%d got %d\n", A_number, n, want[n], s); } void check_A378375_ways (int n, uint64_t ways) { static int want[] = { -1, /* no n=0 */ 1, 1, 1, 1, 2, 1, 2, 1, 1, 1, 1, 1, 2, 2, 1, 1, 22, 2, 2, 2, 3, 1, 1, 1, 15, 1, 2, 2, 1, 1, 3, 1, 15, 34, 2, 4, 8, 2, 1, 2, 2, 7, 1, 1, 2, 2, 3, 1, 11, 22, 2, 2, 1, 2, 8, 2, 5, 1, 1, 1, 1, 3, 4, 1, 10, 30, 2, 1, 1, 2, 1, 4, 6, 15, 2, 2, 1, 1, 1, 2, 1, 3, 6, 11, 3, 1, 34, 1, 2, 2, 1, 3, 2, 3, 1, 1, 5, 15, 33, 2, 12, 2, 1, 2, 1, 1, 2, 2, 4, 8, 9, 2, 3, 5, 1, 1, 70, 1, 4, 1, 1, 2, 5, 5, 4, 4, 3, 1, 5, 15, 41, 67, 16, 2, 26, 1, 1, 2, 2, 3, 1, 1, 2, 4, 4, 11, 12, 29, 7, 2, 10, 2, 43, 1, 1, 1, 1, 1, 4, 2, 1, 1, 4, 8, 2, 7, 5, 12, 1, 5, 1, 1, 29, 40, 57, 1, 21, 2, 2, 2, 1, 1, 2, 4, 459, 3, 4, 4, 1, 1, 2, 1, 5, 10, 19, 33, 15, 34, 15, 2, 8, 19, 1, 2, 1, 1, 1, 2, 84, 1, 1, 1, 2, 2, 2, 2, 4, 6, 8, 10, 7, 9, 7, 2, 8, 5, 7, 1, 25, 1, 55, 1, 60, 117, 1, 1, 4, 4, 2, 1, 1, 1, 2, 4, 1, 5, 2, 6, 1, 4, 4, 5, 1, 3, 2, 1, 4, 11, 29, 34, 17, 43, 43, 75, 15, 2, 2, 2, 16, 26, 1, 1, 149, 1, 1, 3, 272, 2, 1, 3, 12, 1, 1, 1, 2, 2, 5, 4, 3, 5, 13, 19, 4, 15, 11, 33, 6, 11, 2, 2, 6, 12, 7, 1, 29, 77, 1, 1, 54, 1, 160, 1, 73, 1, 1, 1, 4, 4, 2, 2, 254, 1, 3, 4, 2, 4, 3, 8, 1129, 4, 7, 9, 4, 6, 4, 13, 2, 3, 5, 7, 13, 1, 21, 1, 33, 34, 46, 74, 49, 70, 55, 1, 8, 25, 4, 4, 19, 2, 30, 2, 1, 1, 2, 2, 1, 2, 1, 5, 345, 620, 2, 4, 1, 4, 1, 4, 6, 1, 1, 3, 2, 2, 8, 1, 1526, 5, 14, 24, 5, 19, 10, 41, 10, 29, 37, 42, 11, 16, 2, 2, 5, 13, 16, 26, 64, 1, 102, 2, 124, 1, 1, 1, 1, 1, 244, 2, 74, 104, 1, 2, 8, 1, 4, 1, 2, 2, 2, 2, 2, 2, 3, 2, 1, 4, 5, 11, 5, 8, 4, 12, 4, 9, 14, 13, 5, 7, 4, 2, 2, 2, 6, 10, 2, 7, 37, 1, 18, 43, 1, 1, 34, 74, 50, 1, 22, 76, 112, 180, 57, 1, 31, 1, 4, 4, 4, 4, 2, 2, 44, 1, 112, 1, 1, 3, 520, 2, 1, 4, 1, 2, 5, 5, 809, 2, 4, 7, 1, 2, 4, 5, 1, 4, 1, 5, 2, 1, 3, 5, 2, 2, 11, 1, 2352, 10, 15, 28, 19, 29, 21, 40, 6, 36, 46, 57, 24, 47, 39, 84, 8, 21, 20, 4, 1, 2, 15, 2, 6, 23, 26, 1, 1, 1, 152, 1, 99, 169, 1, 2, 118, 1, 1, 3, 254, 459, 2, 2, 1, 1, 1, 4, 10, 12, 1, 1, 4, 1, 4, 1, 4, 2, 2, 2, 2, 5, 4, 4, 676, 5, 8, 12, 5, 13, 5, 19, 2, 6, 18, 20, 8, 15, 14, 37, 1, 7, 11, 15, 2, 2, 1, 2, 5, 8, 10, 14, 7, 7, 61, 1, 26, 45, 91, 1, 62, 1, 60, 1, 52, 108, 1, 1, 82, 181, 140, 1, 58, 84, 1, 1, 48, 1, 8, 1, 2, 4, 4, 6, 2, 2, 70, 2, 209, 520, 2, 2, 2, 3, 2, 4, 1, 4, 4, 5, 2, 4, 5, 8, 699, 2, 2, 9, 3, 7, 2, 10, 32, 4, 5, 7, 4, 4, 1, 13, 2, 8, 3, 6, 2, 5, 2, 7, 8, 20, 1, 1, 12, 25, 18, 1, 20, 38, 34, 55, 19, 46, 24, 94, 15, 60, 73, 98, 47, 58, 50, 1, 8, 12, 39, 46, 2, 4, 17, 4, 11, 33, 2, 2, 16, 37, 26, 1, 71, 1, 1, 1, 289, 2, 417, 2, 165, 1, 2, 2, 1, 1, 1, 5, 179, 499, 785, 2, 202, 2, 1, 4, 1, 1, 4, 4, 16, 1, 1, 4, 4, 8, 1, 1, 2, 1, 2, 3, 3, 2, 2, 2, 9, 8, 7, 1, 655, 4, 11, 13, 6, 14, 17, 29, 5, 13, 19, 23, 4, 17, 21, 46, 2, 12, 20, 43, 4, 37, 8, 43, 8, 15, 18, 19, 2, 2, 1, 2, 2, 5, 14, 25, 4, 16, 7, 26, 67, 1, 1, 1, 58, 121, 156, 1, 71, 149, 1, 1, 74, 1, 112, 1, 118, 1, 1, 1, 202, 272, 211, 2, 19, 106, 1, 1, 58, 1, 77, 2, 6, 12, 1, 1, 4, 4, 4, 1, 42, 2, 2, 2, 2, 2, 1, 2, 1, 2, 2, 5, 3, 3, 3, 2, 243, 3, 6, 8, 1, 6, 4, 13, 5, 5, 8, 8, 2054, 4, 2, 17, 3, 4, 8, 11, 3, 14, 4, 13, 1, 6, 7, 11, 4, 5, 2, 2, 114, 2, 1, 5, 2, 6, 6, 10, 2, 4, 7, 7, 24, 46, 1, 1, 20, 29, 41, 65, 34, 1, 48, 1, 10, 39, 55, 105, 52, 54, 66, 1, 15, 28, 79, 160, 51, 126, 59, 230, 50, 73, 1, 1, 23, 41, 45, 1, 2, 4, 4, 4, 2, 4, 38, 4, 19, 2, 2, 2, 30, 44, 26, 1, 103, 254, 2, 2, 273, 1, 1, 3, 380, 1, 3, 2, 1, 2, 2, 4, 483, 1, 2, 3, 1, 5, 831, 5, 580, 1129, 2, 2, 2, 4, 3, 7, 100, 1, 2, 6, 14, 4, 16, 5 }; const char *func = "check_A378375_ways"; const char *A_number = "A378375"; int min_n = 1; if (! (n >= min_n)) error("oops, %s() is for n>=%d\n", func, min_n); if (n >= numberof(want)) { if (option_verbose) printf("%s(%d) beyond sample data\n", A_number,n); return; } if (ways == want[n]) { if (option_verbose) printf("%s(%d) good vs sample data\n", A_number,n); return; } error("oops %s(%d) want x=%d got %zu\n", A_number, n, want[n], ways); } void check_A378376_smallest_of_s (int n, size_t x) { /* OEIS_data("A378376",'bfile) b261870.txt */ static int want[] = { 1, 2, 4, 8, 3, 6, 11, 22, 43, 15, 29, 10, 20, 7, 14, 5, 9, 18, 35, 13, 23, 46, 91, 31, 61, 21, 41, 81, 161, 55, 109, 37, 73, 25, 49, 17, 33, 65, 129, 257, 87, 173, 341, 117, 225, 455, 153, 305, 607, 209, 405, 809, 273, 543, 185, 369, 721, 1433, 481, 961, 321, 641, 1281, 2561, 865, 1713, 577, 1153, 385, 769, 1537, 513, 1041, 2049, 4049, 8049, 2769, 5397, 10773, 3609, 7217, 2457, 4817, 9623, 3281, 6497, 2193, 4385, 8529, 16961, 5825, 11417, 22743, 7713, 15201, 30401, 10353, 20297, 6929, 13857, 27033, 53889, 18401, 35937, 71873, 24545, 48513, 96081, 32769, 64145, 127761, 43713, 85345, 170369, 57217, 113793, 38785, 76289, 25857, 51713, 103809 }; const char *func = "check_A378376_smallest_of_s"; const char *A_number = "A378376"; int min_n = 0; if (! (n >= 0)) error("oops, %s() is for n>=%d\n", func, min_n); if (n >= numberof(want)) { if (option_verbose) printf("%s(%d) beyond sample data\n", A_number,n); return; } if (x == want[n]) { if (option_verbose) printf("%s(%d) good vs sample data\n", A_number,n); return; } error("oops %s(%d) want x=%d got %zu\n", A_number, n, want[n], x); } /* -------------------------------------------------------------------------- */ /* value_t */ /* GNU C type __int128 available on a 64 bit target. Change value_t to something bigger or smaller if desired. */ typedef unsigned __int128 uint128_t; typedef uint128_t value_t; #define VALUE_BITS ((int) (sizeof(value_t) * CHAR_BIT)) #define VALUE_MAX ((value_t) (~ (value_t) 0)) /* VALUE_DIV3_MUL = multiplicative inverse of 3 mod value_t size, so 3*VALUE_DIV3_MUL == 1. Then with q = x*VALUE_DIV3_MUL, q <= VALUE_DIV3_CMP_0 iff x == 0 mod 3 q >= VALUE_DIV3_CMP_1 iff x == 1 mod 3 This requires VALUE_BITS even (see check_configuration()). */ #define VALUE_DIV3_MUL (2*(VALUE_MAX/3) + 1) #define VALUE_DIV3_CMP_0 (VALUE_MAX/3) #define VALUE_DIV3_CMP_1 VALUE_DIV3_MUL /* Return a string which is the decimal digits of x. Explicit conversion since __int128 not in printf() etc. */ char * value_str (value_t x) { /* 3 decimals per byte is enough, plus \0 terminator */ static char buf[3*sizeof(value_t) + 1]; int i = sizeof(buf)-1; do { if (--i < 0) error("value_str() buffer too small\n"); int digit = x % 10; x /= 10; buf[i] = digit + '0'; } while (x); return (buf + i); } /* Return the (x+1)/3 predecessor of x, or error() if not an integer. */ static inline value_t predecessor_third (value_t x) { if (x%3 != 2) error("predecessor_third() only exists for x==2 mod 3\n"); return ((x-2)/3 + 1); } /* (VALUE_MAX + 1) / 3, largest x for which successor 3x-1 <= VALUE_MAX */ #define SUCCESSOR_THIRD_MAX \ ((VALUE_MAX / 3) + (VALUE_MAX % 3 == 2)) static inline value_t successor_third (value_t x) { if (x > SUCCESSOR_THIRD_MAX) error("successor_third() x too big\n"); return (3*x - 1); } static inline value_t successor_third_or_cap (value_t x) { return (x > SUCCESSOR_THIRD_MAX ? VALUE_MAX : successor_third (x)); } /* -------------------------------------------------------------------------- */ /* State */ struct stack { value_t x; int s; /* number of steps which has reached x */ }; struct best { #if WANT_WAYS uint32_t ways : 24; uint32_t min_steps : 8; #else uint8_t min_steps; #endif }; /* best[x].min_steps = STEPS_UNKNOWN when no way to reach x yet seen, so minimum steps unknown */ #define STEPS_UNKNOWN 0xFF /* in 8 bit field */ static struct state { struct best best[MAX_WATCH+1]; struct stack stack[MAX_STEPS + 3]; size_t saved_sp; value_t saved_x; int saved_s; uint64_t num_x; uint64_t num_best_changes; double cputime; } state; static uint64_t prev_num_x; static uint64_t prev_num_best_changes; void init_state (void) { size_t i; for (i=0; i <= MAX_WATCH; i++) state.best[i].min_steps = STEPS_UNKNOWN; state.saved_x = 1; state.saved_s = 0; } /* -------------------------------------------------------------------------- */ /* Save and Restore */ /* will_prune_x() takes a stack entry x at s steps. Return 1 if it's no improvement over best[x] already known and so will be pruned when its time comes. */ static int will_prune_x (value_t x, int s) { return (x <= MAX_WATCH && state.best[x].min_steps < s + WANT_WAYS); } static volatile sig_atomic_t flag_interrupt = 0; static volatile sig_atomic_t flag_terminate = 0; const char saved_state_filename[] = "saved-state.data"; const char temp_filename[] = "saved-state.data.tempfile"; const char progress_filename[] = "progress.txt"; void save_state (value_t x, int s, struct stack *sp) { flag_interrupt = 0; assert (0 <= s && s <= MAX_STEPS); state.saved_x = x; state.saved_s = s; state.saved_sp = sp - state.stack; double this_cputime = elapsed_cputime(); state.cputime += this_cputime; /* Write first to temp_filename and then rename, so atomic replacement. */ FILE *fp = fopen_or_die (temp_filename, "wb"); fwrite_or_die (&state, sizeof(state), 1, fp, temp_filename); fclose_or_die (fp, temp_filename); rename_or_die (temp_filename, saved_state_filename); /*--------------*/ /* progress.txt */ fp = fopen_or_die (temp_filename, "w+"); fprintf(fp, "MAX_STEPS %d, MAX_WATCH %d, WANT_WAYS %d\n", MAX_STEPS, MAX_WATCH, WANT_WAYS); fprintf(fp, "at x=%s s=%d, stack %zu entries\n", value_str(state.saved_x), state.saved_s, state.saved_sp); fprintf(fp, "num_x %"PRIu64", best[] changes %"PRIu64"\n", state.num_x, state.num_best_changes); fprintf(fp, "total CPU time %.2lf hours\n", state.cputime / 3600.0); uint64_t this_num_x = state.num_x - prev_num_x; uint64_t this_num_best_changes = state.num_best_changes - prev_num_best_changes; prev_num_x = state.num_x; prev_num_best_changes = state.num_best_changes; fprintf(fp, "this chunk %.1lf seconds, %"PRIu64" x, %"PRIu64" changes\n", this_cputime, this_num_x, this_num_best_changes); double x_per_second = (this_cputime ? this_num_x / this_cputime : -1); double million_x_per_second = (x_per_second >= 0 ? x_per_second / 1e6 : -1); fprintf(fp, " is %.2lf million x/second\n", million_x_per_second); fprintf(fp, "\n"); fprintf(fp, "stack contents (%zu entries)\n", state.saved_sp); size_t i; for (i = 0; i < state.saved_sp; i++) { const char *note = ""; if (will_prune_x (state.stack[i].x, state.stack[i].s)) note = " (will prune)"; fprintf(fp, " x=%s s=%d%s\n", value_str(state.stack[i].x), state.stack[i].s, note); } fclose_or_die (fp, temp_filename); rename_or_die (temp_filename, progress_filename); if (flag_terminate) exit(0); } void restore_state (void) { const char *filename = saved_state_filename; if (option_verbose) printf("resume from %s\n", filename); FILE *fp = fopen_or_die (filename, "rb"); fread_or_die (&state, sizeof(state), 1, fp, filename); fexpect_eof_or_die (fp, filename); fclose_or_die (fp, filename); prev_num_x = state.num_x; prev_num_best_changes = state.num_best_changes; } void sig_handler_save_state (int signum) { if (option_verbose) 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. */ static void init_save (void) { 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 = 15*60; /* 15 minutes */ t.it_interval.tv_usec = 0; t.it_value = t.it_interval; setitimer_or_die(ITIMER_REAL, &t, NULL); elapsed_cputime(); /* dummy first elapsed */ } /* ------------------------------------------------------------------------ */ /* Searching */ int check_consistency (void) { size_t x; for (x = 1; x <= MAX_WATCH; x++) { #if WANT_WAYS assert (state.best[x].ways >= (state.best[x].min_steps==STEPS_UNKNOWN ? 0 : 1)); #endif /* doublings ensure 2*x takes <= s+1 steps */ size_t y = 2*x; if (y <= MAX_WATCH) assert (state.best[y].min_steps <= state.best[x].min_steps + 1); } return (1); } /* Return the smallest x which requires s steps, so first where a[x] = s. Or return 0 if no such x. */ static size_t smallest_of_s (int s) { size_t x; for (x=1; x <= MAX_WATCH; x++) if (state.best[x].min_steps == s) return(x); return (0); } static void output (void) { printf("MAX_STEPS %d, MAX_WATCH %d\n", MAX_STEPS, MAX_WATCH); state.cputime += elapsed_cputime(); printf("CPU time %.1lf seconds, %.1lf million x /second\n", state.cputime, (state.cputime >= 1 ? state.num_x / state.cputime /1e6 : -1)); printf("considered %"PRIu64" x values\n", state.num_x); printf(" made %"PRIu64" state.best[] changes, mean %.4lf per x\n", state.num_best_changes, (state.num_x >= 1 ? (double) state.num_best_changes / state.num_x : -1)); { const char *filename = "b261870.txt"; FILE *fp = fopen_or_die (filename, "w"); unsigned max_steps = 0; size_t n; for (n = 1; n <= MAX_WATCH && state.best[n].min_steps <= MAX_STEPS; n++) { fprintf(fp,"%zu %d\n", n, state.best[n].min_steps); check_A261870_num_steps (n, state.best[n].min_steps); max_steps = MAX(max_steps, state.best[n].min_steps); } printf("%s = min steps, terms to n=%zu, max %u\n", filename, n-1, max_steps); ferror_die (fp, filename, "writing"); fclose_or_die (fp, filename); } { #if WANT_WAYS const char *filename = "b378375.txt"; FILE *fp = fopen_or_die (filename, "w"); uint64_t max_ways = 0; size_t n; for (n = 1; n <= MAX_WATCH && state.best[n].min_steps <= MAX_STEPS; n++) { fprintf(fp,"%zu %d\n", n, state.best[n].ways); check_A378375_ways (n, state.best[n].ways); max_ways = MAX(max_ways, state.best[n].ways); } printf("%s = num ways, terms to n=%zu, max %"PRIu64"\n", filename, n-1, max_ways); if (n == MAX_WATCH + 1) printf(" which is all of MAX_WATCH\n"); ferror_die (fp, filename, "writing"); fclose_or_die (fp, filename); #endif } { const char *filename = "b378376.txt"; FILE *fp = fopen_or_die (filename, "w"); int s; for (s = 0; s <= MAX_STEPS; s++) { size_t x = smallest_of_s(s); if (x==0) break; fprintf(fp, "%d %zu\n", s, x); check_A378376_smallest_of_s(s,x); } printf("%s = smallest of steps, terms to n=%u\n", filename, s-1); if (s > MAX_STEPS ) printf(" which is all of MAX_STEPS\n"); else printf(" which is < MAX_STEPS (try increasing MAX_WATCH)\n"); ferror_die (fp, filename, "writing"); fclose_or_die (fp, filename); } } static value_t can_reach_watch_table[MAX_STEPS+1]; static void init_can_reach_watch_table (void) { value_t x = MAX_WATCH; int more; for (more = 0; more < numberof(can_reach_watch_table); more++) { int s = MAX_STEPS - more; can_reach_watch_table[s] = x; x = successor_third_or_cap(x); } if (option_verbose >= 2) { int s; for (s=MAX_STEPS; s >= 0; s--) { value_t x = can_reach_watch_table[s]; if (x != VALUE_MAX) printf("can_reach_watch_table[%d] = %s\n", s, value_str(x)); } } } /* can_reach_watch(x,s) takes x value reached after s steps. Return 1 if x is small enough that further steps up to MAX_STEPS might be able to drop by enough to become <= MAX_WATCH. */ static inline int can_reach_watch (value_t x, int s) { assert (0 <= s && s < numberof(can_reach_watch_table)); return (x <= can_reach_watch_table[s]); } /* fill_doublings() takes x is reached at s steps b = &state.best[x], old_s = b->min_steps Have already checked s <= old_s or s <= old_s (per WANT_WAYS), so this b changes. Then look at further steps 2*x,s+1, etc. */ static inline void fill_doublings (size_t x, int s, struct best *b, int old_s) { do { DEBUG(printf(" consider best x=%s s=%d vs %d\n", value_str(x), s, old_s)); assert (x <= MAX_WATCH); assert (s < STEPS_UNKNOWN); assert (b == &state.best[x]); assert (old_s == b->min_steps); assert (WANT_WAYS ? (s <= old_s) : (s < old_s)); #if WANT_WAYS if (s == old_s) { if (UNLIKELY (++b->ways == 0)) error("state.best[].ways data type overflow, at x=%zu\n", x); if (UNLIKELY (option_verbose >= 2)) printf("another x=%zu steps %d, to ways=%d\n", x, s, b->ways); } else #endif { if (UNLIKELY (option_verbose >= 2)) { if (b->min_steps == STEPS_UNKNOWN) printf("initial x=%zu steps %d\n", x, s); else printf("improve x=%zu steps %d was %d\n", x, s, old_s); } b->min_steps = s; #if WANT_WAYS b->ways = 1; #endif } state.num_best_changes++; size_t new_x = x << 1; b = &state.best[new_x]; s++; if (UNLIKELY (x > MAX_WATCH/2)) break; x = new_x; old_s = b->min_steps; /* Want s < old_s (or <= for WANT_WAYS). Anything bigger means no improvement at this x. Also means no improvement at 2*x etc since those 2*x etc were filled at successive s+1 etc here previously. */ } while (WANT_WAYS ? (s <= old_s) : (s < old_s)); } /* push_predecessor_thirds() pushes onto state.stack[] all predecessors (x*2^e + 1) / 3 which are integers. x is at s steps, and the predecessors are at s+1+e steps. Return the new stack pointer sp. */ static inline struct stack * push_predecessor_thirds (value_t x, int s, struct stack *sp) { assert (x >= 1); assert (s <= MAX_STEPS); if (UNLIKELY (s == MAX_STEPS)) return (sp); s++; value_t y = x * VALUE_DIV3_MUL; assert ((y <= VALUE_DIV3_CMP_0) == (x%3==0)); assert ((y >= VALUE_DIV3_CMP_1) == (x%3==1)); assert ((VALUE_DIV3_CMP_0 < y && y < VALUE_DIV3_CMP_1) == (x%3==2)); if (UNLIKELY (y <= VALUE_DIV3_CMP_0)) { DEBUG(printf(" x==0 mod 3, no thirds predecessors\n")); return (sp); } if (y >= VALUE_DIV3_CMP_1) { DEBUG(printf(" x==1 mod 3, commence at 2*x=%s, s=%d\n", value_str(2*x),s+1)); x <<= 1; y <<= 1; s++; assert (y == x * VALUE_DIV3_MUL); } y += VALUE_DIV3_MUL; /* to y = (x+1)/3 or (2x+1)/3 */ DEBUG(printf(" x=%s",value_str(x)); printf(" predecessor y=%s\n",value_str(y))); for (;;) { assert (y == predecessor_third(x)); /* y = (x+1)/3 */ assert (successor_third(y) == x); /* 3y-1 = x */ if (UNLIKELY (s > MAX_STEPS)) { break; } if (UNLIKELY(! can_reach_watch(y,s))) { DEBUG(printf(" stop for y=%s,s=%d cannot reach MAX_WATCH\n", value_str(y), s)); break; } DEBUG(printf(" push x=%s", value_str(x)); printf(" predecessor y=%s s=%d\n", value_str(y), s)); assert (state.stack <= sp && sp < state.stack + numberof(state.stack)); sp->x = y; sp->s = s; sp++; /* y increase += x since (4*P*x+1)/3 - (P*x+1)/3 = P*x */ y += x; x <<= 2; s += 2; } return (sp); } static void check_configuration (void) { if (VALUE_BITS % 2 != 0) error("Oops, expected VALUE_BITS to be even.\n" "Is it a highly unusual data type with an odd number of bits?\n" "VALUE_DIV3_MUL etc are only setup for even number of bits.\n"); if (MAX_STEPS + 2 > VALUE_BITS) error ("MAX_STEPS %d too big for VALUE_BITS %d,\n" "change value_t to a bigger type and re-compile\n", MAX_STEPS, VALUE_BITS); { struct best b; b.min_steps = STEPS_UNKNOWN; if (b.min_steps != STEPS_UNKNOWN) error("STEPS_UNKNOWN apparently bigger than best[].min_steps field."); b.min_steps++; if (b.min_steps != 0) error("STEPS_UNKNOWN apparently no maximum for best[].min_steps field."); if (MAX_STEPS >= STEPS_UNKNOWN) error ("MAX_STEPS %d too big for STEPS_UNKNOWN = %d.\n" "Increase best[].min_steps field, set STEPS_UNKNOWN, and re-compile\n", MAX_STEPS, (int) STEPS_UNKNOWN); } { int fill_s = 0; size_t x; for (x=1; x <= MAX_WATCH; x <<= 1) fill_s++; if (MAX_STEPS + fill_s >= STEPS_UNKNOWN) error ("oops, MAX_STEPS + fill implied by MAX_WATCH overflows state.best[].min_steps field\n"); } } static void search (void) { /* sp points at the next unused element of the stack[] */ struct stack *sp = state.stack + state.saved_sp; value_t x = state.saved_x; int s = state.saved_s; for (;;) { if (UNLIKELY(flag_interrupt)) save_state(x,s,sp); state.num_x++; DEBUG(printf("at x=%s s=%d, num_x %"PRIu64"\n", value_str(x), s, state.num_x)); assert (check_consistency()); if (UNLIKELY(x <= MAX_WATCH)) { struct best *b = &state.best[x]; int old_s = b->min_steps; if (LIKELY(WANT_WAYS ? (s > old_s) : (s >= old_s))) { DEBUG(printf(" seen before %d steps, prune\n", old_s)); goto POP_STACK; } fill_doublings (x,s, b,old_s); } sp = push_predecessor_thirds(x,s,sp); POP_STACK: if (UNLIKELY (sp == state.stack)) break; sp--; x = sp->x; s = sp->s; DEBUG(printf("pop to stack %tu entries\n", sp - state.stack)); } if (state.stack[numberof(state.stack)-1].x != 0) error("oops, stack[] array used beyond intended extent\n"); } static void show_configuration (void) { static int done = 0; if (done) return; done = 1; if (WANT_ASSERT || option_verbose) printf("assert()s %s\n", WANT_ASSERT ? "enabled" : "disabled"); if (option_verbose) { printf("value_t %d bits, %zu bytes\n", VALUE_BITS, sizeof(value_t)); printf("MAX_STEPS %d, stack[] %zu bytes each is %zu bytes\n", MAX_STEPS, sizeof(state.stack[0]), sizeof(state.stack)); printf("MAX_WATCH %d, best[] %zu bytes each is %.2lf mbytes\n", MAX_WATCH, sizeof(struct best), sizeof(state.best)/1024.0/1024.0); struct best b; memset(&b, '\xFF', sizeof(b)); int bits; for (bits=0; b.min_steps != 0; bits++, b.min_steps >>= 1); printf(" best[].min_steps bit field %d bits\n", bits); #if WANT_WAYS for (bits=0; b.ways != 0; bits++, b.ways >>= 1); printf(" best[].ways bit field %d bits\n", bits); #endif printf("state size %zu bytes\n", sizeof(state)); } } int main (int argc, char **argv) { setbuf(stdout,NULL); DEBUG (option_verbose = 2); int option_resume = 0; int i; for (i = 1; i < argc; i++) { const char *arg = argv[i]; if (strcmp(arg,"-v") == 0) { option_verbose++; } else if (strcmp(arg,"start")==0) { option_resume = 0; } else if (strcmp(arg,"resume")==0) { option_resume = 1; } else { error("unrecognised command line option: %s\n", arg); } } show_configuration(); init_can_reach_watch_table(); check_configuration(); if (option_resume) restore_state(); else init_state(); init_save(); search(); output(); return (0); }