/* A226595 Lengths of maximal nontouching increasing paths in n X n grids. Search algorithm by Bert Dobbelaere Results: N a(N) Longest path (not unique!) 1 0 Trivial (not supported) 2 2 [(0,1),(1,1),(0,0)] 3 4 [(0,2),(1,2),(0,1),(2,1),(0,0)] 4 7 [(2,0),(2,1),(3,0),(3,2),(1,1),(3,3),(0,3),(1,0)] 5 9 [(0,4),(1,4),(2,3),(0,3),(2,2),(4,4),(4,1),(1,2),(4,0),(0,0)] 6 12 [(3,0),(2,0),(3,1),(1,1),(3,2),(5,0),(5,3),(2,2),(5,4),(1,4),(5,5),(0,5),(1,0)] 7 15 [(2,0),(3,0),(4,1),(2,1),(4,2),(6,0),(6,3),(3,2),(0,0),(0,4),(4,5),(1,2),(5,5),(0,6),(6,5),(1,1)] 8 17 [(4,6),(3,6),(4,7),(2,7),(3,5),(1,7),(1,4),(0,1),(3,3),(6,6),(2,4),(6,7),(7,2),(2,0),(6,4),(1,0),(6,5),(0,0)] 9 20 [(3,8),(3,7),(2,8),(0,8),(2,7),(0,5),(0,2),(1,5),(1,1),(2,5),(5,8),(3,4),(6,8),(4,3),(8,7),(8,1),(2,0),(7,4),(1,0),(7,5),(0,0)] 10 24 [(5,6),(4,6),(5,5),(5,3),(4,5),(2,7),(5,7),(2,8),(0,5),(0,1),(1,5),(4,2),(2,6),(5,2),(6,7),(1,9),(7,9),(6,3),(8,9),(8,2),(9,9),(9,1),(1,2),(9,0),(0,0)] 11 27 [(6,3),(5,3),(6,4),(6,6),(5,4),(3,2),(6,2),(3,1),(1,4),(1,8),(2,4),(5,7),(3,3),(6,7),(7,2),(2,0),(8,0),(7,6),(9,0),(9,7),(10,0),(10,8),(2,7),(10,9),(1,9),(10,10),(0,10),(1,0)] 12 29 [(3,1),(4,1),(5,2),(5,0),(6,2),(8,0),(11,0),(8,1),(11,3),(7,3),(3,2),(6,5),(10,7),(6,4),(11,5),(9,10),(5,6),(8,11),(2,11),(3,5),(7,10),(3,4),(1,11),(1,3),(0,11),(1,2),(8,8),(1,1),(9,8),(0,0)] 13 33 [(8,2),(7,2),(6,3),(6,5),(5,3),(7,1),(4,1),(5,4),(7,7),(7,3),(8,7),(11,10),(9,6),(12,10),(11,5),(9,0),(3,0),(9,1),(11,7),(8,1),(1,0),(5,6),(0,0),(0,8),(7,12),(1,6),(8,12),(1,5),(9,12),(1,4),(10,12),(1,3),(11,12),(2,1)] 14 36 [(4,11),(5,11),(6,10),(6,8),(7,10),(5,12),(8,12),(7,9),(5,6),(5,10),(4,6),(1,3),(3,7),(0,3),(1,8),(3,13),(9,13),(3,12),(1,6),(4,12),(11,13),(7,7),(13,12),(13,4),(6,0),(12,6),(5,0),(12,7),(4,0),(12,8),(3,0),(12,9),(2,0),(12,10),(1,0),(12,11),(0,0)] 15 39 [(10,0),(11,0),(12,1),(10,1),(8,0),(10,2),(13,2),(14,5),(11,3),(7,3),(11,4),(14,7),(10,5),(13,9),(14,14),(12,9),(9,14),(3,14),(9,13),(11,7),(8,13),(1,13),(8,12),(2,8),(9,10),(2,7),(10,7),(2,6),(0,14),(0,5),(9,6),(0,4),(10,6),(0,3),(11,5),(0,2),(12,3),(0,1),(13,3),(0,0)] Benchmark: Single thread execution on Intel Core i7-4770K, 64 bit Linux, gcc 5.4.0 "-O4". 36" for N=12 7'58" for N=13 70'26" for N=14 1001' for N=15 */ #define N 12 ///< Grid dimension, code below computes a226595(N) #define MAXDEPTH 100 ///< Limit for trace length - good enough. #include <cstdlib> #include <cstdio> #include <cassert> #include <vector> #include <cstring> #if __cplusplus < 201103L typedef unsigned long long u64; typedef unsigned u32; typedef int s32; typedef unsigned short u16; #else #include <cstdint> typedef uint64_t u64; typedef uint32_t u32; typedef int32_t s32; typedef uint16_t u16; #endif typedef std::vector<u32> LengthList; /* How many 64 bit words to allocate in class below */ #define NMASKWORDS ((N*N*N*N+63)/64) /** * Option keeper class * For each potential source and destination point, keep track of whether the transition is allowed. * Options are stored in densely packed array to reduce memory bandwidth and increase cache performance. */ class Opt_T { public: Opt_T(); Opt_T(const Opt_T &); void markbad(u32 x1, u32 y1, u32 x2, u32 y2) {u32 bitno=N*(N*((N*x1)+y1)+x2)+y2;o[bitno>>6]&=~(1ULL<<(bitno&0x3F));} void filter(const Opt_T &f); bool canGoBetween(u32 p1, u32 p2) const; private: u64 o[NMASKWORDS]; }; /** * Consult the option * p1 and p2 are source and destination candidates (p=N*x+y encoded) */ bool Opt_T::canGoBetween(u32 p1, u32 p2) const { u32 bitno=(N*N)*p1+p2; return ((o[bitno>>6]>>(bitno&0x3F))&1) != 0; } /** * Make a copy */ Opt_T::Opt_T(const Opt_T &f) { memcpy(o,f.o,sizeof(o)); } /** * Filter out invalid options, bit parallel implementation. */ void Opt_T::filter(const Opt_T &f) { for(u32 k=0;k<NMASKWORDS;k++) o[k]&=f.o[k]; } /** * New option set, everything allowed */ Opt_T::Opt_T() { memset(o,0xFF,sizeof(o)); } /** * Large array of filter pointers. * For each transition (from 1st to 2nd index), keeps a filter * to be applied to all future potential transitions. */ Opt_T *excludeMasks[N*N][N*N]; /** * Determinant, to check whether points are collinear or switched sides. */ s32 det(s32 x1, s32 y1, s32 x2, s32 y2, s32 x3, s32 y3) { x2-=x1;x3-=x1; y2-=y1;y3-=y1; return (x2*y3)-(x3*y2); } /** * Determine if a transition from (x1,y1) to (x2,y2) can be followed by a transition from (x3,y3) to (x4,y4) */ bool canhappen(s32 x1, s32 y1, s32 x2, s32 y2, s32 x3, s32 y3, s32 x4, s32 y4) { s32 d12=(x2-x1)*(x2-x1) + (y2-y1)*(y2-y1); s32 d34=(x4-x3)*(x4-x3) + (y4-y3)*(y4-y3); if(d12==0) return false; // No displacement if(d34==0) return false; // No displacement if(d34>=d12) return false; // Increasing vs diminishing distance. Algorithm used decreasing lengths and reverses the order at the end. s32 dt123=det(x1,y1,x2,y2,x3,y3); s32 dt124=det(x1,y1,x2,y2,x4,y4); s32 dt134=det(x1,y1,x3,y3,x4,y4); s32 dt234=det(x2,y2,x3,y3,x4,y4); if((dt123*dt124)>0) return true; // No touching: OK if((dt134*dt234)>0) return true; // No touching: OK if((x2==x3)&&(y2==y3)) // 3->4 as continuation of 1->2 { if(dt124!=0) return true; // not collinear: OK s32 d14=(x4-x1)*(x4-x1) + (y4-y1)*(y4-y1); if(d14>d34) return true; // 180deg angle: OK (?) } return false; // Collision! } u32 maxlength=0; ///< Longest path found so far LengthList ll; ///< List of increasing square lengths (all options in N*N grid) u32 trace[MAXDEPTH]; ///< Keeps track of visited points in recursive search. /** * Display a trace. Reverses order to get again an increasing length sequence. */ void dumptrace(u32 depth) { printf("N=%u segs=%u [",N,depth); for(u32 i=0;i<=depth;i++) { printf("(%d,%d)",trace[depth-i]/N,trace[depth-i]%N); if(i<depth)printf(","); } printf("]\n"); } /** * "Distance class" to get from one position to the other * Both indices are positions encoded as p=N*x+y * The distance class is an integer between 0 and OEIS a047800(N-1)-1 * We don't need to measure distances, we only need to be able to compare them. */ u16 distclass[N*N][N*N]; /** * For each starting position, keep a list of target positions sorted by increasing distance (nextpdc). * "p" component is the destination point (encoded p=N*x+y), "dc" its distance class from the starting point * startidx tracks the start index in the nextpdc array for each distance class. */ struct NBD{ struct {u16 p;u16 dc;}nextpdc[N*N]; // sorted by increasing distance. p= destination (N*x+y), dc="distance class" [ 0 ... a047800(N-1)-1 ] u32 startidx[N*N]; // nextp start index for distance class N } nodesbydistance[N*N]; /** * Initialize list of increasing square lengths */ void composeLengthList() { ll.push_back(0); bool found=true; while(found) { found=false; u32 minval=ll[ll.size()-1]; u32 maxval=100000000; for(u32 x=0;x<N;x++) for(u32 y=0;y<N;y++) { u32 s=x*x+y*y; if(s>minval) { found=true; if(s<maxval) maxval=s; } } if(found) { ll.push_back(maxval); } } printf("ll.size()=%lu\n",ll.size()); // Should be in OEIS A047800. } /** * Initialize distance tables. */ void init_distances() { composeLengthList(); for(s32 x1=0;x1<N;x1++) for(s32 y1=0;y1<N;y1++) for(s32 x2=0;x2<N;x2++) for(s32 y2=0;y2<N;y2++) { u32 dist=(x2-x1)*(x2-x1)+(y2-y1)*(y2-y1); u32 idx=0; while(ll[idx]<dist) idx++; u32 p1=N*x1+y1; u32 p2=N*x2+y2; distclass[p1][p2]=idx; } for(u32 p1=0;p1<(N*N);p1++) { u32 writeidx=0; for(u32 d=0;d<ll.size();d++) { nodesbydistance[p1].startidx[d]=writeidx; for(u32 p2=0;p2<(N*N);p2++) { if(distclass[p1][p2]==d) { nodesbydistance[p1].nextpdc[writeidx].p=p2; nodesbydistance[p1].nextpdc[writeidx].dc=d; writeidx++; } } } nodesbydistance[p1].startidx[ll.size()]=writeidx; assert(writeidx==(N*N)); } } /** * The core operation: recursive search. * @param state Cumulative filter containing all currently allowed options * @param depth Number of line segments (in decreasing length order) already followed in current path. Positions kept in "trace" array * @param last_dc "distance class" of last transition */ void search(const Opt_T &state, u32 depth, u16 last_dc) { if(depth > maxlength) // We improved the maximum path length found so far. Deserves to be printed. { dumptrace(depth); maxlength=depth; } // If we want to have any hope of finding a path that is 1 longer than found so far, we need to be able to reduce the distances this many times, // so we need to start high enough u32 mindistclass=maxlength+1-depth; // This is where we are now. (p = N*x+y encoded) u32 p1=trace[depth]; // We start at the minimum distance that has a change to find a longer path. u32 minidx=nodesbydistance[p1].startidx[mindistclass]; // We stop when we encounter the same distance as last transition u32 maxidx=nodesbydistance[p1].startidx[last_dc]; for(u32 k=minidx;k<maxidx;k++) { u16 p2=nodesbydistance[p1].nextpdc[k].p; u16 dc=nodesbydistance[p1].nextpdc[k].dc; if(state.canGoBetween(p1,p2)) // Check legality of transition before going one level deeper { Opt_T newstate=state; newstate.filter(*excludeMasks[p1][p2]); trace[depth+1]=p2; search(newstate,depth+1, dc); } } } int main() { // Start with the distance tables init_distances(); // Initialize transition filters for(u32 x1=0;x1<N;x1++) for(u32 y1=0;y1<N;y1++) for(u32 x2=0;x2<N;x2++) for(u32 y2=0;y2<N;y2++) { excludeMasks[N*x1+y1][N*x2+y2] = new Opt_T(); for(u32 x3=0;x3<N;x3++) for(u32 y3=0;y3<N;y3++) for(u32 x4=0;x4<N;x4++) for(u32 y4=0;y4<N;y4++) if(!canhappen(x1,y1,x2,y2,x3,y3,x4,y4)) { excludeMasks[N*x1+y1][N*x2+y2]->markbad(x3,y3,x4,y4); } } // Fill the first segment of the path (2 points), apply the first filter and start the recursive search. // Symmetry used: chose starting point in 1st quadrant, on or below y=x // If starting point is on the diagonal, the 2nd point should be on or below it. for(u32 x1=0;x1<((N+1)/2);x1++) for(u32 y1=0;y1<=x1;y1++) for(u32 x2=0;x2<N;x2++) for(u32 y2=0;y2<N;y2++) if((y1<x1) || (y2<=x2)) if((x1!=x2) || (y1!=y2)) { Opt_T state; state.filter( *excludeMasks[N*x1+y1][N*x2+y2]); trace[0]=N*x1+y1; trace[1]=N*x2+y2; search(state,1, distclass[N*x1+y1][N*x2+y2]); } printf("\na(%u) = %u\n\n",N,maxlength); return 0; }