//==========================================================================
// pick.cc -- For which triangles is Pick's formula almost correct?
// Copyright (C) 2006 Michael Eisermann, all rights reserved
// Institut Fourier, Université Joseph Fourier, Grenoble, France
//
// Description:
// This little program verifies a conjecture of Casson and Gordon. 
// For details see http://www.arxiv.org/abs/math.GT/0602395.html
//
// Algorithm:
// The idea is to carry out a brute force count of lattice points.
// Within this approach the following program is reasonably well
// optimized.  We point out, however, that there may be different 
// and more efficient ways to obtain such lattice point sums.
//
// Complexity:
// Memory cost is fixed and negligible.  Time cost is O(p^4) for fixed p:
// we have to check q=2,...,p^2-1 and for each q a loop of length p^2.
// 
// Compilation:
// g++ pick.cc -o pick -O4 -W -Wall 
//==========================================================================
#include <iostream>
using namespace std;

// We assume to work with 64bit integers.  This allows at least
// for a range p = 3,...,59999 which seems more than sufficient.
typedef long int Integer;

// Notice that 32bit integers would only suffice for p=3,...,249.
// On a 32bit processor you might thus try the following:
// typedef long long int Integer;

//--------------------------------------------------------------------------
// Calculate the gcd of two integers p and q 
//--------------------------------------------------------------------------
Integer gcd( Integer p, Integer q )
{
  while( q != 0 ) { p%= q; swap(p,q); };
  return abs(p);
}

//--------------------------------------------------------------------------
// This functions checks whether Pick's formula is almost correct for 
// the triangle (p,q) and magnifications by a factor r=1,...,p-1.
//--------------------------------------------------------------------------
Integer correct( Integer p, Integer q )
{
  // We only consider coprime positive integers p,q with 
  // 1 < q < p^2, p odd, q even.  All others are rejected.
  // if( p<1 || q<2 || q>=p*p || q%2!=0 || gcd(p,q)!=1 ) return false;

  // Interior points count 2, boundary points count 1
  Integer count= 0;   // count lattice points
  Integer pp= p*p;    // store p^2 once and for all
  Integer rest= q;    // add q for each new column
  Integer height= 0;  // current position is  height/2 + rest/p^2

  // Loop through all magnification factors r=1,...,p-1
  for( Integer r=1; r<p; ++r )
    {
      // Count columns  k = 1,...,p-1
      for( Integer k= 1; k<p; ++k, rest+= q )
	{
	  // We carry out a division with remainder; since this is done 
	  // in each column we can replace this by iterated subtraction.
	  if( rest >= pp ) 
	    { rest-= pp; height+= 2; }

	  // Add interior points
	  count+= height;
	  
	  // Add the boundary point at the bottom.  This is superfluous 
	  // whenever  rest==0  because then the additional boundary 
	  // point at the top has already been counted twice.
	  if( rest ) ++count;
	}
      
      // Count the last column  k = p as boundary points
      if( rest >= pp ) { rest-= pp; height+= 2; }
      count+= 1 + height/2;
      
      // Check difference: it has to be 0 or 1.  If not then we can
      // short-circuit this case and abandon further counting.
      Integer diff= count - r*r*q;
      if( diff != 0 && diff != 1 ) return false;
      
      // Continue counting
      count+= height/2;
      rest+= q;
    }
  return true;
}

//--------------------------------------------------------------------------
// Test whether (p,q) belongs to one of the four Casson-Gordon families.
// For details see  http://www.arxiv.org/abs/math.GT/0602395.html
//--------------------------------------------------------------------------
bool family( Integer p, Integer q )
{
  // First family
  if( (q-1)%p == 0 )
    {
      Integer n= (q-1)/p;
      if( gcd(n,p) == 1 ) return true;
    }
  if( (q+1)%p == 0 )
    {
      Integer n= (q+1)/p;
      if( gcd(n,p) == 1 ) return true;
    }

  // Second and third family combined
  if( q%(p+1) == 0 )
    { 
      Integer n= q/(p+1);
      if( (2*p-1)%n == 0 ) return true;
      if( (p+1)%n == 0 && n%2 != 0 ) return true; 
    }
  if( q%(p-1) == 0 )
    { 
      Integer n= q/(p-1);
      if( (2*p+1)%n == 0 ) return true;
      if( (p-1)%n == 0 && n%2 != 0 ) return true; 
    }
  
  // Fourth family
  if( q%(2*p+1) == 0 )
    { 
      Integer n= q/(2*p+1);
      if( (p-1)%n == 0 )
	if( ((p-1)/n)%2 != 0 ) return true;
    }
  if( q%(2*p-1) == 0 )
    { 
      Integer n= q/(2*p-1);
      if( (p+1)%n == 0 )
	if( ((p+1)/n)%2 != 0 ) return true;
    }
  return false;
}      

//--------------------------------------------------------------------------
// The function main() searches the interval  pstart,...,pstop 
//--------------------------------------------------------------------------
int main()
{
  // Welcome message
  cout << endl
       << "pick -- For which triangles is Pick's formula almost correct?\n"
       << "Copyright (C) 2006 Michael Eisermann, all rights reserved\n"
       << "Institut Fourier, Université Joseph Fourier, Grenoble, France\n"
       << endl
       << "This little program verifies a conjecture of Casson and Gordon.\n" 
       << "For details see http://www.arxiv.org/abs/math.GT/0602395.html\n" 
       << endl;

  // Security check
  if( sizeof(Integer) < 8 ) 
    {
      cout << "*** Warning ***\n"
	   << "It seems that the current compilation is not working with\n"
	   << "64bit integers.  This will be insufficient and will lead\n"
	   << "to erroneous results.  Please check compiler options and\n"
	   << "the source code, and then recompile for 64bit integers.\n"
	   << endl;
      exit(1);
    }

  // Enter the interval to search
  cout << "Please enter the interval for p : ";
  Integer pstart, pstop;
  cin >> pstart >> pstop;
  if ( pstart<3 ) pstart= 3;
  if ( pstart%2==0 ) ++pstart;
  if ( pstop>59999 ) pstop= 59999;
  if ( pstop%2==0 ) --pstop;
  if ( pstart > pstop ) swap( pstart, pstop );
  cout << "Searching for solutions with p=" << pstart << ",...," << pstop << endl;

  // The main loop verifies if all solutions with p=pstart,...,pstop
  // belong to one of the Casson-Gordon families, as conjectured.
  for( Integer p=pstart; p<=pstop; p+=2 )
    {
      cout << endl << "p=" << p << ": ";
      for( Integer q=2; q<p*p; q+=2 )
	{
	  // We only consider coprime pairs (p,q), all others are rejected.
	  if( gcd(p,q)!=1 ) continue;
	  
	  // If (p,q) belongs to one of the Casson-Gordon families, then 
	  // we already know that it almost satisfies Pick's formula.
	  // There is thus no point in recouting.
	  if( family(p,q) ) continue;

	  // If we arrive here and (p,q) almost satisifes Pick's formula
	  // then it is a counter-example to the conjecture.
	  if( correct(p,q) )
	    {
	      // Conclusion in the case of a counter-example.
	      cout << endl << "The conjecture is false!" << endl
		   << "Exceptional pair (" << p << "," << q << ")" << endl;
	      exit(0);
	    }

	  // Notice that the order of the different conditions is not 
	  // arbitrary.  The function correct(p,q) is most expensive 
	  // for pairs that almost satisfies Pick's formula, because
	  // these cannot be short-circuited and the count has to go 
	  // all the way through all loops.  It is thus wise to remove
	  // the Casson-Gordon families *before* calling correct(p,q).
	}
      cout << "ok";
    }

  // Conclusion at the end of the program
  cout << endl << "I have finished searching the interval p=" << pstart << ",...," << pstop
       << endl << "without finding any counter-example.\n" << endl;
}