/* Kevin Ryde, August 2024 Usage: ./a.out [-v] [start] [resume] This is some C code, using GNU MP and best on a >= 64 bit build, to find terms of A372729(n) = smallest k for which Fibonacci(k) begins and ends with n decimal digits making a partial palindrome zyxw ... wxyz Output is a free-form report for each n >= 3, showing k and the matching starting and ending digits of Fibonacci(k), n=9 found k=306312948 677718655 ... 556817776 Optional command line parameter "-v" prints some more verbose notes on the search and its progress. The method is to step through k keeping some low digits and some upper and lower bounds on high digits of Fibonacci(k). GNU MP is used to recalculate high digits when the bounds become too wide, and to resume the search at a particular k. Save and Resume --------------- The target n, and next k to consider, are saved in a file saved-state.txt at 5 minute intervals and on kill (SIGINT, SIGTERM, SIGHUP). This file can be re-read and the search resumed by ./a.out resume Only the n and k parts of saved-state.txt are used. The rest is progress information. The n and k can be edited to resume at a particular point if desired. n = 12 k = 441183729637 Compile Time ------------ Compile-time FG_MINIMUM_DIGITS is the minimum number of decimal digits to calculate at the high and low ends of Fibonacci(k). Digits are packed 5 bits each in machine words (unsigned long). WIDTH is the number of words needed for FG_MINIMUM_DIGITS. The supplied FG_MINIMUM_DIGITS = 24 is WIDTH = 2 words of 12 digits each on a 64 bit system. A 32 bit system can be used too and in that case 24 digits is 4 words of 6 digits each. (Reducing FG_MINIMUM_DIGITS to 18, so 3 words, might be a speedup for small n in that case.) WIDTH_DIGITS is the actual number of digits used in the WIDTH words. Must have WIDTH_DIGITS > the target n digits, and in general want some headroom above the target n. Each step widens the bounds by about +1, so say 7 digits headroom might expect a recalculate of the high digits after maybe 10^6 steps. More WIDTH is slower steps but recalculate less often. Algorithm --------- fmod is the least significant Fibonacci digits, so for instance in the intended WIDTH_DIGITS = 24 digits, fmod = Fibonacci(k) mod 10^24 fu,fl are upper and lower bounds on the most significant 24 digits fu*10^x >= Fibonacci(k) >= fl*10^x with 10^24 > fu >= 10^23, so fu is 24 digits long The scale factor x is not calculated as such, just ensure fu has the most significant digits of Fibonacci(k) in the most significant digit of fu. The high digits where fu and fl agree are known with certainty. Some of those digits, reversed, may match the digits in fmod. If >= n of them match then a(n) = k. All of these Fibonacci numbers are stepped to k+1 by maintaining also g = Fibonacci(k-1) in similar form gu,gl,gmod, so that [fu,gu] <- [fu + gu, gu] [fl,gl] <- [fl + gl, gl] [fmod,gmod] <- [(fmod + gmod) mod 10^24, gmod] For the low digits, (fmod + gmod) mod 10^24 just means ignoring any carry out when adding. For the high digit bounds, if a step makes fu >= 10^24, then all of fu,fl,gu,gl are reduced by dividing /10. Fibonacci growth means at most one /10 is needed. In that division, upper bounds fu and gu round upwards and lower bounds fl and gl round downwards. These roundings widen the interval between fu and fl and eventually fu,fl don't have enough common digits. In that case, recalculate fu,fl,gu,gl by a higher precision double-and-step calculation with GNU MP, which normally results in fu = fl + 1. Digits Form ----------- Decimal digits are held as 5 bits each packed into "fivedec" words. See comments "Fivedec Decimal Digits" for how that helps addition. Quantities fu,fl,gu,gl,fmod,gmod are all WIDTH many words each. Only terms n >= 3 are sought, so that the code doesn't have to notice that Fibonacci number 55 is a full palindrome. The code does notice when the number of matching digits exceeds the known digits, even after recalculate. If that happens then increase FG_MINIMUM_DIGITS to get more known digits. ensure_high_3() checks that fu,fl bounds are tight enough to know at least 3 high digits with certainty, and it re-calculates if not. This means match_is_ge_3() can examine fu alone. match_is_ge_3() compares 3 high digits of fu against 3 low of fmod. reverse3_table[] does the digit reversal needed to compare. Presuming digits are fairly random, only about 1/1000 of k go to the full match check. match_compare_with_recalc() is the full match check. It re-calculates fu,fl,gu,gl if it finds fu,fl are too far apart to find the point where high and low digits differ. Run with "-v" verbose to report when recalculates happen. Decimal Boundaries ------------------ Occasionally high digits of Fibonacci(k) are close to a decimal digit boundary. For example 12 digits, k = 739530013 Fibonacci(k) = 914999999999.2 ... * 10^154552620 fu = 915000000000 fl = 914999999999 If running a reduced FG_MINIMUM_DIGITS = 12 then fu,fl bounds are 12 digits long but only know 2 digits "91" with certainty. ensure_high_3() stops the search. In this particular case, the low digits are fmod = ... 233 which reversed to 332... does not match the high "91" known, so this k is not a(10) (which is sought at that point). There's no attempt to automatically go past this type of case. It should be unlikely in 24 digits and is left to manual attention if it does happen. Speed ----- The key operations are stepping the three f,g forms to k+1, and inspecting the bits. These are simple additions, bit shifts, and some compares. The estimate would be about 50 to 100 cycles per k. A compiler like GNU C with good loop unrolling is important. The fivedec_vector_... functions operate on WIDTH = 2 length vectors and that size is known at compile time and inlined, so should result in two copies of simple operations. Things Not Done --------------- GNU C on a 64 bit system has a __int128 type which could be a single WORD of 25 digits. But current CPUs (eg. x86 family) tend not to have native 128 bit addition, and limited 128 shifts. On a modest x86, two 64 bit WORDs measured faster than one 128 bit. Try 128 by "typedef unsigned __int128 WORD;" if desired. */ /* FG_MINIMUM_DIGITS is the minimum number of decimal digits to have in the f,g Fibonacci numbers. In the intended 64 bit build, 24 is 2 WORDs of 12 each. */ #define FG_MINIMUM_DIGITS 24 /* 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 #include <assert.h> #include <errno.h> #include <inttypes.h> #include <signal.h> #include <stdarg.h> #include <stdio.h> #include <stdlib.h> #include <string.h> #include <stdnoreturn.h> #include <time.h> #include <sys/time.h> #include <gmp.h> /* ------------------------------------------------------------------------ */ /* 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])) /* evaluate to ceil(n/d), with no risk of overflow */ #define DIV_CEIL(n,d) (((n)/(d)) + (((n)%(d)) != 0)) 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 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 (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/1e9; } static int uint64_bit_length (uint64_t k) { /* for progress info, so don't need to be fast */ int len = 0; while (k) { k>>=1; len++; } return (len); } static inline int uint64_count_leading_zeros (uint64_t k) { #ifdef __GNUC__ if (sizeof(unsigned) >= sizeof(k)) return (__builtin_clz(k)); if (sizeof(unsigned long) >= sizeof(k)) return (__builtin_clzl(k)); if (sizeof(unsigned long long) >= sizeof(k)) return (__builtin_clzll(k)); #endif uint64_t mask = (uint64_t)1 << 63; int ret = 0; while (! (k & mask)) { mask>>=1; ret++; } return (ret); } /* ------------------------------------------------------------------------ */ /* Fivedec Decimal Digits. Decimal digits are held in 5 bits each in a WORD of type unsigned long. WORD is expected to be the native CPU word size. At least one bit is spare above the 5s, ready to receive a carry. For example, a 64 bit WORD high low carry five five ... five five \---/ \---------------------/ 4 bits 12*5 = 60 bits When adding two decimal digits x+y, an extra offset +22 turns 0..9 into 22..31 and the latter is ready to have any carry propagate up to the next 5 bit position by ordinary binary addition. After such an offset addition, 5 bit digits are 0..9 if wrap-around (a carry went up) 22..31 if not The 22..31 must drop back down to 0..9. This is indicated by the "16" bit (high bit of the 5). fivedec_adjust() clears that 16 and shifts that 16 to become possible further 6 to subtract. A vector of fivedec words has x[0] least significant up to x[len-1] most significant. */ typedef unsigned long WORD; #define WORD_BITS ((int) (CHAR_BIT * sizeof(WORD))) #define WORD_MAX (~ (WORD)0) #define FIVEDEC_NUMB_DIGITS ((WORD_BITS - 1) / 5) #define FIVEDEC_NUMB_BITS (5*FIVEDEC_NUMB_DIGITS) #define FIVEDEC_ABOVE_BITS (WORD_BITS - FIVEDEC_NUMB_BITS) #define FIVEDEC_NUMB_MASK (WORD_MAX >> FIVEDEC_ABOVE_BITS) #define FIVEDEC_ABOVE_MASK (~ FIVEDEC_NUMB_MASK) /* all digits 1 */ #define FIVEDEC_ONES (FIVEDEC_NUMB_MASK / 31) #define FIVEDEC_ADJUST_MASK (16 * FIVEDEC_ONES) #define FIVEDEC_ADD_OFFSET (22 * FIVEDEC_ONES) /* mask for most significant decimal digit */ #define FIVEDEC_MASK_HIGHEST \ (FIVEDEC_NUMB_MASK ^ (FIVEDEC_NUMB_MASK >> 5)) /* mask for most significant 3 decimal digits */ #define FIVEDEC_MASK_HIGH_3 \ (FIVEDEC_NUMB_MASK ^ (FIVEDEC_NUMB_MASK >> (3*5))) /* mask for least significant 3 digits */ #define FIVEDEC_MASK_LOW_3 (~ (WORD_MAX << (3*5))) /* Check that all digits in x are in the range 0..9, and 0 bits above the FIVEDEC_NUMB_BITS of them (no carry left). */ int fivedec_validate (WORD x) { int i; for (i=0; i < FIVEDEC_NUMB_DIGITS; i++) { int d = (x >> (5*i)) & 31; if (! (0 <= d && d <= 9)) { fprintf(stderr, "oops, i=%d digit %d out of range\n", i,d); return (0); } } WORD above = x >> FIVEDEC_NUMB_BITS; if (above) { fprintf(stderr, "oops, above bits %d\n", (int) above); return (0); } return (1); } /* Check digits are in range. Carry part is ignored. */ int fivedec_validate_noc (WORD x) { return (fivedec_validate (x & FIVEDEC_NUMB_MASK)); } /* Check vector of WORDs x[0..len-1]. If noc=1 then no check of the carry part of high x[len-1] (is allowed to be anything). */ int fivedec_vector_validate (const WORD x[], int len, int noc) { assert (len >= 0); int ret = 1; int i; for (i=0; i < len; i++) { int r = (i==len-1 && noc ? fivedec_validate_noc (x[i]) : fivedec_validate (x[i])); if (! r) ret = 0; } return (ret); } /* Return the least significant digit from x. */ static inline WORD fivedec_low_digit (WORD x) { /* this is used with some carry part still in, so the noc relaxed validate */ assert (fivedec_validate_noc(x)); return (x & 15); } /* Return the most significant digit from x. Must have carry part = 0. */ static inline WORD fivedec_high_digit (WORD x) { assert (fivedec_validate(x)); return (x >> (5*(FIVEDEC_NUMB_DIGITS-1))); } /* Adjust digits of x back into 0..9 range after an addition. Any carry in x is unchanged. */ static inline WORD fivedec_adjust (WORD x) { /* Each 5 bits in x are 0..9 unchanged, or 22..31 -> 0..9. mask is bit value 16 and that bit is always cleared. Then possible subtract 3*(16/8) = 6. */ WORD mask = FIVEDEC_ADJUST_MASK; WORD flags = x & mask; x &= ~mask; return (x - 3*(flags >> 3)); } /* Adjust digits of x back into 0..9 range after an addition. Clear any carry in x too. */ static inline WORD fivedec_adjust_and_clear_carry (WORD x) { WORD flags = x & FIVEDEC_ADJUST_MASK; x &= FIVEDEC_NUMB_MASK - FIVEDEC_ADJUST_MASK; return (x - 3*(flags >> 3)); } /* Return the carry in x after an addition. This is 0 or 1. x digits don't have to be adjusted yet. */ static inline WORD fivedec_carry (WORD x) { WORD c = x >> FIVEDEC_NUMB_BITS; assert (c==0 || c==1); return (c); } /* Return x+y, with any carry left above the digits. Any carries already in x or y add into the resulting carry too (but there's only certain to be 1 bit space there). */ static inline WORD fivedec_add (WORD x, WORD y, WORD cin) { assert (fivedec_validate_noc(x)); assert (fivedec_validate_noc(y)); assert ((cin)==0 || (cin)==1); return (fivedec_adjust (x + y + cin + FIVEDEC_ADD_OFFSET)); } /* Set *dst = x+y+cin digits part, and return any carry. */ static inline WORD fivedec_addc (WORD *dst, WORD x, WORD y, WORD cin) { assert (fivedec_validate(x)); assert (fivedec_validate(y)); assert (cin==0 || cin==1); WORD s = (x + y + cin + FIVEDEC_ADD_OFFSET); *dst = fivedec_adjust_and_clear_carry(s); return (fivedec_carry(s)); } /* Add vectors of decimal digits x[0..len-1] + y[0..len-1]. Store result in dst[0..len-1]. Leave any carry in high word dst[len-1]. */ static inline void fivedec_vector_add (WORD dst[], const WORD x[], const WORD y[], int len) { int i; assert (len >= 0); WORD c = 0; for (i=0; i < len-1; i++) c = fivedec_addc (&dst[i], x[i],y[i], c); dst[i] = fivedec_add (x[i], y[i], c); } /* x += add for a single fivedec word "add". Leave any carry in the high WORD x[len-1]. */ static inline void fivedec_vector_add_single (WORD x[], int len, WORD add) { assert (fivedec_vector_validate(x,len, 1)); assert (len >= 1); assert (0 <= add && add <= 9); /* Always add into low word x[0]. For len>=2, any carry beyond that is very unlikely. */ int i; for (i=0; i < len-1; i++) { add = fivedec_addc (&x[i], x[i], add, 0); if (LIKELY (add==0)) return; } assert (i == len-1); x[i] = fivedec_add (x[i], add, 0); } /* x[0..len-1] /= 10, meaning shift decimal digits down. Any carry in high word x[len-1] shifts down to become new high digit. If round_up=1 then division is round upwards by adding +1 if low digit had been !=0. */ static inline void fivedec_vector_div10 (WORD x[], int len, int round_up) { assert (fivedec_vector_validate(x,len, 1)); assert (round_up==0 || round_up==1); WORD d = 0; /* low digit from the word above */ int i; for (i=len-1; i >= 0; i--) { WORD w = x[i]; x[i] = (w>>5) + (d << (5*(FIVEDEC_NUMB_DIGITS-1))); d = fivedec_low_digit(w); } if (round_up) { WORD c = LIKELY(d!=0); /* possible +1 for ceil */ fivedec_vector_add_single (x, len, c); assert (fivedec_carry(x[len-1]) == 0); } } static inline WORD fivedec_vector_get_digit (const WORD x[], int len, int pos) { assert (0 <= pos && pos < len*FIVEDEC_NUMB_DIGITS); int wpos = pos / FIVEDEC_NUMB_DIGITS; int dpos = pos % FIVEDEC_NUMB_DIGITS;; return ((x[wpos] >> (5*dpos)) & 31); } /* Return 1 if x > y, 0 if x = y, or -1 if x < y. */ int fivedec_vector_cmp (WORD x[], WORD y[], int len) { int i; for (i=len-1; i>=0; i--) { if (x[i] > y[i]) return(1); if (x[i] < y[i]) return(-1); } return (0); } static inline void fivedec_vector_copy (WORD *dst, const WORD *src, int len) { memcpy (dst, src, len*sizeof(WORD)); } /* Set dst[0..len-1] to the decimal digits of z. The value in z is changed. */ void fivedec_vector_from_mpz (WORD dst[], int len, mpz_t z) { /* this is not critical for speed, so a simple form is enough */ int i, dpos; assert (len >= 0); for (i=0; i<len; i++) { WORD w = 0; for (dpos=0; dpos < FIVEDEC_NUMB_DIGITS; dpos++) { WORD digit = mpz_fdiv_q_ui (z, z, 10); w += digit << (dpos*5); } dst[i] = w; } assert (fivedec_vector_validate(dst,len, 0)); if (mpz_sgn(z) != 0) error ("fivedec_vector_from_mpz() mpz too big for %d words\n", len); } void fivedec_print_digit (int d) { if (d <= 9) printf("%d", d); else printf("[%d]", d); } void fivedec_vector_print (const WORD x[], int len) { int i, nz=0; for (i = len*FIVEDEC_NUMB_DIGITS - 1; i>=0; i--) { int digit = fivedec_vector_get_digit(x,len,i); if (!nz && i!=0 && digit==0) continue; /* skip leading 0s */ fivedec_print_digit(digit); nz = 1; } } /* ------------------------------------------------------------------------ */ /* Save State */ static int option_verbose = 0; static int option_maxn = INT_MAX; static const char saved_state_filename[] = "saved-state.txt"; static const char temp_filename[] = "saved-state.txt.tempfile"; static double state_last_cputime; static uint64_t state_last_count; static uint64_t state_last_recalculates; void save (int n, uint64_t k) { /* Write first to temp_filename and then rename, so atomic replacement. */ FILE *fp = fopen_or_die (temp_filename, "w"); fprintf(fp, "n = %d\n", n); fprintf(fp, "k = %"PRIu64"\n", k); /* extra info */ fprintf(fp, "\n"); fprintf(fp, "k is %d bits\n", uint64_bit_length(k)); double t = cputime(); double duration = t - state_last_cputime; state_last_cputime = t; fprintf(fp, "last chunk %.1f seconds CPU time\n", duration); double per_second = (duration >= 1 ? state_last_count / duration : -1.0); fprintf(fp, " %"PRIu64" many k is %.2f million/second\n", state_last_count, per_second/1e6); fprintf(fp, " %"PRIu64" recalculates\n", state_last_recalculates); state_last_count = 0; ferror_die (fp, temp_filename, "writing"); fclose_or_die (fp, temp_filename); rename_or_die (temp_filename, saved_state_filename); if (option_verbose >= 2) printf("save to %s n=%d k=%"PRIu64"\n", saved_state_filename, n, k); } 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); state_last_cputime = cputime(); } /* ------------------------------------------------------------------------ */ void check_A372729_first_match_n (int n, uint64_t k) { static uint64_t want[] = { 0, /* no n=0 */ 0, /* n=1 */ 10, /* n=2 */ 317, /* n=3 */ 1235, /* n=4 */ 28898, /* n=5 */ 120742, /* n=6 */ 1411753, /* n=7 */ 201095722, /* n=8 */ 306312948, /* n=9 677718655 ... 556817776 */ 12306316582, /* n=10 1958298604 ... 4068928591 */ 32679761048, /* n=11 10120617339 ... 93371602101 */ 806327047899, /* n=12 622440432113 ... 311234044226 */ 9600042921304, /* n=13 3060350410283 ... 3820140530603 */ 172192972068022 /* n=14 11718839590715 ... 51709593881711 */ }; if (! (n >= 1)) error("oops, check_A372729_first_match_n() is for n>=1\n"); if (n >= numberof(want)) { if (option_verbose) printf(" (beyond A372729 sample data)\n"); return; } if (k == want[n]) { if (option_verbose) printf(" (k good vs A372729 sample data)\n"); return; } error("oops n=%d want k=%"PRIu64" got %"PRIu64"\n", n, want[n], k); } /* ------------------------------------------------------------------------ */ /* number of WORDs per f,g Fibonacci numbers */ #define WIDTH DIV_CEIL(FG_MINIMUM_DIGITS, FIVEDEC_NUMB_DIGITS) #define WIDTH_BYTES (WIDTH * sizeof(WORD)) #define WIDTH_DIGITS (WIDTH * FIVEDEC_NUMB_DIGITS) static mpz_t z_10pow_table[5 * WIDTH_DIGITS]; static mpz_t z_fmod, z_gmod; static mpz_t z_f_squared, z_g_squared; static mpz_t z_fu, z_fl, z_gu, z_gl; static mpz_t z_temp; /* return a pointer to an mpz equal to 10^e */ const mpz_t * z_10pow (int e) { if (UNLIKELY (! (0 <= e && e < numberof(z_10pow_table)))) error ("10^%d out of range 0 .. %zu\n", e, numberof(z_10pow_table) - 1); return (&z_10pow_table[e]); } /* fg_double() takes f = Fibonacci(k), g = Fibonacci(k-1), and sets f2 = Fibonacci(2*k), g2 = Fibonacci(2*k-1). f,g are read before f2,g2 are set so can have for instance f2=f and g2=g for in-place operation. */ void fg_double (mpz_t f2, mpz_t g2, const mpz_t f, const mpz_t g) { mpz_mul (z_f_squared, f, f); mpz_mul (z_g_squared, g, g); mpz_mul (z_temp, f, g); /* f2 = 2*F(k)*F(k-1) + F(k)^2) */ mpz_mul_2exp (f2, z_temp, 1); mpz_add (f2, f2, z_f_squared); /* g2 = F(k)^2 + F(k-1)^2 */ mpz_add (g2, z_f_squared, z_g_squared); } /* fg_double_and_step() takes f = Fibonacci(k), g = Fibonacci(k-1), and sets f2 = Fibonacci(2*k + bit), g2 = Fibonacci(2*k-1 + bit). So 2k,2k-1 if bit=0 or 2k+1,2k if bit=1. f,g are read before f2,g2 are set so can have for instance f2=f and g2=g for in-place operation. */ void fg_double_and_step (mpz_t f2, mpz_t g2, const mpz_t f, const mpz_t g, uint64_t bit) { if (bit == 0) { fg_double (f2,g2, f,g); } else { fg_double (g2,f2, f,g); mpz_add (f2, f2,g2); } /* As a note, the two squares method for double and step is F(2k+1) = 4*F(k)^2 - F(k-1)^2 + 2*(-1)^k F(2k-1) = F(k)^2 + F(k-1)^2 F(2k) = F(2k+1) - F(2k-1) But here k_to_fg_high_z() runs with f,g as upper or lower bounds instead of full values, so doesn't know the low place where +/-2 should go and whether any carry or borrow would propagate up. */ } void init_fg_calculate (void) { mpz_inits (z_fmod, z_gmod, z_f_squared, z_g_squared, z_fu, z_fl, z_gu, z_gl, z_temp, NULL); int i; for (i=0; i < numberof(z_10pow_table); i++) { mpz_init (z_10pow_table[i]); mpz_ui_pow_ui (z_10pow_table[i], 10, i); } } /* Set z_fmod = Fibonacci(k) mod M z_gmod = Fibonacci(k-1) mod M, where modulus M = 10^WIDTH_DIGITS. */ void k_to_fgmod_z (uint64_t k) { if (UNLIKELY( ! (k >= 2) )) error("k_to_fgmod_z() is for k>=1\n"); mpz_set_ui (z_fmod, 1); /* F(1) = 1 */ mpz_set_ui (z_gmod, 0); /* F(0) = 0 */ int shift = 62 - uint64_count_leading_zeros(k); assert (shift >= 0); uint64_t mask; for (mask = (uint64_t)1 << shift; mask != 0; mask>>=1) { DEBUG (printf("mask=%"PRIu64, mask); gmp_printf(" fmod %Zd gmod %Zd\n", z_fmod, z_gmod)); uint64_t bit = k & mask; fg_double_and_step (z_fmod,z_gmod, z_fmod,z_gmod, bit); mpz_mod (z_fmod, z_fmod, *z_10pow(WIDTH_DIGITS)); mpz_mod (z_gmod, z_gmod, *z_10pow(WIDTH_DIGITS)); } DEBUG(gmp_printf("mpz fmod %Zd gmod %Zd\n", z_fmod, z_gmod)); } /* Set fivedec vectors fmod[] = Fibonacci(k) mod M, and gmod[] = Fibonacci(k-1) mod M, where M = 10^WIDTH_DIGITS, which means the least significant decimal digits which fit into WIDTH words. */ void k_to_fgmod (uint64_t k, WORD fmod[WIDTH], WORD gmod[WIDTH]) { k_to_fgmod_z (k); fivedec_vector_from_mpz (fmod,WIDTH, z_fmod); fivedec_vector_from_mpz (gmod,WIDTH, z_gmod); DEBUG(printf("fivedec fmod "); fivedec_vector_print(fmod,WIDTH); printf(" gmod "); fivedec_vector_print(gmod,WIDTH); printf("\n")); } /* Decimal digits length of z. Must have z >= 1. */ int z_decimals_length (const mpz_t z) { if (mpz_sgn(z) <= 0) error("z_decimals_length() of <= 0\n"); /* sizeinbase is correct or 1 too big. If z < 10^(len-1) then len-1 is the length. */ int len = mpz_sizeinbase(z,10); if (len > 0 && mpz_cmp(z, *z_10pow(len-1)) < 0) len--; DEBUG(gmp_printf("z = %Zd sizeinbase = %d to length %d\n", z, (int) mpz_sizeinbase(z,10), len)); assert (mpz_cmp(*z_10pow(len-1), z) <= 0 && mpz_cmp(z, *z_10pow(len)) < 0); return (len); } /* Divide down all of z_fu,z_fl,z_gu,z_gl by 10^e. Rounded up for "u" and down for "l" upper and lower bounds respectively. */ void fg_high_div_10pow (int e) { const mpz_t *p = z_10pow(e); DEBUG(gmp_printf(" from fu = %Zd\n", z_fu); gmp_printf(" divide %Zd\n", *p)); mpz_cdiv_q (z_fu, z_fu, *p); mpz_fdiv_q (z_fl, z_fl, *p); mpz_cdiv_q (z_gu, z_gu, *p); mpz_fdiv_q (z_gl, z_gl, *p); } /* z_fu,z_fl,z_gu,z_gl all multiplied up by 10^e */ void fg_high_mul_10pow (int e) { const mpz_t *p = z_10pow(e); mpz_mul (z_fu, z_fu, *p); mpz_mul (z_fl, z_fl, *p); mpz_mul (z_gu, z_gu, *p); mpz_mul (z_gl, z_gl, *p); } /* Set global z_fu etc to bounds fu*10^x >= Fibonacci(k) >= fl*10^x gu*10^x >= Fibonacci(k-1) >= gl*10^x where exponent x leaves fu having exactly WIDTH_DIGITS. (Most significant digit at most significant place in fu.) x itself is not calculated, but is either positive or negative as needed to have fu length WIDTH_DIGITS. */ void k_to_fg_high_z (uint64_t k) { /* The upper bounds here are f = (Fibonacci(k) + something)/10^x g = (Fibonacci(k-1) + something)/10^x fg_double() and fg_double_and_step() use formulas which are increasing functions of both f and g so upper bounds f and g at k result in upper bounds at new 2k+1,2k,2k-1. Similarly the lower bounds f = (Fibonacci(k) - something)/10^x etc. */ if (UNLIKELY( ! (k >= 2) )) error("k_to_fg_high_z() is for k>=1\n"); mpz_set_ui (z_fu, 1); /* F(1) = 1 */ mpz_set_ui (z_gu, 0); /* F(0) = 0 */ mpz_set (z_fl, z_fu); mpz_set (z_gl, z_gu); int shift = 62 - uint64_count_leading_zeros(k); assert (shift >= 0); uint64_t mask; for (mask = (uint64_t)1 << shift; mask != 0; mask>>=1) { DEBUG(printf("at mask=0x%"PRIx64"\n", mask); gmp_printf(" fu = %Zd\n", z_fu)); assert (mpz_cmp (z_fu,z_fl) >= 0); assert (mpz_cmp (z_gu,z_gl) >= 0); assert (mpz_cmp_ui (z_fl, 0) > 0); assert (mpz_cmp_ui (z_gl, 0) >= 0); uint64_t bit = k & mask; fg_double_and_step (z_fu,z_gu, z_fu,z_gu, bit); fg_double_and_step (z_fl,z_gl, z_fl,z_gl, bit); /* target about 3*WIDTH_DIGITS */ int discard = mpz_sizeinbase(z_fu,10) - 3*WIDTH_DIGITS; if (discard > 0) fg_high_div_10pow (discard); } /* increase or decrease fu so exactly WIDTH_DIGITS long */ int discard = z_decimals_length(z_fu) - WIDTH_DIGITS; DEBUG(printf(" final length %d discard %d\n", z_decimals_length(z_fu), discard)); if (discard > 0) fg_high_div_10pow (discard); else if (discard < 0) fg_high_mul_10pow (-discard); if (option_verbose) { printf("k=%"PRIu64" interval gap ", k); mpz_sub (z_temp, z_fu,z_fl); gmp_printf("%Zd\n", z_temp); gmp_printf(" fu = %Zd\n", z_fu); gmp_printf(" fl = %Zd\n", z_fl); } } void k_to_fg_high (uint64_t k, WORD fu[WIDTH], WORD fl[WIDTH], WORD gu[WIDTH], WORD gl[WIDTH]) { k_to_fg_high_z (k); assert (mpz_cmp(z_fu,z_fl) >= 0); assert (mpz_cmp(z_gu,z_gl) >= 0); assert (mpz_cmp(z_fu,z_gl) >= 0); fivedec_vector_from_mpz (fu,WIDTH, z_fu); fivedec_vector_from_mpz (fl,WIDTH, z_fl); fivedec_vector_from_mpz (gu,WIDTH, z_gu); fivedec_vector_from_mpz (gl,WIDTH, z_gl); if (fivedec_vector_get_digit(fu,WIDTH, WIDTH_DIGITS-1) == 0) error("oops, fu high digit = 0\n"); } void output (int n, uint64_t k, const WORD fu[WIDTH], const WORD fl[WIDTH], const WORD fmod[WIDTH]) { printf("n=%d found k=%"PRIu64"\n", n, k); int i; printf(" "); for (i = WIDTH_DIGITS - 1; i >= WIDTH_DIGITS - n; i--) fivedec_print_digit(fivedec_vector_get_digit(fu,WIDTH, i)); printf(" ... "); for (i = n-1; i >= 0; i--) fivedec_print_digit(fivedec_vector_get_digit(fmod,WIDTH, i)); printf("\n"); if (option_verbose) { printf(" fu = "); fivedec_vector_print(fu,WIDTH); printf("\n"); printf(" fl = "); fivedec_vector_print(fl,WIDTH); printf("\n"); printf(" fmod = "); fivedec_vector_print(fmod,WIDTH); printf("\n"); } check_A372729_first_match_n (n,k); printf("\n"); } /* Return 1 if fu,fl have 3 high digits in common, so those digits known with certainty. */ static inline int have_high_3 (const WORD fu[WIDTH], const WORD fl[WIDTH]) { assert (fivedec_vector_validate (fu,WIDTH, 0)); assert (fivedec_vector_validate (fl,WIDTH, 0)); return ((fu[WIDTH-1] & FIVEDEC_MASK_HIGH_3) == (fl[WIDTH-1] & FIVEDEC_MASK_HIGH_3)); } /* Ensure fu,fl have 3 high digits in common by a recalculate if necessary. */ static inline void ensure_high_3 (int n, uint64_t k, WORD fu[WIDTH], WORD fl[WIDTH], WORD gu[WIDTH], WORD gl[WIDTH]) { assert (FIVEDEC_NUMB_DIGITS >= 3); if (UNLIKELY(! have_high_3(fu,fl))) { if (option_verbose) printf("recalculate for 3 digits precision at n=%d k=%"PRIu64"\n", n,k); state_last_recalculates++; k_to_fg_high (k, fu,fl, gu,gl); if (UNLIKELY(! have_high_3(fu,fl))) { printf("WIDTH_DIGITS = %d exceeded\n", WIDTH_DIGITS); printf(" at k=%"PRIu64" still don't have 3 known digits after recalculate\n", k); printf (" fu = "); fivedec_vector_print(fu,WIDTH); printf("\n"); printf (" fl = "); fivedec_vector_print(fl,WIDTH); printf("\n"); exit(1); } } } /* Set *matchp to the number of matching digits between fu,fmod. Return 0 if this was successful. Return 1 if the bounds fu,fl don't know enough to find where fu,fmod differ (and in that case *matchp is how far reached). */ static int match_compare (int *matchp, uint64_t k, const WORD fu[WIDTH], const WORD fl[WIDTH], const WORD fmod[WIDTH]) { /* w_fu is a word from fu and is successively shifted up "<<= 5" so the relevant digit to consider is at the highest position. w_fl similarly from fl. w_fmod is a word from fmod and is successively shift down ">>= 5" so the relevant digit to consider is at the lowest position. */ assert (fivedec_vector_validate (fu,WIDTH, 0)); assert (fivedec_vector_validate (fl,WIDTH, 0)); assert (fivedec_vector_validate (fmod,WIDTH, 1)); int wi,i, match = 0; for (wi=0; LIKELY(wi < WIDTH); wi++) { WORD w_fu = fu[WIDTH-1 - wi]; WORD w_fl = fl[WIDTH-1 - wi]; WORD w_fmod = fmod[wi]; for (i=0; LIKELY (i < FIVEDEC_NUMB_DIGITS); i++) { DEBUG(printf(" digits fmod=%d fu=%d fl=%d\n", (int) (w_fmod & 15), (int) ((w_fu >> (5*(FIVEDEC_NUMB_DIGITS-1))) & 15), (int) ((w_fl >> (5*(FIVEDEC_NUMB_DIGITS-1))) & 15))); if (UNLIKELY( (w_fu ^ w_fl) & FIVEDEC_MASK_HIGHEST )) { /* fu,fl different digits, so don't know enough */ *matchp = match; return (1); } WORD w_fu_high_digit = w_fu >> (5*(FIVEDEC_NUMB_DIGITS-1)); WORD diff = w_fmod ^ w_fu_high_digit; if (diff & 15) { DEBUG(printf(" different, match length is %d\n", match)); *matchp = match; return (0); } w_fmod >>= 5; w_fu <<= 5; w_fl <<= 5; match++; } DEBUG(printf(" matching in next WORD\n")); } error (">= %d matching digits\n", WIDTH_DIGITS); } /* reverse_table[t] takes t = x*2^10 + y*2^5 + z and have value which is their 5 bits reversal z*2^10 + y*2^5 + x. Should have 0 <= x,y,z <= 9, though y and z can be up to their full 5 bits values. */ static uint16_t reverse3_table[10 << 10]; static void init_reverse3_table (void) { int i; for (i = 0; i < numberof(reverse3_table); i++) { uint16_t high = (i >> 10) % 32; uint16_t mid = (i >> 5) % 32; uint16_t low = (i >> 0) % 32; reverse3_table[i] = (low<<10) + (mid<<5) + high; } } /* Return 1 if 3 digits of fu and fmod match, so partial palindromes of length >=3. */ static inline int match_is_ge_3 (const WORD fu[WIDTH], const WORD fmod[WIDTH]) { assert (fivedec_vector_validate (fu,WIDTH, 0)); assert (fivedec_vector_validate (fmod,WIDTH, 1)); WORD u = fu[WIDTH-1] >> (5*(FIVEDEC_NUMB_DIGITS - 3)); assert (u < numberof(reverse3_table)); return (reverse3_table[u] == (fmod[0] & FIVEDEC_MASK_LOW_3)); } static int match_compare_with_recalc (int n, uint64_t k, WORD fu[WIDTH], WORD fl[WIDTH], WORD gu[WIDTH], WORD gl[WIDTH], const WORD fmod[WIDTH]) { int match; if (UNLIKELY (match_compare(&match, k, fu,fl, fmod))) { if (option_verbose) { printf("precision exhausted at k=%"PRIu64", recalculate\n", k); printf (" match >= %d digits\n", match); printf (" fu = "); fivedec_vector_print(fu,WIDTH); printf("\n"); printf (" fl = "); fivedec_vector_print(fl,WIDTH); printf("\n"); printf (" fmod = "); fivedec_vector_print(fmod,WIDTH); printf("\n"); } state_last_recalculates++; k_to_fg_high (k, fu,fl, gu,gl); if (UNLIKELY (match_compare(&match, k, fu,fl, fmod))) { printf ("precision still exceeded after recalc at n=%d k=%"PRIu64"\n", n,k); printf (" match >= %d digits\n", match); printf (" fu = "); fivedec_vector_print(fu,WIDTH); printf("\n"); printf (" fl = "); fivedec_vector_print(fl,WIDTH); printf("\n"); printf (" fmod = "); fivedec_vector_print(fmod,WIDTH); printf("\n"); exit(1); } } return (match); } static inline void fmod_step (WORD fmod[WIDTH], WORD gmod[WIDTH]) { WORD e[WIDTH]; fivedec_vector_add (e, fmod, gmod, WIDTH); fivedec_vector_copy (gmod, fmod, WIDTH); fivedec_vector_copy (fmod, e, WIDTH); } static inline void fhigh_step (WORD fu[WIDTH], WORD fl[WIDTH], WORD gu[WIDTH], WORD gl[WIDTH]) { { WORD el[WIDTH]; fivedec_vector_add (el, fl,gl, WIDTH); fivedec_vector_copy (gl, fl, WIDTH); fivedec_vector_copy (fl, el, WIDTH); } { WORD eu[WIDTH]; fivedec_vector_add (eu, fu,gu, WIDTH); fivedec_vector_copy (gu, fu, WIDTH); fivedec_vector_copy (fu, eu, WIDTH); } WORD c = fu[WIDTH-1] & FIVEDEC_ABOVE_MASK; if (UNLIKELY (c)) { DEBUG(printf(" rshift\n")); fivedec_vector_div10 (fu,WIDTH, 1); fivedec_vector_div10 (fl,WIDTH, 0); fivedec_vector_div10 (gu,WIDTH, 1); fivedec_vector_div10 (gl,WIDTH, 0); } } void show_configuration (void) { static int done = 0; if (option_verbose && ! done) { done = 1; #ifdef __GNUC__ printf("Have __GNUC__\n"); #else printf("Not __GNUC__\n"); #endif printf("WORD is %zu bytes, %d bits\n", sizeof(WORD), WORD_BITS); printf("fivedec %d digits per WORD, 5 bits each, %d bits above\n", FIVEDEC_NUMB_DIGITS, WORD_BITS - FIVEDEC_NUMB_BITS); printf("FG_MINIMUM_DIGITS = %d becomes WIDTH = %d words and WIDTH_DIGITS = %d\n", FG_MINIMUM_DIGITS, WIDTH, WIDTH_DIGITS); printf("reverse3_table[] %zu entries, %zu bytes\n", numberof(reverse3_table), sizeof(reverse3_table)); printf("assert()s %s\n", WANT_ASSERT ? "enabled" : "disabled"); } else { assert (printf("assert()s enabled\n")); } } /* Search for terms a(n) onward, starting at k as next candidate. */ void search_from (int n, uint64_t k) { WORD fu[WIDTH], fl[WIDTH], gu[WIDTH], gl[WIDTH]; WORD fmod[WIDTH], gmod[WIDTH]; show_configuration (); /* Only handle search for n>=3, and that means k>=11. These are past single-digit matches, and the 2 digits whole match Fibonacci(10) = 55. */ if (n < 3) n = 3; if (k < 11) k = 11; k_to_fg_high (k, fu,fl, gu,gl); k_to_fgmod (k, fmod, gmod); if (option_verbose) printf("searching n=%d at k = %"PRIu64"\n", n, k); for (;;) { if (UNLIKELY(flag_interrupt)) { flag_interrupt = 0; save (n, k); if (flag_terminate) return; } DEBUG(printf("k=%"PRIu64"\n", k)); DEBUG(printf(" fu ");fivedec_vector_print(fu,WIDTH);printf("\n")); DEBUG(printf(" fl ");fivedec_vector_print(fl,WIDTH);printf("\n")); DEBUG(printf(" fmod ");fivedec_vector_print(fmod,WIDTH);printf("\n")); assert (fivedec_vector_validate(fu,WIDTH, 0)); assert (fivedec_vector_validate(fl,WIDTH, 0)); assert (fivedec_vector_validate(gu,WIDTH, 0)); assert (fivedec_vector_validate(gl,WIDTH, 0)); assert (fivedec_vector_validate(fmod,WIDTH, 1)); assert (fivedec_vector_validate(gmod,WIDTH, 1)); assert (fivedec_vector_cmp(fu,fl,WIDTH) >= 0); assert (fivedec_vector_cmp(gu,gl,WIDTH) >= 0); ensure_high_3 (n,k, fu,fl, gu,gl); if (UNLIKELY(match_is_ge_3(fu,fmod))) { int match = match_compare_with_recalc(n,k, fu,fl,gu,gl, fmod); while (match >= n) { output (n,k, fu,fl, fmod); if (n >= option_maxn) return; n++; } } fmod_step (fmod, gmod); fhigh_step (fu,fl, gu,gl); if (UNLIKELY(k==UINT64_MAX)) error("exceeded 64-bit k\n"); k++; state_last_count++; } } void resume (void) { if (option_verbose) printf("resume from %s\n", saved_state_filename); FILE *fp = fopen_or_die (saved_state_filename, "r"); int n; uint64_t k; if (fscanf(fp, "n = %d\nk = %"SCNu64"\n", &n, &k) != 2) error("%s format unrecognised\n", saved_state_filename); ferror_die (fp, saved_state_filename, "reading"); fclose_or_die (fp, saved_state_filename); search_from (n, k); } int main (int argc, char *argv[]) { int i, endpos; DEBUG (option_verbose = 1); setbuf(stdout,NULL); init_save(); init_reverse3_table(); init_fg_calculate(); state_last_cputime = cputime(); int option_resume = 0; 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 if (sscanf(arg, "maxn=%d%n", &option_maxn, &endpos) == 1 && endpos == strlen(arg)) { } else { error("Unrecognised command line option: %s\n", arg); } } if (option_resume) resume(); else search_from (3, 11); return(0); }