/* Explorer for the graph produced by the exchange relation on Markoff-
   type triples on the field modulo p
   Version of Dec. 4 2003
   Copyright 2003 Martin Hock (mhock@cs.wisc.edu)
   A work in progress - provided without warranty.
   Permission granted to read and modify.
   Permission required to redistribute.
*/
/* #define NDEBUG */
/* #define DEBUG */
/* #define SORT */
#define NODECONJECTURE
/* #define FINDCONNECTED */
#define FINDTRIANGLE
/* #define OUTPUTGRAPH */
/* #define SACOUNT */
/* #define NNCOUNT */
#define W 1L
#define X 2L
#define Y 3L
#define Z 6L

#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
long primes[] = 
  {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 
   59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 
   139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 
   227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 
   311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 
   401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 
   491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 
   599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 
   683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 
   797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 
   887, 907, 911, 919, 929, 937, 941, 947};

#define NUMPRIMES (sizeof(primes)/sizeof(long))

typedef long node[3];

inline node *findNode(node *nodes, int *marks, long numNodes, 
		      long a, long b, long c, int doMark) {
  long i = 0, j = numNodes-1;
#ifdef DEBUG
  printf("%li %li %li\n",a,b,c);
#endif
  while (i <= j) {
    long k = (i+j)/2;
#ifdef DEBUG
    printf("%li %li %li (%li %li %li)\n",i,j,k,
	   (nodes[k][0]),(nodes[k][1]),(nodes[k][2]));
#endif

    if ((nodes[k][0]) == a && 
	(nodes[k][1]) == b && 
	(nodes[k][2]) == c) {
      if (marks[k]) return NULL;
      else {
	if (doMark) marks[k] = 1; 
	return nodes+k;
      }
    }
    if ((nodes[k][0]) != a) {
      if ((nodes[k][0]) > a) j = k - 1;
      else i = k + 1;
    } else if ((nodes[k][1]) != b) {
      if ((nodes[k][1]) > b) j = k - 1;
      else i = k + 1;
    } else {
      if ((nodes[k][2]) > c) j = k - 1;
      else i = k + 1;
    }
  }
  assert(0);
#ifdef NDEBUG
  printf("Unfound triple %li %li %li, recompile in debug mode\n",
	 a,b,c);
#endif
  exit(0);
  return NULL;
}

inline int isAdjacent(node *n1, node *n2, long base) {
  long a;
  long b;
  long c;
  a = n1[0][0];
  b = n1[0][1];
  c = n1[0][2];
  /* Can't be adjacent if too different */
  if ((a == n2[0][0]) + (b == n2[0][1]) + (c == n2[0][2]) < 2) return 0;
  if (a != n2[0][0]) {
    if (((((Z/W)*b % base) *c % base) + (base-a)) % base == n2[0][0]) return 1;
  } else if (b != n2[0][1]) {
    if (((((Z/X)*a % base) *c % base) + (base-b)) % base == n2[0][1]) return 1;
  } else { /* c != n2[0][2] */ 
    if (((((Z/Y)*a % base) *b % base) + (base-c)) % base == n2[0][2]) return 1;
  }
  return 0;
}

inline int fillQueue(node **queue, int i, int j, 
		     node *nodes, int *marks, long numNodes, long base,
		     int doMark) {
  long a;
  long b;
  long c;
  long x;
  node *idx;
  a = ((*(queue[i]))[0]);
  b = ((*(queue[i]))[1]);
  c = ((*(queue[i]))[2]);
#ifdef DEBUG
  printf("Examining node (%li %li %li)\n",
	 ((*(queue[i]))[0]),
	 ((*(queue[i]))[1]),
	 ((*(queue[i]))[2]));
#endif
      
#ifdef DEBUG
  printf("Swap a: ");
#endif
  x = ((((Z/W)*b % base) *c % base) + (base-a)) % base;
#ifdef SORT      
  if (x < c) {
    idx = findNode(nodes, marks, numNodes, b, c, x, doMark); 
  } else if (x < b) {
    idx = findNode(nodes, marks, numNodes, b, x, c, doMark);
  } else     
#endif
    idx = findNode(nodes, marks, numNodes, x, b, c, doMark);

  if (idx != NULL) queue[j++] = idx;
      
#ifdef DEBUG
  printf("Swap b: ");
#endif
  x = ((((Z/X)*a % base) *c % base) + (base-b)) % base;
#ifdef SORT
  if (a < x) {
    idx = findNode(nodes, marks, numNodes, x, a, c, doMark);
  } else if (x < c) {
    idx = findNode(nodes, marks, numNodes, a, c, x, doMark);
  } else
#endif
    idx = findNode(nodes, marks, numNodes, a, x, c, doMark);

  if (idx != NULL) queue[j++] = idx;
      
#ifdef DEBUG
  printf("Swap c: ");
#endif
  x = ((((Z/Y)*a % base) *b % base) + (base-c)) % base;
      
#ifdef SORT
  if (a < x) {
    idx = findNode(nodes, marks, numNodes, x, a, b, doMark);
  } else if (b < x) {
    idx = findNode(nodes, marks, numNodes, a, x, b, doMark);
  } else
#endif
    idx = findNode(nodes, marks, numNodes, a, b, x, doMark);

  if (idx != NULL) queue[j++] = idx;
  return j;
}

#ifdef SORT
#define MAXB a
#define MAXC b
#else
#define MAXB (base-1)
#define MAXC (base-1)
#endif

int main (int argc, char *argv[]) {
  long base; /* We are working in the field Z mod base (base must be prime) */
  node *nodes; /* A sorted array of Markoff triples mod base */
  int *marks; /* A boolean array to mark whether we have visited triples */
  node **queue; /* A moving queue of triples to visit */
  long a;
  long i;
  long p;
  long numNodes;
#ifdef SACOUNT
  long selfAdj;
#endif
#ifdef SORT
  if (W != X || X != Y || Y != W) {
    printf("W, X, Y not all equal: Recompile without #define SORT\n");
    return 0;
  }
#endif
  if (argc != 2) {
    printf("Usage: %s <prime>\n", argv[0]);
    return 0;
  }
  base = atol(argv[1]);
  if (base <= primes[NUMPRIMES] && base != 0) {
    int i = 0, j = NUMPRIMES-1;
    while (i <= j) {
      int k = primes[(i+j)/2];
      if (k == base) break;
      else if (k < base) 
	/* We know that k is in the range [(i+j)/2+1,j] */
	i = (i+j)/2+1;
      else // k > base
	/* We know that k is in the range [i,(i+j)/2-1] */
	j = (i+j)/2-1;
    }
    if (i > j) {
      printf("Given base is not a prime\n");
      return 0;
    }
  }
  for (p = 0; p < NUMPRIMES; p++) {
    if (p == 0 && base != 0) p = sizeof(primes)/sizeof(long);
    else base = primes[p];
    i = 0; numNodes = 0;
#ifdef SACOUNT
    selfAdj = 0;
#endif
#ifdef NODECONJECTURE
    /* The conjecture is that numNodes <= base*base + 3*base */
    if (!(W == 1 && X == 1 && Y == 1 && Z == 3) &&
	!(W == 1 && X == 1 && Y == 2 && Z == 4) &&
	!(W == 1 && X == 2 && Y == 3 && Z == 6)) {
      printf("Conjecture untested for %lia^2 + %lib^2 + %lic^2 == %liabc\n",
	     W,X,Y,Z);
      return 0;
    }
    numNodes = base*base+3*base;
#else
    for (a = 0; a < base; a++) {
      long b = 0;
      for (; b <= MAXB; b++) {
	long c = 0;
	for (; c <= MAXC; c++) {
	  if (a+b+c == 0) continue;
	  if (((W*a % base)*a % base + 
	       (X*b % base)*b % base + 
	       (Y*c % base)*c % base) % base == 
	      (((Z*a % base)*b % base)*c) % base) 
	    numNodes++;
	}
      }
    }
#endif
    nodes = malloc(sizeof(node)*numNodes);
    queue = malloc(sizeof(node *)*(numNodes+1)); /* +1 for girth */
    marks = calloc(numNodes,sizeof(int));
    if (!nodes || !queue || !marks) return 0;
    for (a = 0; a < base; a++) {
      long b = 0;
      for (; b <= MAXB; b++) {
	long c = 0;
	for (; c <= MAXC; c++) {
	  if (a+b+c == 0) continue;
	  if (((W*a % base)*a % base + 
	       (X*b % base)*b % base + 
	       (Y*c % base)*c % base) % base == 
	      (((Z*a % base)*b % base)*c) % base) {
#ifdef DEBUG
	    printf("%li %li %li\n",a,b,c);
#endif
	    nodes[i][0] = a;
	    nodes[i][1] = b;
	    nodes[i][2] = c;
	    i++;
#ifdef SACOUNT
	    if ((((((Z/W)*b % base) *c % base) + (base-a)) % base == a) +
		(((((Z/X)*a % base) *c % base) + (base-b)) % base == b) +
		(((((Z/Y)*a % base) *b % base) + (base-c)) % base == c) > 1)
	      selfAdj++;
#endif /*SACOUNT*/
	  }
	}
      }
    }
#ifdef NODECONJECTURE
    assert (i <= numNodes);
    if (i > numNodes) {
      printf("Conjecture is false for base %li\n",base);
      return 0;
    }
    numNodes = i;
#else
    assert (i == numNodes);
#endif
#ifdef SACOUNT
    printf("%li %li\n",base,selfAdj);
    fflush(stdout);
#endif
#ifdef NNCOUNT
    printf("%li %li\n",base,numNodes);
    fflush(stdout);
#endif
#ifdef FINDCONNECTED
#ifdef DEBUG
    printf("-CONN-\n");
#endif
    // Mark node 0 as seen
    marks[0] = 1;
    for (i = 0, j = 1, queue[0] = nodes+(0); i < j; i++) {
      j = fillQueue(queue, i, j, nodes, marks, numNodes, base, 1);
      assert(j <= numNodes);
    }
    if (j != numNodes) {
      printf("Base %li disconnected (reached %li of %li nodes)\n",
	     base,j,numNodes);
      for (i = 0; i < numNodes; i++)
	if (!marks[i]) 
	  printf("(%li, %li, %li)\n",nodes[i][0],nodes[i][1],nodes[i][2]);
    }
#endif /* FINDCONNECTED */
#ifdef FINDTRIANGLE
    /* This code does not work at the moment */
#ifdef DEBUG
    printf("-TRIANGLE-\n");
#endif    
    /* For each node, find a triangle */
    for (a = 0; a < numNodes; a++) {
#ifdef DEBUG
      printf("Searching from node %li\n",a);
#endif
      marks[a] = 1;
      queue[0] = nodes+(a);
      i = fillQueue(queue, 0, 1, nodes, marks, numNodes, base, 0);
      if (i >= 3) { /* Check for adjacency between (first two) nodes */
	if (isAdjacent(queue[1],queue[2], base)) break;
      }
      if (i >= 4) {
	if (isAdjacent(queue[2],queue[3], base) || 
	    isAdjacent(queue[3],queue[1], base))
	  break;
      }
    }
    if (a < numNodes) {
      printf("Triangle found\n");
    } else {
      printf("Triangle not found\n");
    }
#endif /* FINDTRIANGLE */
    free(nodes);
    free(queue);
    free(marks);
#ifdef DEBUG
    printf("\n");
#endif
  }
  return 0;
}