/* Kevin Ryde, October 2023 Usage: ./a.out This is some C code to find terms A366312(t) = n which is the index n where A366311(n) = t first occurs, for t = 1..9. Output is free-form like seen term 6 at n=47 seen term 7 at n=1031 ... A366311 places terms t on a square spiral with no palindromes length >= 3 anywhere in lines North-South, East-West, and diagonals. The strategy is an explicit calculation of those terms, watching for when each t first occurs. This is fairly specific to the problem, but makes it possible to find t=9 at n=115318369671 in about 30 minutes (on a modest CPU at the time of writing). Implementation Notes -------------------- As described in A366311, when considering a new term in A366311 it's enough to check only for palindromes length 3 or 4 involving the new term. Any longer palindrome would have a length 3 or 4 within it from existing terms, so does not occur. The code keeps state.buf[] as a sliding window of the most recent BUFLEN many terms. Length 3 or 4 requires up to the 3rd earlier spiral loop. BUFLEN is a power of 2 to allow "&" for wraparound. 2^22 = 4 mbytes is enough to find the first occurrence of t=9. Along the bottom of the spiral, going to the right, the lines examined for possible palindromes are (line 1) (line 2) (line 3) A A A B B B C C C (line 0) A B C P ----> spiral direction In each line "A B C P", must have P != B as otherwise a palindrome B C B, and if B=C then must have P != A as otherwise a palindrome A B B A. Disallowed values for P are accumulated by flags which are bits in an integer "disallow", where bit 1 << (X-1) means X disallowed. Table ABC_to_disallow[A][B][C] makes the "disallow" for given line of terms A,B,C. Table disallow_to_P[disallow] is the resulting smallest allowed P. Positions A_pos, B_pos, C_pos are indices into state.buf[] along the "diagonal back" which is line 1 shown above, A_pos ----> B_pos ----> C_pos ----> P_pos ----> The 4 lines are terms at buf[] at positions line 0 = [ P_pos - 3, P_pos - 2, P_pos - 1 ] line 1 = [ A_pos, B_pos, C_pos ] line 2 = [ A_pos + 3 B_pos + 2, C_pos + 1 ] line 3 = [ A_pos + 6 B_pos + 4, C_pos + 2 ] All these positions, and increments of A_pos,B_pos,C_pos,P_pos to step to the next term, are mod BUFLEN so they wrap-around within state.buf[]. Lines where the spiral turns a corner are not as easy. The following diagram shows P at the last general case, and each of 0..10 are the special cases before P' after the corner. * * * P' after corner | | | | * * * 10 | | | | * * * 9 11 corner cases A_pos | | | | corners[0..10] data table A--*--*--A--*--*--A * * 8 | | | --B--*--B--*--B--*--* * 7 | | --C--C--C--*--*--*--* 6 | --P--0--1--2--3--4--5 P = last general case before corner The 11 cells marked 0..10 are handled in a corners[] table. For that table, P_pos is at the "0" cell shown, having incremented from the last general case. A_pos,B_pos,C_pos are the "diagonal back" relative to there still. The terms to be fetched for each line are relative to one of A_pos,B_pos,C_pos,P_pos. In "struct w_offset", w=0..3 is which one of those and "offset" is the amount to add to it (mod BUFLEN). At Some of the lines go outside the part of the spiral filled so-far. w = -1 marks a no-such-term. w_offset_get_term() returns 0 for that term. ABC_to_disallow[][][] knows 0 means no-such-term and nothing to disallow for that. Parts of buf[] before "A_pos" are not used and under WANT_ASSERT are marked buf[i] = -1 as a consistency check. But not for full-speed running. */ /* Set WANT_ASSERT to 1 to enable self-checks (slower). Set WANT_DEBUG to 1 for some rough development prints. */ #define WANT_ASSERT 0 #define WANT_DEBUG 0 #define WANT_PROGRESS 0 #if ! WANT_ASSERT #define NDEBUG #endif #include #include #include #include #include #include #include #include #include /* ------------------------------------------------------------------------ */ /* 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))) #define PREFETCH_FOR_READ(addr) \ do { __builtin_prefetch (addr, 0, 0); } while(0) #define PREFETCH_FOR_RW(addr) \ do { __builtin_prefetch (addr, 1, 0); } while(0) #else #define LIKELY(cond) (cond) #define UNLIKELY(cond) (cond) #define ATTRIBUTE_PRINTF #define PREFETCH_FOR_READ(addr) #define PREFETCH_FOR_RW(addr) #endif #define MIN(x,y) ((x)<=(y) ? (x) : (y)) #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); } /* ------------------------------------------------------------------------ */ /* BUFLEN must be a power of 2, ie. "1 << somehting". To reach n, need BUFLEN about 11.8*sqrt(n). Conversely a given BUFLEN is roughly n = BUFLEN^2/139. */ #define BUFLEN (1 << 22) #define BUFLEN_MASK (BUFLEN - 1) static int option_verbose = 0; static struct state { uint64_t n; char buf[BUFLEN]; int side_length, side_increment; uint64_t first_occurrence[10]; /* for t=1..9 */ int saved_A_pos; int saved_B_pos; int saved_C_pos; int saved_P_pos; } state; /* T_pos is an index into state.buf[]. Return T_pos + offset with wraparound. */ static inline int add (int T_pos, int offset) { assert (0 <= T_pos && T_pos < BUFLEN); return (T_pos + offset) & BUFLEN_MASK; } /* T_pos is a position in state.buf[]. Return its "n". P_pos is the place in the buffer of state.n and the distance back to T_pos (with wraparound) is the distance back from state.n. */ uint64_t T_pos_to_n (int T_pos, int P_pos) { assert (0 <= T_pos && T_pos < BUFLEN); assert (0 <= P_pos && P_pos < BUFLEN); return (state.n - ((P_pos - T_pos) & BUFLEN_MASK)); } /* When WANT_ASSERT, state.buf[] entries before A_pos are set to -1 to mark those outside the intended extent. */ int zap_buf_entry_prev_T = -1; int zap_buf_entry (int T_pos) { DEBUG(printf(" zap buf at %d\n", T_pos)); assert (0 <= T_pos && T_pos < BUFLEN); assert (zap_buf_entry_prev_T == -1 || T_pos == add(zap_buf_entry_prev_T, 1)); assert (state.buf[T_pos] != -1); state.buf[T_pos] = -1; zap_buf_entry_prev_T = T_pos; return (1); } /* ------------------------------------------------------------------------ */ /* Spiral Coordinates */ /* Set *x = A174344(n), *y = A274923(n) */ void n_to_xy (uint64_t n, int *x, int *y) { assert (n >= 1); int d = ((int) sqrt(n-1) + 1)/2; uint64_t Nbase = 4*(uint64_t)d*(uint64_t)d - 4*(uint64_t)d + 2; int Nrem = n - Nbase; if (0) { DEBUG(printf("n=%"PRIu64" d=%d Nrem=%d\n", n,d,Nrem)); } if (Nrem <= 2*d-1) { *x = d; *y = -(int)d + 1 + Nrem; } else if (Nrem <= 4*d-1) { *x = 3*d-1-Nrem; *y = d; } else if (Nrem <= 6*d-1) { *x = -d; *y = 5*d-1-Nrem; } else { *x = -7*d + 1 + Nrem; *y = -d; } } uint64_t xy_to_n (int x, int y) { if (x > abs(y)) { return ((int64_t) (4*x-2)*x - (x-y) + 1); /* right vertical */ } else if (-x > abs(y)) { return ((int64_t) (4*x-2)*x + (x-y) + 1); /* left vertical */ } else if (y>0) { return ((int64_t) (4*y-2)*y - (x-y) + 1); /* top horizontal */ } else { return ((int64_t) (4*y-2)*y + (x-y) + 1); /* bottom horizontal */ } } void checks (void) { if (option_verbose) printf("checks()\n"); { /* BUFLEN is a power of 2 */ int b = BUFLEN; while (b > 1 && (b&1)==0) b >>= 1; assert (b == 1); } { /* n_to_xy() -> xy_to_n() */ uint64_t n; for (n=1; n <= 63; n++) { int x,y; n_to_xy(n,&x,&y); uint64_t n_again = xy_to_n(x,y); if (0) { DEBUG(printf("n=%"PRIu64" xy=%d,%d back again %"PRIu64"\n", n, x,y, n_again)); } if (n != n_again) error("check_n_to_xy()"); } } { /* xy_to_n() -> n_to_xy() */ int x, y; for (x = -20; x <= 20; x++) { for (y = -20; y <= 20; y++) { uint64_t n = xy_to_n(x,y); int x_again, y_again; n_to_xy(n,&x_again,&y_again); if (! (x==x_again && y==y_again)) error("check_n_to_xy()"); } } } } /* Check that A_pos,B_pos,C_pos are the diagonal back positions shown in the comments above. (This is a check of the initial values, and then of the increments after going around a corner.) */ int check_ABCP_relative_pos (int A_pos, int B_pos, int C_pos, int P_pos) { int x,y; n_to_xy(state.n, &x, &y); int x2,y2; n_to_xy(state.n+1,&x2,&y2); int dx=x2-x, dy=y2-y; int lx = - dy; /* left 90 degrees */ int ly = dx; int ox = -dx + lx; /* diagonally back South-West from dx,dy */ int oy = -dy + ly; DEBUG(printf(" xy=%d,%d dxdy=%d,%d oxoy=%d,%d\n", x,y, dx,dy, ox,oy)); uint64_t A_n = T_pos_to_n(A_pos,P_pos); uint64_t B_n = T_pos_to_n(B_pos,P_pos); uint64_t C_n = T_pos_to_n(C_pos,P_pos); if (! (xy_to_n(x+ ox, y+ oy) == C_n && xy_to_n(x+2*ox, y+2*oy) == B_n && xy_to_n(x+3*ox, y+3*oy) == A_n)) error("oops, A,B,C,P relative positions"); return(1); } void dxdy_turn_left (int *dx, int *dy) { int new_dx = - *dy; int new_dy = *dx; *dx = new_dx; *dy = new_dy; } /* ------------------------------------------------------------------------ */ /* Expected Data, for consistency checks */ static const char want_data[] = { -1, /* no term n=0 */ /* terms to n=10000 */ 1,1,1,1,2,2,3,2,4,3,3,4,2,3,4,4,2,3,3,4,2,3,4,4,2, /* n=1..25 */ 2,2,5,3,2,2,2,4,5,3,2,4,1,1,5,4,1,1,1,5,5,6,1,1,2, /* n=26..50 */ 1,5,4,1,1,3,3,1,5,5,3,2,1,3,3,1,4,5,1,2,4,3,1,2,3, /* n=51..75 */ 5,1,1,3,3,4,1,3,6,4,1,1,3,4,4,1,4,3,1,1,4,2,3,1,4, /* n=76..100 */ 2,2,3,6,2,2,3,3,5,4,3,2,2,6,4,3,2,4,4,3,2,1,2,4,5, /* n=101..125 */ 3,3,2,2,5,4,1,1,2,2,3,1,4,2,6,4,1,2,4,1,2,5,3,2,4, /* n=126..150 */ 3,1,2,4,3,1,2,3,4,6,2,3,4,2,2,5,1,2,3,3,1,4,5,2,2, /* n=151..175 */ 4,3,2,2,1,3,3,2,2,4,3,2,6,4,2,3,1,4,5,3,1,1,3,5,1, /* n=176..200 */ 4,3,1,1,4,5,1,1,2,2,4,3,1,4,5,5,3,1,5,2,3,5,5,2,3, /* n=201..225 */ 3,1,5,6,2,1,1,4,4,5,3,3,4,2,3,3,1,4,4,3,1,4,2,3,3, /* n=226..250 */ 1,2,3,5,4,1,3,4,2,5,3,1,1,5,4,2,1,5,4,2,1,4,3,1,5, /* n=251..275 */ 4,1,1,3,5,1,1,3,4,4,1,2,4,4,5,4,1,2,5,3,1,1,5,4,6, /* n=276..300 */ 5,4,1,2,4,4,1,1,3,4,1,1,3,4,1,3,2,6,5,2,2,3,4,2,2, /* n=301..325 */ 2,3,4,6,2,2,4,3,2,2,4,3,6,5,3,1,4,2,1,3,2,2,1,6,2, /* n=326..350 */ 1,3,2,4,1,1,2,4,1,1,2,2,2,4,5,1,2,6,3,3,2,2,1,1,4, /* n=351..375 */ 2,1,5,4,2,1,3,2,2,1,3,2,5,1,4,5,1,2,4,5,3,1,4,2,5, /* n=376..400 */ 1,1,4,3,1,2,4,3,2,4,3,1,2,3,3,1,2,3,4,5,1,2,2,3,4, /* n=401..425 */ 2,2,3,4,2,2,3,5,5,1,2,3,3,1,1,5,2,3,3,5,2,1,4,4,2, /* n=426..450 */ 2,3,3,1,2,5,3,1,2,5,3,1,2,3,4,5,2,2,4,3,2,2,4,5,3, /* n=451..475 */ 1,1,4,3,1,5,4,1,3,3,1,4,3,1,5,4,1,1,3,4,5,5,1,1,2, /* n=476..500 */ 2,5,3,1,4,3,1,2,5,3,6,4,1,3,3,4,4,5,3,2,4,3,5,2,3, /* n=501..525 */ 3,4,5,1,1,1,3,3,5,1,3,4,2,1,1,4,5,3,3,2,4,3,1,2,4, /* n=526..550 */ 5,1,2,2,3,4,4,1,3,4,4,1,3,2,4,3,1,2,3,1,4,2,3,6,4, /* n=551..575 */ 2,2,4,6,2,2,4,3,1,5,5,6,6,1,2,4,4,5,2,1,4,3,1,2,3, /* n=576..600 */ 3,4,4,1,3,4,1,1,3,4,1,1,3,4,4,1,2,4,4,5,6,6,3,3,2, /* n=601..625 */ 1,2,5,4,1,6,5,2,1,3,5,1,4,4,2,1,3,4,2,1,3,4,4,1,5, /* n=626..650 */ 3,2,5,4,1,1,3,4,1,1,3,4,1,3,2,6,4,2,2,3,4,2,2,1,4, /* n=651..675 */ 2,1,3,2,4,5,2,2,4,3,2,2,4,3,2,2,4,3,3,5,4,1,2,5,5, /* n=676..700 */ 3,1,2,1,4,4,2,1,3,2,5,1,2,2,6,1,2,4,1,1,2,4,1,1,2, /* n=701..725 */ 2,3,3,4,4,3,2,4,4,5,1,2,4,3,5,2,3,5,2,1,1,4,2,1,3, /* n=726..750 */ 4,2,1,3,3,4,1,1,4,3,1,5,2,2,1,3,2,5,1,4,2,1,5,4,2, /* n=751..775 */ 3,1,4,2,3,1,1,5,4,1,1,3,3,5,1,1,3,6,2,2,3,3,4,2,1, /* n=776..800 */ 3,3,2,5,4,3,5,2,1,4,4,2,1,2,2,3,1,2,2,3,4,2,2,3,4, /* n=801..825 */ 2,2,3,1,2,3,3,1,1,2,2,4,4,1,2,4,3,1,4,5,2,1,3,3,2, /* n=826..850 */ 4,3,5,2,2,1,3,5,2,4,3,1,2,5,3,1,2,2,4,3,1,1,2,3,6, /* n=851..875 */ 5,2,2,6,3,2,2,5,5,4,2,3,1,1,4,3,1,1,4,3,3,5,4,2,5, /* n=876..900 */ 3,1,4,5,3,1,1,5,4,1,1,3,4,1,1,3,4,1,1,2,4,5,3,1,2, /* n=901..925 */ 3,1,4,4,3,1,3,3,5,1,2,5,3,1,4,4,3,1,4,3,3,5,2,3,3, /* n=926..950 */ 4,2,6,5,4,4,1,1,2,2,3,1,1,3,2,1,1,3,3,2,1,3,4,1,1, /* n=951..975 */ 6,4,4,2,3,5,4,3,1,2,4,5,1,2,2,4,3,2,4,3,2,2,4,5,1, /* n=976..1000 */ 3,4,4,1,3,2,4,3,6,2,1,5,3,1,4,2,3,5,4,6,2,2,3,3,4, /* n=1001..1025 */ 2,2,6,4,4,7,2,2,4,3,1,1,5,4,2,1,5,4,4,6,1,1,4,2,1, /* n=1026..1050 */ 3,3,2,1,5,4,3,1,4,4,3,5,4,1,3,4,1,1,5,4,1,1,3,4,4, /* n=1051..1075 */ 1,2,4,4,5,5,6,3,1,2,5,3,4,2,1,2,2,4,1,2,5,4,1,2,4, /* n=1076..1100 */ 1,1,3,4,1,3,5,2,1,3,4,2,1,3,4,4,1,5,3,1,2,2,3,2,4, /* n=1101..1125 */ 4,3,2,4,4,1,1,3,4,1,1,3,4,1,3,5,6,2,2,4,3,2,2,4,1, /* n=1126..1150 */ 6,2,1,1,4,3,1,2,2,4,1,2,2,3,5,2,2,4,3,2,2,4,3,2,2, /* n=1151..1175 */ 4,3,3,2,4,5,2,4,5,3,1,2,2,3,4,2,1,4,2,2,1,4,4,2,1, /* n=1176..1200 */ 3,2,5,1,2,2,4,1,1,5,2,1,1,4,2,1,3,3,5,5,4,3,1,4,4, /* n=1201..1225 */ 2,3,5,4,4,3,2,4,4,5,1,2,4,5,2,2,3,3,1,1,2,4,6,1,3, /* n=1226..1250 */ 4,2,1,3,3,4,1,1,3,4,1,1,4,3,1,1,4,3,1,5,2,2,1,3,2, /* n=1251..1275 */ 4,1,1,2,4,4,3,2,4,1,3,2,4,1,1,5,4,1,1,2,4,1,4,3,1, /* n=1276..1300 */ 1,5,3,1,1,5,4,3,5,2,2,3,3,1,2,3,3,6,2,2,3,3,4,2,1, /* n=1301..1325 */ 4,4,2,6,3,1,2,2,2,5,3,1,2,2,3,1,2,2,3,4,2,2,3,4,2, /* n=1326..1350 */ 2,3,1,2,3,3,1,1,2,2,6,3,1,2,4,3,1,5,4,2,1,3,3,4,1, /* n=1351..1375 */ 3,5,2,1,3,3,2,4,3,2,2,4,1,2,4,6,1,3,5,2,1,3,2,2,6, /* n=1376..1400 */ 4,3,1,4,4,3,1,1,2,3,1,1,2,3,5,5,2,2,4,3,2,2,4,5,2, /* n=1401..1425 */ 1,3,4,1,1,3,4,1,5,3,3,2,4,3,5,5,6,2,1,3,1,4,3,3,1, /* n=1426..1450 */ 4,5,3,1,1,5,4,1,1,3,4,1,1,3,4,1,5,2,4,1,1,3,2,6,1, /* n=1451..1475 */ 3,4,5,1,1,4,3,1,3,3,4,1,3,3,5,1,2,5,3,1,4,4,3,1,4, /* n=1476..1500 */ 3,3,5,2,4,3,3,2,4,5,1,2,2,4,3,1,5,6,2,1,3,3,1,4,2, /* n=1501..1525 */ 3,1,1,2,3,1,5,3,2,1,3,3,5,1,4,4,5,2,3,3,1,2,4,3,1, /* n=1526..1550 */ 2,5,4,2,2,3,4,2,2,3,4,4,3,2,2,4,3,2,2,4,5,1,3,4,4, /* n=1551..1575 */ 1,3,2,4,3,3,2,1,1,3,2,4,1,3,2,4,3,2,2,4,3,3,2,5,3, /* n=1576..1600 */ 2,2,6,3,4,2,2,3,4,2,2,6,3,1,1,5,4,4,1,2,4,4,1,1,5, /* n=1601..1625 */ 4,1,1,2,5,3,1,2,5,5,1,3,4,4,1,1,1,3,5,4,1,3,4,4,5, /* n=1626..1650 */ 3,1,5,4,1,1,5,4,1,1,3,4,4,1,2,4,4,5,3,4,4,5,2,1,3, /* n=1651..1675 */ 4,2,1,3,3,2,1,2,4,4,1,2,5,6,1,2,5,4,1,2,4,1,1,3,4, /* n=1676..1700 */ 1,3,3,2,4,3,1,4,4,3,1,4,5,1,3,2,2,1,3,2,2,1,3,2,4, /* n=1701..1725 */ 4,3,2,4,4,3,2,4,4,1,1,3,4,1,1,3,4,1,2,4,3,2,2,4,3, /* n=1726..1750 */ 2,1,4,2,1,1,4,2,1,3,6,4,4,3,1,2,2,4,1,2,2,4,1,2,2, /* n=1751..1775 */ 3,5,4,3,2,2,4,5,2,2,4,6,5,1,4,6,5,5,4,3,1,5,2,2,4, /* n=1776..1800 */ 5,2,3,4,2,2,3,2,4,5,2,1,4,2,2,1,4,4,2,1,3,2,5,1,2, /* n=1801..1825 */ 2,4,1,1,5,2,1,1,4,2,1,3,3,4,5,1,2,6,3,4,2,1,4,4,2, /* n=1826..1850 */ 3,1,4,2,3,5,4,2,3,6,4,1,2,4,4,1,2,2,5,1,3,2,1,4,4, /* n=1851..1875 */ 2,1,3,4,2,1,3,3,4,1,1,3,4,1,1,3,4,1,2,2,3,1,1,4,3, /* n=1876..1900 */ 1,1,4,3,1,5,2,2,1,3,2,4,1,1,2,4,4,3,2,4,1,3,5,5,1, /* n=1901..1925 */ 1,3,4,1,1,2,4,1,1,3,4,1,4,4,5,1,4,3,1,1,4,3,1,1,4, /* n=1926..1950 */ 4,2,2,3,5,1,2,3,5,5,2,2,3,3,4,2,5,4,1,2,4,3,6,2,1, /* n=1951..1975 */ 3,2,5,3,3,2,3,4,6,2,2,1,3,2,2,1,3,2,2,4,3,2,2,4,3 /* n=1976..2000 */ }; /* A366312 */ uint64_t want_first_occurrence[] = { -1, /* no d=0 */ /* 1 2 3 4 5 6 7 8 9 */ /**/ 1, 5, 7, 9, 28, 47, 1031, 745771, 115318369671 }; /* ------------------------------------------------------------------------ */ /* A,B,C to disallow Line of terms ending ...,A,B,C,P where P is a prospective term. Cannot have P=B or that would be palindrome B,C,B. If B=C then cannot have P=A or that would be palindrome A,B,B,A. ABC_to_disallow[A][B][C] = disallow is bits 1<<(B-1) and possible OR of 1<<(A-1) which are values to disallow at P. Values A=0 or B=0 are taken to mean no such term. This happens when P is near a corner of the spiral so there may be only 1 or 2 terms in a given direction relative to P. disallow_to_P[disallow] is the smallest P value which "disallow" bits don't prevent. So the smallest bit 1<<(P-1) which is a 0 in disallow. There's at most 8 1-bits set in disallow so that disallow=0x1FF never occurs. That disallow_to_P[] entry is a dummy P=0; */ /* ABC_to_disallow[] is 9 bits each disallow_to_P[] is only values 0 .. 9 */ static int16_t ABC_to_disallow[10][10][10]; static int8_t disallow_to_P[0x200]; static struct w_offset { int w; int offset; } corners[11][4][3]; void make_tables (void) { if (option_verbose) { printf("make_tables()\n"); printf("sizes in bytes: state %zu\n", sizeof(state)); printf(" ABC_to_disallow %zu, disallow_to_P %zu\n", sizeof(ABC_to_disallow), sizeof(disallow_to_P)); } { /* disallow has bits set 1<<(B-1) bit-OR 1<<(A-1) if/when A or B terms are to be disallowed because they cause a palindrome */ int A,B,C; for (A=0; A<=9; A++) { for (B=0; B<=9; B++) { for (C=0; C<=9; C++) { int disallow = 0; if (B != 0) disallow |= 1<<(B-1); if (A != 0 && B==C) disallow |= 1<<(A-1); ABC_to_disallow[A][B][C] = disallow; assert (disallow < numberof(disallow_to_P)); assert (ABC_to_disallow[A][B][C] == disallow); /* fits data type */ } } } } { /* disallow_to_P[diallow] = smallest P for which bit 1<<(P-1) is a 0-bit in "disallow" */ assert (numberof(disallow_to_P) == 1<<9); int disallow; for (disallow=0; disallow < numberof(disallow_to_P); disallow++) { int P; for (P=1; P<=9; P++) { if ((disallow & (1 << (P-1))) == 0) { disallow_to_P[disallow] = P; break; } } } } { /* corners[] table. dxdy[0] are steps dx,dy in the directions of the 4 lines shown in the comments above, which is for going right-ward along the bottom side of the spiral dxdy[1] are those lines turned +90 degrees, for after turning +90 around the corner. Spiral coordinates x=2,y=-6 are an example of a "0" place shown in the comments above. This is after the side traversal has incremented A_pos,B_pos,C_pos,P_pos. Those pos are held fixed for considering corner lines. Combinations of xy_to_n() and n_to_xy() find cell locations and numbers n and which "w" pos variable to be relative to, and offset from it. */ static const int dxdy[2][4][2] = { { {-1,0}, /* backward (West) */ {-1,1}, /* back diagonal (North-West) */ { 0,1}, /* up perpendicular (North) */ { 1,1} /* forward diagonal (North-East) */ }, { { 0,-1}, /* backward down (South) */ {-1,-1}, /* back diagonal (South-West) */ {-1, 0}, /* left perpendicular (West) */ {-1, 1} /* forward diagonal (North-West) */ } }; int P_x=2, P_y=-6; int A_pos = xy_to_n(P_x-3, P_y+3); int B_pos = xy_to_n(P_x-2, P_y+2); int C_pos = xy_to_n(P_x-1, P_y+1); int P_pos = xy_to_n(P_x, P_y); int ABCP_pos[4]; ABCP_pos[0]=A_pos; ABCP_pos[1]=B_pos; ABCP_pos[2]=C_pos; ABCP_pos[3]=P_pos; DEBUG(printf("corner ABCP_pos %d %d %d %d\n", A_pos, B_pos, C_pos, P_pos)); /* 11 corner cells */ assert (numberof(corners) == 11); int i; for (i=0; i < numberof(corners); i++) { int d = (i <= 5 ? 0 : 1); /* which of dxdy[] */ int x,y; n_to_xy(P_pos+i, &x,&y); /* 4 lines of potential palindrome */ int line; for (line=0; line<4; line++) { int dx = dxdy[d][line][0]; int dy = dxdy[d][line][1]; /* 3 terms in a line, index e=0..2 */ int e; for (e=0; e<3; e++) { int distance = 3 - e; /* e=0 is "A" at distance 3 */ int n = xy_to_n(x + distance*dx, y + distance*dy); DEBUG(printf(" e=%d xy=%d,%d n=%d\n", e,x,y,n)); corners[i][line][e].w = -1; /* for no such element */ /* at or after P_pos+i means no term in the spiral yet */ if (n < P_pos + i) { int w; for (w=3; w>=0; w--) { int T_pos = ABCP_pos[w]; if (n >= T_pos - 12) { corners[i][line][e].w = w; corners[i][line][e].offset = n - T_pos; break; } } } } } } if (option_verbose) { for (i=0; i= 139); assert (BUFLEN >= 139); assert (xy_to_n(3,3) == state.saved_A_pos); assert (xy_to_n(2,4) == state.saved_B_pos); assert (xy_to_n(1,5) == state.saved_C_pos); assert (xy_to_n(0,6) == state.saved_P_pos); /* dummy values marking unused parts of buf[] */ memset (state.buf, -1, sizeof(state.buf)); /* n=1 to n=138 */ for (state.n=1; state.n < state.saved_P_pos; state.n++) { int term = want_data[state.n]; if (! state.first_occurrence[term]) output_first_occurrence (term); if (state.n >= state.saved_A_pos) state.buf[state.n] = term; } assert (state.buf[138] == 4); state.side_length = 1; state.side_increment = 0; } /* Return the "disallow" bit flags for terms located at A_pos, B_pos, C_pos. */ static inline int disallow_at_ABC_pos (int A_pos, int B_pos, int C_pos) { assert (0 <= A_pos && A_pos < BUFLEN); assert (0 <= B_pos && B_pos < BUFLEN); assert (0 <= C_pos && C_pos < BUFLEN); int A = state.buf[A_pos]; int B = state.buf[B_pos]; int C = state.buf[C_pos]; assert (1 <= A && A <= 9); assert (1 <= B && B <= 9); assert (1 <= C && C <= 9); DEBUG(printf(" line %d,%d,%d,p at pos %d,%d,%d disallow=0x%X\n", A,B,C, A_pos,B_pos,C_pos, ABC_to_disallow[A][B][C])); return (ABC_to_disallow[A][B][C]); } /* Get new term from "disallow" bits and store at P_pos in buf[]. state.n is its index n. */ static inline void set_term (int disallow, int P_pos) { assert (0 <= disallow && disallow < numberof(disallow_to_P)); int term = disallow_to_P[disallow]; DEBUG(printf(" disallow 0x%X set %d\n", disallow, term)); assert (1 <= term && term <= 9); assert (state.buf[P_pos] == -1); /* no overwrite */ state.buf[P_pos] = term; assert (term < numberof(state.first_occurrence)); if (UNLIKELY(state.first_occurrence[term] == 0)) output_first_occurrence(term); if (0) { if (UNLIKELY (term == 8)) printf("n=%"PRIu64" term %d\n", state.n, term); } #if WANT_ASSERT if (state.n < numberof(want_data)) { if (term != want_data[state.n]) { error("n=%"PRIu64" want %d got %d\n", state.n, want_data[state.n], term); } } #endif } static inline int w_offset_get_term (const int ABCP_pos[4], const struct w_offset *wp) { assert (-1 <= wp->w && wp->w < 4); if (wp->w < 0) return(0); /* no term */ int T = add(ABCP_pos[wp->w], wp->offset); return (state.buf[T]); } /* Return the "disallow" bit flags for corners[i] line number "line". */ static int corner_disallow_at_line (int i, int line, const int ABCP_pos[4]) { int A = w_offset_get_term(ABCP_pos, &corners[i][line][0]); int B = w_offset_get_term(ABCP_pos, &corners[i][line][1]); int C = w_offset_get_term(ABCP_pos, &corners[i][line][2]); DEBUG(printf(" i=%d line %d %d,%d,%d,p from %d,%d disallow=0x%X\n", i,line, A,B,C, corners[i][line][0].w, corners[i][line][0].offset, ABC_to_disallow[A][B][C])); assert (0 <= A && A <= 9); assert (0 <= B && B <= 9); assert (0 <= C && C <= 9); return (ABC_to_disallow[A][B][C]); } void search (void) { int A_pos = state.saved_A_pos; int B_pos = state.saved_B_pos; int C_pos = state.saved_C_pos; int P_pos = state.saved_P_pos; uint64_t progress_prev = 0; for (;;) { if (WANT_PROGRESS) { uint64_t now = state.n >> 28; if (UNLIKELY(now != progress_prev)) { printf("at n=%"PRIu64"\n", state.n); progress_prev = now; } } /* terms along one side of the square spiral */ int i, this_side = state.side_length; for (i=0; LIKELY(i < this_side); i++) { DEBUG(printf("at n=%"PRIu64" A=%d,B=%d,C=%d,P=%d i=%d of side=%d side_inc=%d\n", state.n, A_pos,B_pos,C_pos,P_pos, i, state.side_length, state.side_increment)); assert (P_pos == (state.n & BUFLEN_MASK)); assert (check_ABCP_relative_pos(A_pos,B_pos,C_pos,P_pos)); int disallow = 0; disallow |= disallow_at_ABC_pos(add(P_pos, -3), /* back */ add(P_pos, -2), add(P_pos, -1)); disallow |= disallow_at_ABC_pos(A_pos, /* back diagonal */ B_pos, C_pos); disallow |= disallow_at_ABC_pos(add(A_pos, 3), /* perpendicular */ add(B_pos, 2), add(C_pos, 1)); disallow |= disallow_at_ABC_pos(add(A_pos, 6), /* forward diagonal */ add(B_pos, 4), add(C_pos, 2)); set_term (disallow, P_pos); PREFETCH_FOR_READ (state.buf + A_pos + 256); PREFETCH_FOR_READ (state.buf + B_pos + 256); PREFETCH_FOR_READ (state.buf + C_pos + 256); PREFETCH_FOR_RW (state.buf + P_pos + 256); assert (zap_buf_entry(A_pos)); /* not needed again */ A_pos = add(A_pos, 1); B_pos = add(B_pos, 1); C_pos = add(C_pos, 1); P_pos = add(P_pos, 1); state.n++; } if (UNLIKELY(add(A_pos, - P_pos) <= 30)) error("buffer size exceeded at n=%"PRIu64"\n", state.n); int ABCP_pos[4]; ABCP_pos[0] = A_pos; ABCP_pos[1] = B_pos; ABCP_pos[2] = C_pos; ABCP_pos[3] = P_pos; DEBUG(printf("corners ABCP_pos %d %d %d %d\n", A_pos, B_pos, C_pos, P_pos)); for (i=0; i