// This search routine searches for stable billiard paths
// whose word in the group of reflections is a square 
// of an odd word

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#define Pi 3.141592653589
#include "plist.c"

typedef struct
{ word w;            // Stores the word
  int count[3];      // signed sum of occurances of each letter
  int mod2;          //
  matrix m; 
  double interval[2];
  Pletter lexi;      // Pointer for the lexitest
  int lexi_count;
} oss_word;

void oss_word_destroy(w)
     oss_word *w;
{
  w->lexi=NULL;
  word_destroy(&(w->w));
}

void oss_word_empty(w)
     oss_word *w;
{
  w->count[0]=w->count[1]=w->count[2]=0;
  w->w.first=w->w.last=NULL;
  w->w.l=0;
  w->mod2=1;
  w->m.m[0][0]=1; w->m.m[0][1]=0; w->m.m[0][2]=0;
  w->m.m[1][0]=0; w->m.m[1][1]=1; w->m.m[1][2]=0;
  w->interval[0]=0;
  w->interval[1]=1;
}

void oss_addLetter(w,i)
     oss_word *w;
     int i;
{
  if (w->w.last!=NULL){   // theres is at least one letter
    w->w.last->next=(Pletter) malloc (sizeof(struct pletter));
    w->w.last->next->l=i;
    w->w.last->next->next=NULL;
    w->w.last->next->prev=w->w.last;
    w->w.last=w->w.last->next;
  } else { // the first letter
    w->w.last=(Pletter) malloc (sizeof(struct pletter));
    w->w.last->l=i;
    w->w.first=w->w.last;
    w->w.last->next=NULL;
    w->lexi=NULL;
    w->w.last->prev=NULL;
  }
  w->count[i]+=w->mod2;
  w->mod2*=-1;
  w->w.l++;
}

void oss_word_copy(w,new)
     oss_word *w,*new;
{
  int i;
  Pletter lw;

  // This can no longer be used, because we need to keep track of
  // lexi information
  //  word_copy(&(w->w),&(new->w));


  new->w.l=0;
  new->w.first=new->w.last=NULL;
  i=0;
  lw=w->w.first;
  while (lw!=NULL){
    oss_addLetter(new,lw->l);
    if (lw==w->lexi)
      new->lexi=new->w.last;
    lw=lw->next;
  }
  if (lw==w->lexi)
    new->lexi=new->w.last;
  new->lexi_count=w->lexi_count;

  new->count[0]=w->count[0];
  new->count[1]=w->count[1];
  new->count[2]=w->count[2];
  new->mod2=w->mod2;
  new->m.m[0][0]=w->m.m[0][0];
  new->m.m[0][1]=w->m.m[0][1];
  new->m.m[0][2]=w->m.m[0][2];
  new->m.m[1][0]=w->m.m[1][0];
  new->m.m[1][1]=w->m.m[1][1];
  new->m.m[1][2]=w->m.m[1][2];
  new->interval[0]=w->interval[0];
  new->interval[1]=w->interval[1];  
}
     

void oss_word_addLetter_r(w,i,r)
     oss_word *w;
     int i;
     reflection *r;
{
  if (w->w.last!=NULL){   // theres is at least one letter
    w->w.last->next=(Pletter) malloc (sizeof(struct pletter));
    //pword_alloc++;
    w->w.last->next->l=i;
    w->w.last->next->next=NULL;
    w->w.last->next->prev=w->w.last;
    w->w.last=w->w.last->next;
  } else { // the first letter
    w->w.last=(Pletter) malloc (sizeof(struct pletter));
    //pword_alloc++;
    w->w.last->l=i;
    w->w.first=w->w.last;
    w->w.last->next=NULL;
    w->w.last->prev=NULL;
  }
  w->count[i]+=w->mod2;
  w->mod2*=-1;
  w->w.l++;
  timesEquals(&(w->m), &(r->m[i])); 
}

// returns 0 if passes
// 1 if fails
int oss_lexi_test(w)
     oss_word *w;
{
  Pletter l,m;
  int i;
  if (w->w.l<=2){
    if (w->w.l==2){
      w->lexi=w->w.first;
      w->lexi_count=0;
    }
    return 0;
  }
    
  if (w->lexi->l > w->w.last->l)
    return 1;  //fails

  //passes w.first part of test!
  if (w->lexi->l == w->w.last->l){
    w->lexi=w->lexi->next;
    w->lexi_count++;
  } else {
    // otherwise (w->lexi->l < w->w.last->l)
    w->lexi_count=0; 
    w->lexi=w->w.first;
  }
  
  //start the second part
  l=w->w.first;
  m=w->w.last;
  i=0;
  while ((l->l == m->l)&&(i++<w->w.l/2)){
    l=l->next;
    m=m->prev;
  }
  if (l->l > m->l)
    return 1; //fails
  return 0;
}


int main() {
  // Let w=l_1...l_n be an odd word
  // Let T_0... T_n be the unfolding, so that 
  // T_0 and T_n differ by an odd number of reflections in edges
  // (a,b,c) will be the coordinates of T_n relative to T_0
  // as displayed in the word.t display window
  int a,b,c, dist;

  int lim1,lim2;
  
  double angA,angB,ang; // records +- 2*angle of triangle
  
  pcomplex tri, temp;
  double d;
  int DEPTH; // Maximal length of word

  triangle T1,T2,Tstart;
  reflection R1,Rstart;
  matrix M;

  word_link SUCCEED;
  oss_word *TRY, *w;
  word new;
  Pletter l;
  int TRY_l, i, dont_continue;
  
  fscanf(stdin,"%lf",&(tri.x)); // the first angle of our triangle
  fscanf(stdin,"%lf",&(tri.y)); // the second angle of our triangle
  fscanf(stdin,"%d",&DEPTH);    // the depth to search to

  T1=triangle_make2(tri);         // initial triangle
  reflect_from_triangle(T1,&R1); // array of reflections in edges
  //printf("(P1x,P1y)=(%lf,%lf)\n", T1.z[0].x, T1.z[0].y);
  //printf("(P2x,P2y)=(%lf,%lf)\n", T1.z[1].x, T1.z[1].y);
  //printf("(P3x,P3y)=(%lf,%lf)\n", T1.z[2].x, T1.z[2].y);
  T2.z[0]=eucTimes(&(R1.m[0]),T1.z[0]);
  T2.z[1]=eucTimes(&(R1.m[0]),T1.z[1]);
  T2.z[2]=eucTimes(&(R1.m[0]),T1.z[2]);
  // Now T2 stores the triangle in position (1,0,0)

  angA=-tri.y*Pi;   // records +- 2*angle of triangle
  angB= tri.x*Pi;   // records +- 2*angle of triangle
  // The triangle in position (a,b,c) w/ a+b+c=0 
  // can be obtained from the triangle in position (0,0,0)
  // by rotating by angle a*angA+b*angB

  DEPTH=DEPTH/2; // we are searching for the square root of the word
  DEPTH=DEPTH-(1-DEPTH%2); // if DEPTH is even decrease by 1
  // Now depth stores the maximal length of the square root 

  TRY=(oss_word *)malloc(DEPTH * sizeof(oss_word));
  SUCCEED=NULL;
  
  for (a=-DEPTH/2; a<=DEPTH/2+1; a++) {
    if (a>=0){
      lim1=-DEPTH/2;
      lim2=DEPTH-a+lim1;
    } else { 
      // a<0
      lim1=-(DEPTH/2+a);
      lim2=DEPTH/2+1;
    }
    for (b=lim1; b<=lim2; b++){
      c=1-a-b;
      // Now (a,b,c) has an odd distance <=DEPTH from the origin
      //printf("Coordinates: (%d, %d, %d)  ",a,b,c);

      ang=(a-1)*angA+b*angB; // Rotation angle between (1,0,0) and (a,b,c)
      ang=fmod(ang,2*Pi); // reduce angle mod 2*Pi
      if (ang<0)
	ang+=2*Pi; // now the angle should be between 0 and 2PI
      ang=-ang/2;
      //printf("Angle to rotate by=%lf*Pi\n",ang/Pi);

      // The axis of the transformation taking (0,0,0) to (a,b,c)
      // has angle -ang... this is well defined
           
      // now rotate so that this axis is horizontal:
      M.m[0][0]=M.m[1][1]=cos(ang);
      M.m[1][0]=sin(ang);
      M.m[0][1]=-sin(ang);
      M.m[0][2]=M.m[1][2]=0;
      Tstart.z[0]=eucTimes(&M,T1.z[0]);
      Tstart.z[1]=eucTimes(&M,T1.z[1]);
      Tstart.z[2]=eucTimes(&M,T1.z[2]);
      reflect_from_triangle(Tstart,&Rstart); // array of reflections in edges
      
      oss_word_empty(&TRY[0]);
      oss_word_addLetter_r(&TRY[0],0,&Rstart);
      TRY[0].interval[0]=Tstart.z[1].y;
      TRY[0].interval[1]=Tstart.z[2].y;
      TRY_l=1;
      dont_continue=0;
      
      while (TRY_l>0){
	w=&(TRY[TRY_l-1]);
	//if ((a==1)&&(b==1)){
	//word_print_file(stdout,&(w->w));
	//printf("Interval=[%lf,%lf]\n",w->interval[0],w->interval[1]);
	//}
	// compute the distance from target (a,b,c)
	dist=0;
	if (w->count[0]-a>=0) 
	  dist+=w->count[0]-a;
	else 
	  dist-=w->count[0]-a;
	if (w->count[1]-b>=0) 
	  dist+=w->count[1]-b;
	else 
	  dist-=w->count[1]-b;
	if (w->count[2]-c>=0) 
	  dist+=w->count[2]-c;
	else 
	  dist-=w->count[2]-c;
	// now distance is stored in dist

	if (dist==0) { 
	  // then we have a candidate-- perform the strong test

	  temp=eucTimes(&(w->m),Tstart.z[0]);
	  d=(temp.y+Tstart.z[0].y)/2;
	  // The axis for the glide reflection is the line y=d

	  // The word passes the strong test iff d is in the interval
	  if ((d>w->interval[0])&&(d<w->interval[1])){
	    word_copy(&(w->w),&new);
	    // make from palindrome form into standard form
	    l=new.first;
	    for (i=new.l; i>0; i--) {
	      word_addLetter(&new,l->l);
	      l=l->next;
	    }	  
	    word_simplify(&new);
	    SUCCEED=word_insert(SUCCEED,&new);
	    //printf("HOORRAY!!!\n");
	  
	    oss_word_destroy(w);
	    TRY_l--;
	    dont_continue=1;
	  }
	}
	if (dont_continue)
	  dont_continue=0;
	else if (w->w.l+dist>DEPTH) {
	  // we have strayed too far from the target-- burn the word
	  oss_word_destroy(w);
	  TRY_l--;
	} else if (oss_lexi_test(w)) {
	  // Fails lexi test
	  oss_word_destroy(w);
	  TRY_l--;
	} else {
	  // record the location of the newest vertex in temp
	  temp=eucTimes(&(w->m),Tstart.z[w->w.last->l]);
	  d=temp.y;
	  //if ((a==1)&&(b==1)){
	  //  printf("(x,y)=(%lf,%lf)\n", temp.x, temp.y);
	  //}
	  if (d<=w->interval[0]){
	    if (w->w.l%2==0){
	      oss_word_addLetter_r(w,(w->w.last->l+2)%3,&Rstart);
	    } else {
	      oss_word_addLetter_r(w,(w->w.last->l+1)%3,&Rstart);
	    }
	  } else if (d>=w->interval[1]){
	    if (w->w.l%2==0){
	      oss_word_addLetter_r(w,(w->w.last->l+1)%3,&Rstart);
	    } else {
	      oss_word_addLetter_r(w,(w->w.last->l+2)%3,&Rstart);
	    }	    
	  } else { // d lies in the interval split the word
	    oss_word_copy(w,&(TRY[TRY_l]));
	    
	    if (w->w.l%2==0){
	      oss_word_addLetter_r(w,(w->w.last->l+1)%3,&Rstart);
	      oss_word_addLetter_r(&TRY[TRY_l],
				   (TRY[TRY_l].w.last->l+2)%3,&Rstart);
	    } else {
	      oss_word_addLetter_r(w,(w->w.last->l+2)%3,&Rstart);
	      oss_word_addLetter_r(&TRY[TRY_l],
				   (TRY[TRY_l].w.last->l+1)%3,&Rstart);
	    }
	    w->interval[1]=d;
	    TRY[TRY_l].interval[0]=d;
	    TRY_l++;
	  }
	}	
      }
    }
  }
  word_link_print_file(stdout, SUCCEED);
  word_link_destroy(SUCCEED);

  return 0;
}