/*
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;
}