/* Kevin Ryde, June 2023 Usage: ./a.out [-v] [nk n k] Find terms of A357425, which is the smallest number with sum of digits n in fractional base 4/3. Output is b-file style so a(100)=4142502 is 100 4142502 The strategy is to maintain the base 4/3 expansion of a candidate k. If its sum of digits is below the target n then increase at digit position i where the digits at and above i are too small to reach n, even if all digits below i were 3s (the maximum). Save File --------- Search state n and k are saved in a file saved-state.txt on program interrupt (SIGINT, SIGTERM, SIGHUP), and periodically on a timer. The next program run reads this and resumes the search at that point so a long run can be done in pieces. Remove the file to start again. Command Line ------------ Optional parameter "-v" prints results and progress with greater verbosity. Optional parameter "nk" starts the search at a given target n and prospective k. For example ./a.out nk 100 4000000 k can be any value <= the true a(n) and the search starts at k. The n and k values in saved-state.txt can be edited in the same way. Implementation Notes -------------------- Digits of candidate k>=0 are held in the e[] array. The least significant digit of k is omitted and the rest are e[0] = second least, e[1] = third least, and so on. e[i].digit base 4/3 digit, value 0 .. 3 e[i].value (value of digits e[i..MAX]) /3 e[i].delta 3*(i+1) - sum digits e[i..MAX] In base 4/3, taking just some high digits is a value which is a multiple of 3, ready for (4/3)*high + digit for the next lower digit. e[i].value is without this factor of 3. e[i].delta is how much the sum of digits would exceed n if digits e[0..i-1], and the absent least significant digit, were all 3s. Value Type Size --------------- Set VALUE_BITS to the desired bit width for type value_t which is used for e[i].value. The conditionals select among unsigned unsigned long unsigned long long GNU MP limbs (link with -lgmp) A 64-bit type overflows in about 15 minutes with a good compiler and even a modest CPU (at the time of writing). So look to a 128 bit long long on a 64-bit system, or say 96 bits (3 words) with GNU MP on a 32-bit system. Type selection is by preprocessor conditionals on UINT_WIDTH etc defines from limits.h. Digit Increases --------------- If e[i].delta < 0 then sum of digits n cannot be reached by digits in positions i-1 and below, even all were digit 3. This means something at and above i must change. The highest i position with e[i].delta < 0 is found and an increase is made there. e[i].delta is an increasing function of i and the initial START_SOMEWHERE code uses that to find highest delta < 0. The increase is +3 at position i. This amount ensures the high part remains a multiple of 3 (as noted above). Carries propagate up in the usual way: If a digit is >=4 then -4 there and carry +3 to the digit above it. Carry ends at a 0 digit (since +3 on it becomes 3). e[i].value is held divided /3, so +3 means increase +1 there. On reaching the top of carry propagation, new "delta"s are established back down again with the new digits, e[j].delta = e[j+1].delta + e[j].digit - 3 e[j+1] is always available by having end-most e[ELEMS_LEN-1] "unused", in the sense that carry is not allowed to propagate into it. If some delta < 0 is found, then that position is the new highest delta < 0. Stop filling deltas and increase again there. Digit Down ---------- When an increase finds no new delta < 0, all e[i..MAX] are good digits, values, and deltas > 0. A new digit is set at e[i-1] by value = e[i].value 3*q + r = value divide by 3, remainder 0 <= r <= 2 if r=0 new_digit = 0 new_value = value + q if r!=0 new_digit = 3-r new_value = value + q + 1 e[i-1].digit = new_digit e[i-1].value = new_value new_digit is the smallest consistent with digits i-1 and above having a value which is a multiple of 3. value and new_value are in e[] style with /3 divided out. The calculation with *3 included is 3*new_value = (4/3)*(3*value) + new_digit which simplifies to new_value = value + (value + new_digit) / 3 If value != 0 mod 3 then new_digit has only a single way to ensure an integer "(value + new_digit) / 3". If value == 0 mod 3 then either 0 or 3 would make integer "(value + new_digit) / 3". Taking the smaller, ie. 0, is the smallest increase in k and so will find the smallest k with digit sum n. It can be noted that taking a 3 digit is a bigger digit sum at that point. But sometimes 0 permits digit combinations below with the same (or larger) overall total. Compiler Help ------------- BASE_NUMERATOR and BASE_DENOMINATOR (4/3) are compile-time constants to let the compiler see them for operator strength reduction. e[] is a static global with the idea that its location can be a constant for comparing ptr==e, or ptr==e_second_last. It may be necessary to disable position-independent executable options for this. In recent versions of gcc use -no-pie. (Or just -fno-pie for constant in the code and let the loader relocate when run.) Things Not Done --------------- e[i].delta is based on maximum digit 3s below. Such digits can and do occur, but often the actual maximum is smaller. e[i].value mod 3^x suffices to determine possible digit combinations in the next x many digits. Maybe a table of remainders to maximum digit sum would reduce "delta"s enough to apply increases at higher places. An increase carries up to the next 0 digit, so a higher place has to be above that to actually help. Maybe e[i].value could be held in ternary digits so its "mod 3" is ready for new_digit. Some bit twiddling might do ternary carry propagation by normal word additions. */ /* Set WANT_ASSERT to 1 to enable self-checks (slow). Set WANT_DEBUG to 1 for some rough development prints. */ #define WANT_ASSERT 0 #define WANT_DEBUG 0 #if ! WANT_ASSERT #define NDEBUG #endif /* for UINT_WIDTH etc from limits.h */ #define __STDC_WANT_IEC_60559_BFP_EXT__ #include <assert.h> #include <ctype.h> #include <errno.h> #include <inttypes.h> #include <limits.h> #include <float.h> #include <signal.h> #include <stdarg.h> #include <stdint.h> #include <stdio.h> #include <stdlib.h> #include <string.h> #include <stdnoreturn.h> #include <time.h> #include <unistd.h> #include <sys/time.h> #define BASE_NUMERATOR 4 /* fractional base 4/3 */ #define BASE_DENOMINATOR 3 #define DIGIT_MAX (BASE_NUMERATOR - 1) #define VALUE_BITS 128 /* Uncomment USE_GMP to force GNU MP on any VALUE_BITS, rather than the default use-if-needed. (Run with verbosity "-v" to see final type and sizes.) */ /* #define USE_GMP 1 */ /* ------------------------------------------------------------------------ */ /* Generic */ #if WANT_DEBUG #define DEBUG(expr) do { expr; } while (0) #else #define DEBUG(expr) do { } while (0) #endif #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 #define numberof(array) (sizeof(array)/sizeof((array)[0])) noreturn ATTRIBUTE_PRINTF void error (const char *format, ...) { va_list ap; va_start (ap, format); vfprintf (stderr, format, ap); va_end (ap); exit(1); } 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); } void ferror_die (FILE *fp, const char *filename, const char *action) { if (ferror(fp)) error("Error %s %s: %s\n", action, filename, strerror(errno)); } void fclose_or_die (FILE *fp, const char *filename) { if (fclose(fp) != 0) error("Error closing %s\n", filename); } 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)); } 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)); } 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, in seconds */ double cputime (void) { struct timespec t; if (clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &t) != 0) error("Cannot clock_gettime() CPUTIME: %s\n", strerror(errno)); return t.tv_sec + t.tv_nsec/10e9; } /* ------------------------------------------------------------------------ */ #if ! USE_GMP #ifndef UINT_WIDTH #error Oops, UINT_WIDTH not defined #endif #if VALUE_BITS <= UINT_WIDTH typedef unsigned value_t; const char value_t_name[] = "unsigned"; #define PRIvalue "u" #elif VALUE_BITS <= ULONG_WIDTH typedef unsigned long value_t; const char value_t_name[] = "unsigned long"; #define PRIvalue "lu" #elif VALUE_BITS <= ULLONG_WIDTH typedef unsigned long long value_t; const char value_t_name[] = "unsigned long long"; #define PRIvalue "Lu" #elif VALUE_BITS <= 64 typedef uint64_t value_t; const char value_t_name[] = "uint64_t"; #define PRIvalue PRIu64 #elif ! defined(USE_GMP) #define USE_GMP 1 #else #error No type big enough for VALUE_BITS, and USE_GMP set to no #endif #endif #if USE_GMP #include <gmp.h> #define VALUE_NUM_LIMBS \ ((VALUE_BITS + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS) typedef struct value_t { mp_limb_t limbs[VALUE_NUM_LIMBS]; } value_t; const char value_t_name[] = "mp_limb_t[]"; #else #define VALUE_MAX (~ (value_t) 0) #endif static value_t value_zero; #define ASSERT_DIGIT(d) \ do { assert(0 <= (d) && (d) <= DIGIT_MAX); } while (0) struct elem { value_t value; int delta; char digit; }; static inline void value_overflow (int bool, const char *where) { if (UNLIKELY(bool)) error("value overflow (in %s)\n", where); } /* Print value *p as decimal digits. */ void value_fprint (FILE *fp, const value_t *p) { #if USE_GMP gmp_fprintf(fp, "%Nu", p->limbs, VALUE_NUM_LIMBS); #else fprintf(fp, "%"PRIvalue, *p); #endif } void value_print (const value_t *p) { value_fprint(stdout,p); } /* Add "*p + add" and put the result back in *p. */ static inline void value_add_small (value_t *p, unsigned add) { #if USE_GMP value_overflow (mpn_add_1(p->limbs, p->limbs, VALUE_NUM_LIMBS, add) !=0, "add small"); #else value_overflow (*p > VALUE_MAX - add, "add"); *p += add; #endif } /* Add "*x + *y" and put the result back in *x. */ static inline void value_add (value_t *x, value_t *y) { #if USE_GMP value_overflow (mpn_add_n(x->limbs, x->limbs, y->limbs, VALUE_NUM_LIMBS) !=0, "add"); #else value_overflow (*x > VALUE_MAX - *y, "add"); *x += *y; #endif } /* Multiply "*src * mul" and put the result in *dst. Can have src==dst for an in-place operation. */ static inline void value_mul_small (value_t *dst, const value_t *src, unsigned mul) { #if USE_GMP mp_limb_t c; if (mul == 4) c = mpn_lshift (dst->limbs, src->limbs, VALUE_NUM_LIMBS, 2); else c = mpn_mul_1(dst->limbs, src->limbs, VALUE_NUM_LIMBS, mul); value_overflow (c!=0, "mul"); #else value_overflow (*src > VALUE_MAX / mul, "mul"); *dst = *src * mul; #endif } /* Divide "*src / div". Set *dst to quotient and return remainder. Can have src==dst for an in-place operation. */ static inline unsigned value_divrem (value_t *dst, value_t *src, unsigned div) { #if USE_GMP return (mpn_divrem_1(dst->limbs, 0, src->limbs, VALUE_NUM_LIMBS, div)); #else unsigned r = *src % div; *dst = *src / div; return(r); #endif } /* Return 1 if values *x == *y. */ static inline int value_is_equal (const value_t *x, const value_t *y) { return (memcmp(x, y, sizeof(*x)) == 0); } /* str is a null-terminated character string of decimal digits. Set *p to their value, or error exit if bad string. */ void value_from_str (value_t *p, const char *str) { /* For integer types, don't use sscanf() as it doesn't report overflow. Could use mpz_set_str() or similar when GMP, but may as well same code all value_t types. */ int c; *p = value_zero; while ((c = *str++)) { if (! isdigit(c)) error("value_from_str() invalid digit %d \"%c\"\n", c,c); c -= '0'; assert (0 <= c && c <= 9); value_mul_small(p,p,10); value_add_small(p, c); } } /* Return the sum of digits of k in base 4/3 (A244041). */ int value_sumdigits_43 (value_t k) { int sum_digits = 0; while (! value_is_equal(&k,&value_zero)) { sum_digits += value_divrem(&k, &k, BASE_NUMERATOR); value_mul_small (&k, &k, BASE_DENOMINATOR); } return (sum_digits); } int value_bit_length (value_t k) { int len = 0; while (! value_is_equal(&k,&value_zero)) { value_divrem(&k, &k, 2); len++; } return (len); } /* ------------------------------------------------------------------------ */ /* command line option */ int option_verbose = 0; /* length of the digits[] array */ #define ELEMS_LEN 1000 static int n; static struct elem e[ELEMS_LEN]; static const struct elem * const e_second_last = e + ELEMS_LEN - 2; /* Number of base 4/3 digits represented by e[], inclusive of the least significant digit not stored in e[]. */ int elems_num_digits (void) { int i; for (i=ELEMS_LEN-1; i>=0; i--) if (e[i].digit!=0) return(i+2); return(1); } /* Print the base 4/3 digits in e[] and additional "low_digit". */ void elems_digits_print (int low_digit) { int len = elems_num_digits(); int i; printf("base %d/%d len=%d\n", BASE_NUMERATOR, BASE_DENOMINATOR, len); printf(" digits = "); for (i=len-2; i>=0; i--) { printf ("%d", e[i].digit); ASSERT_DIGIT (e[i].digit); } ASSERT_DIGIT (low_digit); printf ("%d\n", low_digit); } /* elems_print(pos) is a diagnostic print of e[pos.. upwards]. Lower positions e[0..pos-1] are junk and are printed as 'x". */ void elems_print (int pos) { int len = ELEMS_LEN; int i; while (len > 0 && e[len-1].digit==0) len--; if (len <ELEMS_LEN) len++; printf(" n=%d base %d/%d num digits %d\n", n, BASE_NUMERATOR, BASE_DENOMINATOR, len-1); printf(" digits = "); for (i=0; i<len; i++) { if (i>=pos) printf ("%d%s", e[i].digit, i==len-1?"":","); else printf ("x%s", i==len-1?"":","); } printf("\n"); printf(" deltas = "); for (i=0; i<len; i++) { if (i>=pos) printf ("%d%s", e[i].delta, i==len-1?"":","); else printf ("x%s", i==len-1?"":","); } printf("\n"); printf(" values = "); for (i=0; i<len; i++) { if (i>=pos) value_print(&e[i].value); else printf ("x"); if (i < len-1) printf (","); } printf ("\n");; } /* Consistency checks on e[pos..ELEMS_LEN-1] when seeking sum of digits n and where ptr = &e[pos]. Return 1 if good or error exit if bad. */ int elems_validate (const struct elem *ptr) { int i; int delta = DIGIT_MAX*(ELEMS_LEN+1) - n; value_t k = value_zero; /* value represented e[] digits */ if (ptr < e) error ("ptr below e"); if (ptr > e_second_last) error ("ptr above second-last of e"); int pos = ptr - e; for (i = ELEMS_LEN-1; i >= pos; i--) { int digit = e[i].digit; if (! (0 <= digit && digit <= DIGIT_MAX)) { fprintf(stderr, "\ndigit i=%d bad: %d\n", i, digit); elems_print(pos); exit(1); } delta += digit - DIGIT_MAX; if (delta != e[i].delta) { fprintf(stderr, "\ni=%d sum digits got %d want %d\n", i, e[i].delta, delta); elems_print(pos); exit(1); } /* k *= 4/3, and to do that must have k divisible by 3 */ int r = value_divrem (&k, &k, BASE_DENOMINATOR); if (r!=0) error("i=%d not divisible by %d", i, BASE_DENOMINATOR); value_mul_small (&k, &k, BASE_NUMERATOR); value_add_small (&k, digit); /* e[i].value is k/3, and k must be divisible by 3 */ value_t k3 = k; r = value_divrem (&k3, &k3, BASE_DENOMINATOR); if (r!=0) error("i=%d after digit not divisible by %d", i, BASE_DENOMINATOR); if (! value_is_equal(&k3, &e[i].value)) { printf("i=%d value got ", i); value_print(&e[i].value); printf(" want "); value_print(&k3); printf("\n"); exit(1); } } return (1); } /* Fill e[] from the given k. The least significant digit of k is discarded. The rest go in e[], with second lowest at e[0]. e[].delta values use n. */ void elems_from_value (value_t k) { int i; value_divrem(&k, &k, BASE_NUMERATOR); for (i=0; i < ELEMS_LEN; i++) { e[i].value = k; /* before multiply by 3 */ value_mul_small (&k, &k, BASE_DENOMINATOR); e[i].digit = value_divrem(&k, &k, BASE_NUMERATOR); } if (! value_is_equal(&k, &value_zero)) { printf("too big for ELEMS_LEN = %d (leftover ", ELEMS_LEN); value_print(&k); printf(")\n"); exit(1); } int delta = DIGIT_MAX*(ELEMS_LEN+1) - n; for (i=ELEMS_LEN-1; i>=0; i--) e[i].delta = (delta += e[i].digit - DIGIT_MAX); } /* ptr is a pointer into e[] which has a good value etc. Find the smallest next digit (lower position) possible and set *(ptr-1) accordingly. If ptr->value == 0 mod 3, then next digit can be 0 or 3 and 0 is set. Otherwise only one choice. */ static inline int elems_down (struct elem *ptr) { DEBUG(printf("pos=%d down\n", ptr-e)); assert (ptr > e); int delta = ptr->delta - BASE_DENOMINATOR; int r = value_divrem(&(ptr-1)->value, &ptr->value, BASE_DENOMINATOR); DEBUG(printf(" q="); value_print(&(ptr-1)->value); printf(" r=%d\n",r)); value_add (&(ptr-1)->value, &ptr->value); int digit = 0; if (UNLIKELY(r)) { digit = BASE_DENOMINATOR - r; /* digit+r multiple of 3 */ value_add_small (&(ptr-1)->value, 1); } DEBUG(printf(" new digit %d to value ", digit); value_print(&(ptr-1)->value); printf("\n")); (ptr-1)->digit = digit; delta += digit; (ptr-1)->delta = delta; return (delta); } /* ------------------------------------------------------------------------ */ /* Save State */ const char saved_state_filename[] = "saved-state.txt"; const char temp_filename[] = "saved-state.txt.tempfile"; double state_last_cputime; static void save (struct elem *ptr) { /* Write first to temp_filename and then rename, so atomic replacement. */ for ( ; ptr > e; ptr--) elems_down (ptr); value_t k; value_mul_small (&k, &e[0].value, BASE_NUMERATOR); FILE *fp = fopen_or_die (temp_filename, "w"); fprintf(fp, "n = %d\n", n); fprintf(fp, "k = "); value_fprint(fp,&k); fprintf(fp,"\n"); fprintf(fp, "\n(k is %d bits, %d base %d/%d digits, sum digits %d)\n", value_bit_length(k), elems_num_digits(), BASE_NUMERATOR, BASE_DENOMINATOR, value_sumdigits_43(k)); double t = cputime(); double duration = t - state_last_cputime; fprintf(fp, "(%.1f seconds CPU time since last save)\n", duration); state_last_cputime = t; ferror_die (fp, temp_filename, "writing"); fclose_or_die (fp, temp_filename); rename_or_die (temp_filename, saved_state_filename); if (option_verbose) { printf("save to %s n=%d k=", saved_state_filename, n); value_print(&k); printf("\n"); } } static volatile int flag_interrupt = 0; static volatile int flag_terminate = 0; 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 hourly 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 = 3600; /* 1 hour */ t.it_interval.tv_usec = 0; t.it_value = t.it_interval; setitimer_or_die(ITIMER_REAL, &t, NULL); state_last_cputime = cputime(); } /* ------------------------------------------------------------------------ */ /* Some initial values, as a consistency check. */ const char * const want_str[] = { "0", /* n=0 */ "1", "2", "3", "5", "6", "7", "10", "11", "15", "21", /* n=10 */ "22", "23", "31", "39", "43", "54", "55", "74", "75", "101", /* n=20 */ "102", "103", "138", "139", "183", "187", "246", "247", "330", "331", "439", "443", "587", "783", "790", "791", "1047", "1355", "1398", "1399", "1866", "1867", "2487", "2491", "3318", "3319", "4199", "4427", "5903", "5911", "7882", "7883", "9959", "10571", "13991", "14711", "18791", "23607", "24887", "25063", "31479", "33419", "40439", "41975", "45047", "58999", "62007", "77927", "79223", "82679", "110239", "127819", "146987", "170427", "187895", "227238", "227239", "287479", "302987", "331511", "414711", "450550", "450551", "552951", "718198", "718199", "737271", "983029", "983030", "983031", "1310709", "1310710", "1310711", "1747615", "2330151", "2330155", "3106870", "3106871", "4142495", "4142502", /* n=100 */ "4142503", "5523338", "5523339", "6000631", "7364454", "7364455", "9819274", "9819275", "13092343", "17456458", "17456459", "17649399", "23275318", "23275319", "23532535", "31033767", "31376715", "41378358", "41378359", "41835623", "55171147", "59938727", "73561531", "98082038", "98082039", "99165927", "130776054", "130776055", "132221239", "174368075", "176294987", "235059983", "249543479", "309987703", "413316918", "413316919", "443632855", "551089227", "551089255", "734785638", "734785639", "979714186", "979714187", "1051574183", "1356442471", "1417417719", "1741714199", "1858629623", "1889890295", "2347745015", "2600076619", "3215271051", "3466768827", "4173768919", "4287028071", "4622358439", "5716037430", "5716037431", "7420033639", "7621383243", "8217526119", "10161844326", "10161844327", "10956701495", "13549125771", "17397497590", "17397497591", "18065501031", "24087334709", "24087334710", "24087334711", "32116446282", "32116446283", "42821928375", "42821928379", "57095904502", "57095904503", "76127872671", "76127872678", "76127872679", "101503830231", "135338440309", "135338440310", "135338440311", "180451253749", "180451253750", "180451253751", "240601671669", "240601671670", "240601671671", "320802228895", "427736305191", "427736305195", "570315073590", "570315073591", "760420098039", "760420098123", "1013893464054", "1013893464055", "1013893464167", /* n=200 */ "1351857952075", "1411685285878", "1411685285879", "1882247047839", "2403303025911", "2509662730454", "2509662730455", "3204404034551", "3346216973942", "3346216973943", "4461622631925", "4461622631926", "4461622631927", "5948830175903", "7314760036455", "9753013381941", "9753013381942", "9753013381943", "13004017842591", "13503332227943", "17338690456790", "17338690456791", "23118253942389", "23118253942390", "23118253942391", "30824338589855", "32007898614391", "41099134211403", "42677198152439", "54798845615206", /* n=230 */ "54798845615207", /* n=231 */ "58812458406903" /* n=232 */ }; /* k is a number which has sum of digits n. low_digit is the low digit of k. */ void found (value_t k, int low_digit) { if (option_verbose) { printf ("n=%d found a(n)=", n); value_print(&k); printf(" "); elems_digits_print(low_digit); printf("\n"); } else { /* b-file style */ printf ("%d ", n); value_print(&k); printf ("\n"); } elems_validate (e); { int got_sum_digits = value_sumdigits_43(k); if (got_sum_digits != n) error(" oops, sum digits %d != %d = n\n", got_sum_digits, n); } if (n>0 && low_digit==0) error(" oops, low_digit %d cannot be smallest k\n", low_digit); if (n < numberof(want_str)) { value_t want_value; value_from_str (&want_value, want_str[n]); if (! value_is_equal(&k, &want_value)) error(" oops, expected %s\n", want_str[n]); } } /* e[] contains a number which together with some low digit >=1 will reach sum of digits n. If this low digit < 3 then increasing it +1 or +2 reaches sum of digits n+1 or n+2 too. Call found() with the resulting k values and increase n to new target after those now found. */ void elems_found (void) { DEBUG(printf("elems_found()\n")); elems_validate(e); int low_digit = BASE_DENOMINATOR - e[0].delta; if (low_digit < 0) error ("oops, found sum of digits bigger than n\n" " meaning starting k was bigger than a(n)\n"); if (low_digit > DIGIT_MAX) error ("oops, found final digit %d is > %d\n", low_digit, DIGIT_MAX); value_t k; value_mul_small (&k, &e[0].value, BASE_NUMERATOR); value_add_small (&k, low_digit); for ( ; low_digit <= DIGIT_MAX; low_digit++) { found (k, low_digit); n++; int i; for (i=0; i<ELEMS_LEN; i++) e[i].delta--; if (low_digit == DIGIT_MAX) break; value_add_small (&k, 1); } elems_validate(e); } /* Run a search looking for a(n) with global variable n, starting at prospective value k_start. */ void search (value_t k_start) { init_save(); struct elem *ptr = e; elems_from_value (k_start); if (option_verbose) { printf ("search n=%d k=", n); value_print(&k_start); printf(" BASE = %d/%d\n", BASE_NUMERATOR, BASE_DENOMINATOR); #if USE_GMP printf(" VALUE_BITS = %d using %d limbs of %d bits\n", VALUE_BITS, VALUE_NUM_LIMBS, GMP_NUMB_BITS); #else printf(" VALUE_BITS = %d using \"%s\" %d bytes\n", VALUE_BITS, value_t_name, (int) sizeof(value_t)); #endif printf(" e[] at %p, elem struct %d bytes\n", e, (int) sizeof(struct elem)); printf("\n"); } START_SOMEWHERE: /* if ptr->delta good then go down, otherwise search upward for highest bad to increase */ assert (elems_validate (ptr)); if (UNLIKELY(ptr->delta >= 0)) goto DOWN; for (;;) { assert (e <= ptr && ptr <= e_second_last); if ((ptr+1)->delta >= 0) goto INCREASE; if (UNLIKELY(ptr == e_second_last)) /* second last */ error("ELEMS_LEN too small for n=%d\n", n); ptr++; } DOWN: /* ptr->delta is good, go down to next lower digit, and keep going down while still good */ DEBUG(printf("pos=%d down\n", ptr-e); elems_print(ptr-e)); assert (elems_validate (ptr)); for (;;) { assert (ptr->delta >= 0); if (UNLIKELY(ptr == e)) { elems_found (); goto START_SOMEWHERE; } int new_delta = elems_down(ptr); ptr--; if (new_delta < 0) goto INCREASE; } INCREASE: /* ptr->delta is bad (but delta above it is good), increase ptr->digit */ if (UNLIKELY(flag_interrupt)) { flag_interrupt = 0; save(ptr); if (flag_terminate) return; } DEBUG(printf("pos=%d increase\n", ptr-e); elems_print(ptr-e)); assert (elems_validate (ptr)); assert ( ptr ->delta < 0); assert ((ptr+1)->delta >= 0); { struct elem *c = ptr; /* carry pointer */ for (;;) { value_add_small (&c->value, 1); c->digit += BASE_DENOMINATOR; if (UNLIKELY(c->digit < BASE_NUMERATOR)) break; /* carry upwards */ if (UNLIKELY (c >= e_second_last)) error("ELEMS_LEN too small for n=%d\n", n); c->digit -= BASE_NUMERATOR; c++; } /* digits c..ptr inclusive have changed, work down from c to ptr setting new deltas */ DEBUG (printf(" changed %d digits pos %d to %d\n", c-ptr+1, ptr-e, c-e)); assert (c <= e_second_last); int delta = (c+1)->delta; assert (delta >= 0); for (;;) { delta += c->digit - BASE_DENOMINATOR; c->delta = delta; if (UNLIKELY(delta < 0)) { ptr = c; DEBUG(printf("increase more at pos=%d\n", c-e); elems_print(ptr-e)); goto INCREASE; } if (c==ptr) goto DOWN; c--; } } } void resume (void) { if (option_verbose) printf("\nresume from %s\n", saved_state_filename); FILE *fp = fopen_or_die (saved_state_filename, "r"); char *k_str; if (fscanf(fp, "n = %d\nk = %ms\n\n", &n, &k_str) != 2) error("%s format unrecognised\n", saved_state_filename); ferror_die (fp, saved_state_filename, "reading"); fclose_or_die (fp, saved_state_filename); value_t k_start; value_from_str (&k_start, k_str); free(k_str); search(k_start); } int main (int argc, char *argv[]) { int i; setbuf(stdout,NULL); for (i = 1; i < argc; i++) { const char *arg = argv[i]; if (strcmp(arg,"-v")==0) { option_verbose++; } else if (strcmp(arg,"nk")==0) { if (i+2 >= argc) error("missing parameters n and k for nk\n"); n = atoi(argv[i+1]); value_t k_start; value_from_str(&k_start, argv[i+2]); search(k_start); return(0); } else { error("Unrecognised command line option: %s\n", arg); } } if (access(saved_state_filename, F_OK) == -1 && errno == ENOENT) { search(value_zero); /* with n=0 */ } else { resume(); } return(0); }