/* 
 * This is a program to draw fundamantal domains
 * this is version 1.0
 */

/* 
 * This program is written in c++, though actually it
 * is written in c, with just a few c++ commands thrown
 * in, and could easily be converted back to c, should
 * you desire.
 *
 * This program was written on a linux machine.
 * It also works on unix systems I've tried it on.  
 * A very slightly modified version works with windows.
 * You are free to modify it to work on whatever system you're 
 * running.
 *
 */

/*
 * Simple instructions to compile and run under unix type systems:
 * 
 * Give this file the name fundomain.c
 * Compile with: gcc -o drawfd fundomain.c -lm 
 * This creates the file "drawfd"
 * then run by typing: drawfd 
 * or, if that doesn't work try: ./drawfd
 * after pressing return, the program runs, and    
 * you will be prompted for whether you want a picture of
 * Gamma_0(N) or Gamma_1(N), and for the value of N,
 * and for the name of the postscript putput file.
 *
 */

/*
 * This program was written by Helena A. Verrill, 2000
 * Date of most recent modification: June 28 2000.
 */

/* #include </usr/include/math.h> */
#include <math.h>
#include <stdio.h>

#define PI 3.14159265358979
#define infinity 600
/* #define scale 150   */


#define max 50

/* note, in c on windows machine, the following three
 * matrices were declared to be "const"
 * but this does not seem to work on unix.
 */

int T[2][2] = {{1,1},{0,1}};

int Tm[2][2] = {{1,-1},{0,1}};

int S[2][2] = {{0,1},{-1,0}};

int proj_st[max][max][3];
int scale = 150;

void findfund(int N, int m[2][2], int done[max][max][2][3], int dist,int sbsp);
int  find_red_uv(int sign, int m2[2][2], int N, int comp);


FILE           *output;		/* the output file */
FILE           *outputtxt;		/* the output file */

void
do_point(double x, double y, char *s)
{
	fprintf(output, "%.8lg %.8lg %s ", x, y, s);
}

void
my_arc(double thrang, double startang, double radius)
{
	double sign;

	if (radius<0){
		  thrang=thrang*-1;
		  radius=-1*radius;}
	fprintf(output, "/thr-ang  %.8lg def\n",thrang);
	if (thrang<0) {sign=-1;}
	if (thrang>0) {sign=1;}
	fprintf(output, "/sign  %.8lg def\n",sign);
	fprintf(output, "/start-ang  %.8lg def\n",startang);
	fprintf(output, "/radius  %.8lg def\n",radius);
	fprintf(output, "myarc\n");
}

void
closed_line(double x1, double y1, double x2, double y2)
{
	fprintf(output, "newpath ");
	do_point(x1, y1, "moveto");
	do_point(x2, y2, "lineto");
	fprintf(output, "stroke\n");
}

void
make_title(char thescript, char thetype, int Anumber)
{  	
  fprintf(output,"/textfontsize %d def\n",4*Anumber);
  fprintf(output,"/Helvetica-BoldOblique findfont textfontsize scalefont setfont\n");
  fprintf(output,"100 -100 moveto\n");
  fprintf(output,"(Gamma%c%c(%d)) show\n\n",thescript,thetype,Anumber);
} 
   
  

void
open_line(double x1, double y1, double x2, double y2)
{
	do_point(x1, y1, "lineto");
	do_point(x2, y2, "lineto");
}

double
matmulty(double a, double b, double c, double d,double x, double y)
{
	return a*x+b*y;
}

double
matmultx(double a, double b, double c, double d,double x, double y)
{
	return c*x+d*y;
}

void
triA(double ang1, double ang2, double ang3, double r1, double r2, double r3, double x, double y,double red, double green, double blue)
{
	fprintf(output, "gsave\n");
	fprintf(output, "newpath\n");
	fprintf(output, "%.8lg %.8lg  moveto\n",x,y);

	fprintf(output, "/r1 %.6lg abs def\n",r1);
	fprintf(output, "/r2 %.6lg abs def\n",r2);
	fprintf(output, "/r3 %.6lg abs def\n",r3);
	fprintf(output, "/a1 %.6lg %d mul def\n",ang1,2*(r1>0)-1);
	fprintf(output, "/a2 %.6lg %d mul def\n",ang2,2*(r2>0)-1);
	fprintf(output, "/a3 %.6lg %d mul def\n",ang3,2*(r3>0)-1);
	fprintf(output, "%.4lg %.4lg %.4lg %.4lg %.4lg mytriA fill\n",red,green,blue,x,y);
	fprintf(output, "0 0 0 %.4lg %.4lg mytriA stroke grestore\n",x,y);
	fprintf(output, "\n");
}

void
triB
(double ang1, double length2, double ang3, double r1, double r3, double x, double y,double red, double green, double blue)
{
	fprintf(output, "gsave\n");
	fprintf(output, "newpath\n");
	fprintf(output, "%.8lg %.8lg  moveto\n",x,y);
	fprintf(output, "/r1 %.6lg abs def\n",r1);
	fprintf(output, "/r3 %.6lg abs def\n",r3);
	fprintf(output, "/a1 %.6lg %d mul def\n",ang1,2*(r1>0)-1);
	fprintf(output, "/a3 %.6lg %d mul def\n",ang3,2*(r3>0)-1);
	fprintf(output, "%.4lg %.4lg %.4lg 0 %.4lg abs %.4lg 0 mytriB fill\n",red,green,blue,length2,x);
	fprintf(output, "0 0 0 0 %.4lg abs %.4lg 0 mytriB stroke grestore\n",length2,x);
	fprintf(output, "\n");
}

void
triC
(double length1, double ang2, double ang3, double r2, double r3, double x, double y,double red, double green, double blue)
{
	fprintf(output, "gsave\n");
	fprintf(output, "newpath\n");
	fprintf(output, "%.8lg %.8lg  moveto\n",x,y);
	fprintf(output, "/r2 %.6lg abs def\n",r2);
	fprintf(output, "/r3 %.6lg abs def\n",r3);
	fprintf(output, "/a2 %.6lg %d mul def\n",ang2,2*(r2>0)-1);
	fprintf(output, "/a3 %.6lg %d mul def\n",ang3,2*(r3>0)-1);
	fprintf(output, "%.4lg %.4lg %.4lg 0 %.4lg abs %.4lg 0 mytriC fill\n",red,green,blue,length1,x);
	fprintf(output, "0 0 0 0 %.4lg abs %.4lg 0 mytriC stroke grestore\n",length1,x);
	fprintf(output, "\n");

}

void
triD
(double ang1, double ang2, double r1, double r2, double x, double y,double red, double green, double blue)
{

	fprintf(output, "gsave\n");
	fprintf(output, "newpath\n");
	fprintf(output, "%.8lg %.8lg  moveto\n",x,y);
	fprintf(output, "/r1 %.6lg abs def\n",r1);
	fprintf(output, "/r2 %.6lg abs def\n",r2);
	fprintf(output, "/a1 %.6lg %d mul def\n",ang1,2*(r1>0)-1);
	fprintf(output, "/a2 %.6lg %d mul def\n",ang2,2*(r2>0)-1);
	fprintf(output, "%.4lg %.4lg %.4lg %.4lg 0 mytriD fill\n",red,green,blue,x);
	fprintf(output, "0 0 0 %.4lg 0 mytriD stroke grestore\n",x);
	fprintf(output, "\n");

}


void 
triE
(double r2, double ang2, double x, double el,double red, double green, double blue)
{
	fprintf(output, "gsave\n");
	fprintf(output, "newpath\n");
	fprintf(output, "%.8lg 0  moveto\n",x);
	fprintf(output, "/r2 %.6lg abs def\n",r2);
	fprintf(output, "/a2 %.6lg %d mul def\n",ang2,2*(r2>0)-1);
	fprintf(output, "%.4lg %.4lg %.4lg 0 %.4lg abs dup %.4lg mytriE fill\n",red,green,blue,el,x);
	fprintf(output, "0 0 0 0 %.4lg abs dup %.4lg mytriE stroke grestore\n",el,x);
	fprintf(output, "\n");
}


void
transtri(double a,double b,double c,double d, double red, double green, double blue)
{
	double xm1,ym1,xp1,yp1,xmh,ymh,xph,yph,xinf,yinf,m1 ,p1 ,mh ,ph, inf;
	double el1, el2, elI1, elI2, dist1, dist2a, dist2b, dist3;
	double center1, center2, center3;
	double ang1, ang2, ang3;
	double radius1, radius2, radius3;
	
	if(a*d-c*b>0){ 
		xm1  = matmulty(a,b,c,d,-1,1);   /* image of -1 */
		ym1  = matmultx(a,b,c,d,-1,1);
		xp1  = matmulty(a,b,c,d,1,1);    /* image of 1 */
		yp1  = matmultx(a,b,c,d,1,1);
		xmh  = matmulty(a,b,c,d,-1,2);   /* image of -1/2 */
		ymh  = matmultx(a,b,c,d,-1,2);
		xph  = matmulty(a,b,c,d,1,2);    /* image of 1/2 */
		yph  = matmultx(a,b,c,d,1,2);
		xinf = matmulty(a,b,c,d,1,0);    /* image of infinity */
		yinf = matmultx(a,b,c,d,1,0);
		m1=0; p1=0; mh=0; ph=0; inf=0;
		if (ym1!=0){
			m1  =  scale * xm1/ym1;}
		if (yp1!=0){
			p1  =  scale * xp1/yp1;}
		if (ymh!=0){
			mh  =  scale * xmh/ymh;}
		if (yph!=0){
			ph  =  scale * xph/yph;}
		if (yinf!=0){
			inf =  scale * xinf/yinf;}
		/* define real and imag parts of elliptic points */
		el1 =  scale * ((a+2*b)*(2*d+c)+3*a*c)/((2*d+c)*(2*d+c)+3*c*c);
		el2 =  scale * ((-a+2*b)*(2*d-c)+3*a*c)/((2*d-c)*(2*d-c)+3*c*c);
		elI1 =  scale * sqrt(3)*2*(a*d-b*c)/((2*d+c)*(2*d+c)+3*c*c);
		elI2 =  scale * sqrt(3)*2*(a*d-b*c)/((2*d-c)*(2*d-c)+3*c*c);
		
		center1 = (inf+ph)/2;
		center2 = (inf+mh)/2;
		center3 = (m1+p1)/2;
		radius1 = (inf-ph)/2;
		radius2 = (m1-p1)/2;
		radius3 = (mh-inf)/2;
		
		dist1 = el1-center1;   
		dist3 = center2 - el2;
		dist2a = center3-el1;
		dist2b = el2-center3;
		
		ang1 = (acos(dist1/radius1))*180/PI;
		ang3 = (acos(dist3/radius3))*180/PI;
		ang2 = -(asin(dist2a/radius2)+asin(dist2b/radius2))*180/PI;

		if (ym1*yp1*ymh*yph*yinf!=0){
			triA(ang1,ang2,ang3,radius1,radius2,radius3,inf,0,red,green,blue);
		}
		if (ym1*yp1==0){
			triB(ang1,elI2-elI1,ang3,radius1,radius3,inf,0,red,green,blue);
		}
		if (yph==0){
			triC(elI1,ang2,ang3,radius2,radius3,inf,0,red,green,blue);
		}
		if (ymh==0){
			triD(ang1,ang2,radius1,radius2,inf,0,red,green,blue);
		} 
		if (yinf==0){
			triE(radius2,ang2,el1,elI1,red,green,blue);
		}
	}
}


void 
lots_of_triangles()
{
	double          a,b,c,d,e;        
	e=1;
	for (a = 1; a < 7; a = a + 1) {
		for (b = -1; b < 6; b = b + 1) {
			for (c = -1; c < 5; c = c + 1) {
				for (d = -1; d < 4 ; d = d + 1) {
					e=e+0.4;
					transtri(a,b,c,d,1,0.5+cos(e),0.5+sin(e));
	}}}}	
}

void
more_trianglesa()
{
	double          a,e;
	e=0;
	for (a = -5; a < 5; a = a + 1) {
		  e=e+.1;
		transtri(0,1,-1,a,(1-e)*(1-e)+1-(.5-e)*(.5-e),1-e*e,e*e);
		transtri(a,1,-1,0,(1-e)*(1-e)+1-(.5-e)*(.5-e),1-e*e,e*e);
		transtri(1,a,0,1,(1-e)*(1-e)+1-(.5-e)*(.5-e),1-e*e,e*e);
		transtri(1,0,a,1,(1-e)*(1-e)+1-(.5-e)*(.5-e),1-e*e,e*e);
	}
}

void 
more_triangles()
{
	double          a,e;
	e=0;
	for (a = -17; a < 18; a = a + 1) {
		  e=1-e;
		transtri(0,1,-1,a,e,0,1);
		if(a*a!=1){transtri(1,0,a,1,1-e,0,1);}
	}
	e=0;
	for (a = -5; a < 6; a = a + 1) {
		  e=1-e;
		if(a<-2 || a>1) {transtri(1,1,a,a+1,e,0,1);}
		if(a<-1 || a>2) {transtri(1,-1,a,-a+1,e,0,1);}
	}
}

void
row()
{
	 fprintf(output, "/block{\n");
	 more_triangles();

	 fprintf(output, "} def\n");
	 fprintf(output, "\n");
	 fprintf(output, "\n");
	 fprintf(output,"-120 0 translate \n");
	 fprintf(output,"/col 0 def \n");
	 fprintf(output,"4 {\n");
	 fprintf(output,"block\n");
	 fprintf(output,"150 0 translate \n");
	 fprintf(output,"} repeat\n");

}


void
more_trianglesb()
{
	double          a,e;
	e=0;
	for (a = -5; a < 5; a = a + 1) {
		  e=e+.1;
		transtri(0,1,-1,a,e,0,1);
		transtri(a,1,-1,0,1,0,e);
		transtri(1,a,0,1,1,e,0);
		transtri(1,0,a,1,0,1,e);
	}
}


int coprime(int a[3])
{
	int b[3];

	/* this function only works if it is fed non negative integers! */
	if((a[1]*a[1]-1)*(a[0]*a[0]-1)*(a[2]*a[2]-1)==0){return 1;}
	if(a[0]==0)
	 {
		if(a[1]==0)
		{return 0;}
		else{a[0]=a[1];}
	 }
	b[0]=a[1]-(a[1]/a[0])*a[0];
	b[1]=a[2]-(a[2]/a[0])*a[0];
	b[2]=a[0];
	return coprime(b);
}


void
Proj_reps(int N, int mat[10][10])
{
	int j,k,p;

	for(j=1;j<N;++j){
		for(k=1;k<N;++k){
			if(mat[j][k]==1){
				for (p=2;p<N;++p){
					if(mat[0][p]==1){
						mat[j*p-((j*p)/N)*N][k*p-((k*p)/N)*N]=0;
					}
				}
			}
		}
	}
}



void
matinit(int m1[2][2], int a, int b, int c, int d)
{

	m1[0][0]=a;
	 m1[0][1]=b;
	m1[1][0]=c;
	m1[1][1]=d;

}





void
matmult(int c[2][2], int a[2][2], int b[2][2])
{

	c[0][0]=a[0][0]*b[0][0]+a[0][1]*b[1][0];
	c[0][1]=a[0][0]*b[0][1]+a[0][1]*b[1][1];
	c[1][0]=a[1][0]*b[0][0]+a[1][1]*b[1][0];
	c[1][1]=a[1][0]*b[0][1]+a[1][1]*b[1][1];

}


void
printmat(int m1[2][2])
{
	int i,j;
	
	for (i=0;i<2;++i)
	{
		printf("\n");
		for (j=0;j<2;++j)
		{
			printf(" %d",m1[i][j]);
		}
	}
	printf("\n");
}





void
bezout(int pair[2])
{
	int u,v,v1,i,j;
	int newmat[2][2], runningmat[2][2], newrun[2][2];



	u=pair[0];
	v=pair[1];

	matinit(newmat,0,1,1,0);
	printmat(newmat);
	 matinit(runningmat,1,0,0,1);

	printmat(runningmat);

	while(u>0)
	 {


		newmat[1][1]=-(v/u);

		matmult(newrun,runningmat,newmat);

		for (i=0;i<2;++i){
			for (j=0;j<2;++j){
				runningmat[i][j]=newrun[i][j];

			}
		}

		printmat(runningmat);

		v1=u;
		u=v-(v/u)*u;
		v=v1;


	 }

	printmat(runningmat);


	pair[0]=runningmat[1][0];
	pair[1]=runningmat[0][0];

}



void
listpairs(int N, int pairlist[10][10], int matlist[10][10][2])
{
	int i,j, b[3];
	int pair[2];
	
	/* matlist contains the bottom parts of the matrix with top tow [i,j] */
	
	N=6;
	
	for (i=0;i<N;++i){
		for (j=0;j<N;++j){
			pairlist[i][j]=0;
			b[0]=i;
			b[1]=j;
			b[2]=N;
			pairlist[i][j]=coprime(b);
			
		}
	}
	
	
	Proj_reps(N,pairlist);
	
	for (i=0;i<N;++i){
		for (j=0;j<N;++j){
			if (i>0 && j>0 && pairlist[i][j]!=0) 
			{	pair[0]=i;
			pair[1]=j;
			
			bezout(pair);
			matlist[i][j][0]=-pair[1];
			matlist[i][j][1]=pair[0];

			}
			else {
				matlist[i][j][0]=0;
				matlist[i][j][1]=0;
			}
		}
	}
	
}

void
proj_eq(int N, int kind)
{
	
	int i,j,k,u,v;
	
	
	for(i=1;i<N;++i) {
		for(j=1;j<N;++j) {
										  for(k=0;k<2;++k) {
											  proj_st[i][j][k]=0;
	}}}
	
	
	if (kind==0){
		for(i=1;i<N;++i) {
			proj_st[0][i][1]=1;
			proj_st[i][0][0]=1;
		}



		for(i=1;i<N;++i) {
			for(j=1;j<N;++j) {
				if(proj_st[i][j][0]==0 && proj_st[i][j][1]==0){
					proj_st[i][j][0]=i;
					proj_st[i][j][1]=j;
					for(k=1;k<N;++k){
						u=k*i-((k*i)/N)*N;
						v=k*j-((k*j)/N)*N;
						if(u!=0 && v!=0){
							if (proj_st[u][v][0]==0&&proj_st[u][v][1]==0){
								proj_st[u][v][0]=i;
								proj_st[u][v][1]=j;
		}}}}}}

	}
	
	if(kind==1){
		for(i=1;i<N/2;++i) {
			proj_st[0][i][1]=i;
			proj_st[i][0][0]=i;
		}
		for(i=N/2;i<N;++i) {
			proj_st[0][i][1]=N-i;
			proj_st[i][0][0]=N-i;
		}



		for(i=1;i<N;++i) {
			for(j=1;j<N;++j) {
				if(j<=N/2) {
				    proj_st[i][j][0]=i;
					proj_st[i][j][1]=j;}
				else {
					proj_st[i][j][0]=N-i;
					proj_st[i][j][1]=N-j;
				}
			}
	
		}

		if ((N/2)*2==N){
			for(i=1;i<N/2;++i) {
				proj_st[N/2][i][0]=N/2;
				proj_st[N/2][i][1]=i;
				proj_st[i][N/2][0]=i;
				proj_st[i][N/2][1]=N/2;
			}
			for(i=N/2;i<N;++i) {
			    proj_st[N/2][i][0]=N/2;
				proj_st[N/2][i][1]=N-i;
				proj_st[i][N/2][0]=N-i;
				proj_st[i][N/2][1]=N/2;
			}
		}


	}
	
	
}






void 
shorten(int N, int done[max][max][2][3])
{

	int s,i,j, u,v,u1,v1;
	int m[2][2],m1[2][2],m2[2][2],m3[2][2];




	s=0;
	
	for(i=0;i<N;++i){
		for(j=0;j<N;++j){
			if (done[i][j][0][2] > s) {
				s=done[i][j][0][2]; u=i; v=j;
			}
	}}

	if (s<1) return;
	for(i=0;i<2;++i){for(j=0;j<2;++j){m[i][j]=done[u][v][i][j];}}


	matmult(m1,S,m);
		
	u1=find_red_uv(1,m1,N,0);
	v1=find_red_uv(1,m1,N,1);

	if (done[u1][v1][0][0]==m1[0][0] && done[u1][v1][0][1]==m1[0][1] && done[u1][v1][1][0]==m1[1][0] && done[u1][v1][1][1]==m1[1][1]) {
		if (done[u1][v1][0][2]<s-1) {
			for(i=0;i<2;++i){for(j=0;j<2;++j){m2[i][j]=done[u1][v1][i][j];}}
			matmult(m3,S,m2);
			for(i=0;i<2;++i){for(j=0;j<2;++j){done[u][v][i][j]=m3[i][j];}}
			done[u][v][0][2]=done[u1][v1][0][2]+1;
			shorten(N,done);
	}}

	matmult(m1,T,m);
		
	u1=find_red_uv(1,m1,N,0);
	v1=find_red_uv(1,m1,N,1);

	if (done[u1][v1][0][0]==m1[0][0] && done[u1][v1][0][1]==m1[0][1] && done[u1][v1][1][0]==m1[1][0] && done[u1][v1][1][1]==m1[1][1]) {
		if (done[u1][v1][0][2]<s-1) {
			for(i=0;i<2;++i){for(j=0;j<2;++j){m2[i][j]=done[u1][v1][i][j];}}
			matmult(m3,Tm,m2);
			for(i=0;i<2;++i){for(j=0;j<2;++j){done[u][v][i][j]=m3[i][j];}}
			done[u][v][0][2]=done[u1][v1][0][2]+1;
			shorten(N,done);
	}}


		matmult(m1,Tm,m);
		
	u1=find_red_uv(1,m1,N,0);
	v1=find_red_uv(1,m1,N,1);

	if (done[u1][v1][0][0]==m1[0][0] && done[u1][v1][0][1]==m1[0][1] && done[u1][v1][1][0]==m1[1][0] && done[u1][v1][1][1]==m1[1][1]) {
		if (done[u1][v1][0][2]<s-1) {
			for(i=0;i<2;++i){for(j=0;j<2;++j){m2[i][j]=done[u1][v1][i][j];}}
			matmult(m3,T,m2);
			for(i=0;i<2;++i){for(j=0;j<2;++j){done[u][v][i][j]=m3[i][j];}}
			done[u][v][0][2]=done[u1][v1][0][2]+1;
			shorten(N,done);
	}}

	return;



	

}



int 
find_red_uv(int sign, int m2[2][2], int N, int comp)
{
	
	int uu,vv,u,v;
	
	u=sign*(m2[0][0]-(m2[0][0]/N)*N);
	v=sign*(m2[0][1]-(m2[0][1]/N)*N);
	if(u<0) u=N+u;
	if(v<0) v=N+v;
	uu=u;
	vv=v;
	return proj_st[uu][vv][comp];
	
}



void
nearest(int m[2][2], int A[2][2], int N, int done[max][max][2][3], int dist,int sbsp)
{

	int u,v,i,j;
	int dist1,dist3,dist4;
	int m2[2][2],m3[2][2],m4[2][2];
	int u3,v3,u4,v4;
	int C[2][2], D[2][2], Cm[2][2], Dm[2][2];

		dist3=max*max;
	dist4=max*max;


	if (A[0][0]==T[0][0] && A[0][1]==T[0][1] && A[0][1]==T[0][1] && A[0][1]==T[0][1]) {
		for (i=0;i<2;++i){for (j=0;j<2;++j){C[i][j]=S[i][j];}}
		for (i=0;i<2;++i){for (j=0;j<2;++j){D[i][j]=T[i][j];}}
		for (i=0;i<2;++i){for (j=0;j<2;++j){Cm[i][j]=S[i][j];}}
		for (i=0;i<2;++i){for (j=0;j<2;++j){Dm[i][j]=Tm[i][j];}}
	}

		if (A[0][0]==S[0][0] && A[0][1]==S[0][1] && A[0][1]==S[0][1] && A[0][1]==S[0][1]) {
		for (i=0;i<2;++i){for (j=0;j<2;++j){C[i][j]=T[i][j];}}
		for (i=0;i<2;++i){for (j=0;j<2;++j){D[i][j]=Tm[i][j];}}
		for (i=0;i<2;++i){for (j=0;j<2;++j){Cm[i][j]=Tm[i][j];}}
		for (i=0;i<2;++i){for (j=0;j<2;++j){Dm[i][j]=T[i][j];}}
	}

			if (A[0][0]==Tm[0][0] && A[0][1]==Tm[0][1] && A[0][1]==Tm[0][1] && A[0][1]==Tm[0][1]) {
		for (i=0;i<2;++i){for (j=0;j<2;++j){C[i][j]=S[i][j];}}
		for (i=0;i<2;++i){for (j=0;j<2;++j){D[i][j]=Tm[i][j];}}
		for (i=0;i<2;++i){for (j=0;j<2;++j){Cm[i][j]=S[i][j];}}
		for (i=0;i<2;++i){for (j=0;j<2;++j){Dm[i][j]=T[i][j];}}
	}

	matmult(m2,m,A);


	u=find_red_uv(1,m2,N,0);
	v=find_red_uv(1,m2,N,1);

	if (done[u][v][0][0]==0 && done[u][v][0][1]==0) {
		/*&& done[u1][v1][0][0]==0 && done[u1][v1][0][1]==0) {*/

		if (sbsp==0 && m2[1][0]==0) {
			dist1=0;
			findfund(N,m2,done,dist1,sbsp);
		}
		else {



			dist1=dist+1;

			/* note, divide by 2 on following row because even though the distance as a graph increases, the size gets smaller at
			a non linear rate.  Possibly should change this to a different function that just divide by 2 */

			if (sbsp==1 && m2[0][0]==1) {
				dist1 = m2[0][1];
				if (dist1*dist1 < 5) dist1=dist1/2;
				if (dist1<0) dist1=-dist1;
			}

			matmult(m3,m2,C);
			matmult(m4,m2,D);


			u3=find_red_uv(1,m3,N,0);
			v3=find_red_uv(1,m3,N,1);
	/*		u31=find_red_uv(-1,m3,N,0);
			v31=find_red_uv(-1,m3,N,1);*/
			if (done[u3][v3][0][0]!=0 || done[u3][v3][0][1]!=0) {
				dist3=done[u3][v3][0][2];
				for (i=0;i<2;++i){for (j=0;j<2;++j){m3[i][j]=done[u3][v3][i][j];}}
			}
	/*		if (done[u31][v31][0][0]!=0 || done[u31][v31][0][1]!=0)  {
				dist3=done[u31][v31][0][2];
				for (i=0;i<2;++i){for (j=0;j<2;++j){m3[i][j]=done[u31][v31][i][j];}}
			}
	*/


			u4=find_red_uv(1,m4,N,0);
			v4=find_red_uv(1,m4,N,1);
/*			u41=find_red_uv(-1,m4,N,0);
			v41=find_red_uv(-1,m4,N,1);  */
			if (done[u4][v4][0][0]!=0 || done[u4][v4][0][1]!=0) {
				dist4=done[u4][v4][0][2];
				for (i=0;i<2;++i){for (j=0;j<2;++j){m4[i][j]=done[u4][v4][i][j];}}
			}
/*			if (done[u41][v41][0][0]!=0 || done[u41][v41][0][1]!=0)  {
				dist4=done[u41][v41][0][2];
				for (i=0;i<2;++i){for (j=0;j<2;++j){m4[i][j]=done[u4][v4][i][j];}}
			}
			*/


			if (dist3 < dist1 && dist3 < dist4) {
				matmult(m2,m3,Cm); 
				findfund(N, m2, done, dist3+1,sbsp);
			}

			else if (dist4 <= dist3 && dist4<dist1)
			{
				matmult(m2,m4,Dm); 
				findfund(N, m2, done, dist4+1,sbsp);

			}
			else {
				findfund(N,m2,done,dist1,sbsp);
			}


		}
	}

}



void
findfund(int N, int m[2][2], int done[max][max][2][3], int dist,int sbsp)
{
	int u,v,Nu,Nv,u1,v1,i,j,ua,va;


	u=m[0][0]-(m[0][0]/N)*N;
	v=m[0][1]-(m[0][1]/N)*N;
	if(u<0) u=N+u;
	if(v<0) v=N+v;




	ua=proj_st[u][v][0];
	va=proj_st[u][v][1];

	u=ua;
	v=va;


	Nu=(N-u)-((N-u)/N)*N;
	Nv=(N-v)-((N-v)/N)*N;

	u1=proj_st[Nu][Nv][0];
	v1=proj_st[Nu][Nv][1];


	for (i=0;i<2;++i){
		for (j=0;j<2;++j){

			done[u][v][i][j]=m[i][j];
			if (done[u1][v1][i][j]!=done[u][v][i][j]);
			done[u1][v1][i][j]=m[i][j];

		}
	}
	 done[u][v][0][2]=dist;

/*	fprintf(outputtxt," u1 %d, %d, %d %d %d\n",u,v, done[u][v][0][0], done[u][v][0][1],dist); */

	/* CASE 1 */


	nearest(m,T,N,done,dist,sbsp);

			/* CASE 2 */

	nearest(m,S,N,done,dist,sbsp);
	nearest(m,Tm,N,done,dist,sbsp);

}

void
topmatter(int N,int sbsp)
{
	double sc;
	int tra;
	int scale;

	scale = 150; /* scale tells how accurately the circles are drawn */
    if (sbsp == 1) scale=400;

    tra=60; 
	sc=(float)3.5/N;
	if(sbsp==1) {sc = 3; tra=300;}


	fprintf(output, "%%!PS \n");
	if(sbsp==1) {fprintf(output, "0.2 setlinewidth \n");}
	fprintf(output, "%d 300 translate \n",tra);
	fprintf(output, "\n");
	fprintf(output, "%.4lg %.4lg scale\n",sc,sc);
	fprintf(output, "\n");
	fprintf(output, "/infinity 200 def\n");
	fprintf(output, " \n");
	fprintf(output, "/circle {    \n");
	fprintf(output, "   0 0 rlineto \n");
	fprintf(output, "   0 chord rlineto           \n");
	fprintf(output, "{\n");
	fprintf(output, "   rotate                \n");
	fprintf(output, "   0 chord rlineto           \n");
	fprintf(output, "} repeat\n");
	fprintf(output, "}bind def \n");
	fprintf(output, "\n");
	fprintf(output, "\n");
	fprintf(output, "/myarc{\n");
	fprintf(output, "/N %d def \n",scale);
	fprintf(output, "/turn thr-ang N div def\n");
	fprintf(output, "/chord sign 2 radius turn 2 div sin mul mul mul def\n");
	fprintf(output, "start-ang rotate\n");
	fprintf(output, "N {turn}repeat N  circle\n");
	fprintf(output, "180 rotate  \n");
	fprintf(output, "} def\n");
	fprintf(output, "\n");
	fprintf(output, "\n");
	fprintf(output,"/mytriA{\n");
	fprintf(output,"newpath\n");
	fprintf(output,"moveto\n");
	fprintf(output,"/thr-ang  a1 def\n");
	fprintf(output,"/sign  a1 a1 abs div round def\n");
	fprintf(output,"/start-ang  0 def\n");
	fprintf(output,"/radius  r1 def\n");
	fprintf(output,"myarc\n");
	fprintf(output,"/thr-ang  a2 def\n");
	fprintf(output,"/sign  a2 a2 abs div round def\n");
	fprintf(output,"/start-ang  60 def\n");
	fprintf(output,"/radius  r2 def\n");
	fprintf(output,"myarc\n");
	fprintf(output,"/thr-ang  a3 def\n");
	fprintf(output,"/sign  a3 a3 abs div round  def\n");
	fprintf(output,"/start-ang  60 def\n");
	fprintf(output,"/radius  r3 def\n");
	fprintf(output,"myarc\n");
	fprintf(output,"setrgbcolor\n");
	fprintf(output,"} def\n");
	fprintf(output, "\n");
	fprintf(output, "\n");
	fprintf(output,"/mytriB{\n");
	fprintf(output,"newpath\n");
	fprintf(output,"moveto\n");
	fprintf(output,"/thr-ang  a1 def\n");
	fprintf(output,"/sign  a1 a1 abs div round def\n");
	fprintf(output,"/start-ang  0 def\n");
	fprintf(output,"/radius  r1 def\n");
	fprintf(output,"myarc\n");
	fprintf(output,"60 rotate \n");
	fprintf(output,"rlineto \n");
	fprintf(output,"/thr-ang  a3 def\n");
	fprintf(output,"/sign  a3 a3 abs div round  def\n");
	fprintf(output,"/start-ang  240 def\n");
	fprintf(output,"/radius  r3 def\n");
	fprintf(output,"myarc\n");
	fprintf(output,"setrgbcolor\n");
	fprintf(output,"} def\n");
	fprintf(output, "\n");
	fprintf(output, "\n");
	fprintf(output,"/mytriC{\n");
	fprintf(output,"newpath\n");
	fprintf(output,"moveto\n");
	fprintf(output,"rlineto \n");
	fprintf(output,"/thr-ang  a2 def\n");
	fprintf(output,"/sign  a2 a2 abs div round def\n");
	fprintf(output,"/start-ang  240 def\n");
	fprintf(output,"/radius  r2 def\n");
	fprintf(output,"myarc\n");
	fprintf(output,"/thr-ang  a3 def\n");
	fprintf(output,"/sign  a3 a3 abs div round  def\n");
	fprintf(output,"/start-ang  60 def\n");
	fprintf(output,"/radius  r3 def\n");
	fprintf(output,"myarc\n");
	fprintf(output,"setrgbcolor\n");
	fprintf(output,"} def\n");
	fprintf(output, "\n");
	fprintf(output, "\n");
	fprintf(output,"/mytriE{\n");
	fprintf(output,"newpath\n");
	fprintf(output,"1 copy %d moveto\n",infinity);
	fprintf(output,"exch lineto \n");
	fprintf(output,"/thr-ang  a2 def\n");
	fprintf(output,"/sign  a2 a2 abs div round def\n");
	fprintf(output,"/start-ang  60 def\n");
	fprintf(output,"/radius  r2 def\n");
	fprintf(output,"myarc\n");
	fprintf(output,"60 rotate \n");
	fprintf(output,"%d exch sub rlineto \n",infinity);
	fprintf(output,"setrgbcolor\n");
	fprintf(output,"} def \n");
	fprintf(output,"\n");
	fprintf(output,"\n");
	fprintf(output,"/mytriD{\n");
	fprintf(output,"newpath\n");
	fprintf(output,"moveto\n");
	fprintf(output,"/thr-ang  a1 def\n");
	fprintf(output,"/sign  a1 a1 abs div round def\n");
	fprintf(output,"/start-ang  0 def\n");
	fprintf(output,"/radius  r1 def\n");
	fprintf(output,"myarc\n");
	fprintf(output,"/thr-ang  a2 def\n");
	fprintf(output,"/sign  a2 a2 abs div round def\n");
	fprintf(output,"/start-ang  60 def\n");
	fprintf(output,"/radius  r2 def\n");
	fprintf(output,"myarc\n");
	fprintf(output,"closepath\n");
	fprintf(output,"setrgbcolor\n");
	fprintf(output,"240 rotate \n");
	fprintf(output,"} def\n");
	fprintf(output,"\n");
	fprintf(output,"\n");

}



int dopsfile()
{
	int i,j,k,l, sbsp, N, gotsbsp, gotN, gotName;
	char ch, thescript, thetype;
	int kind, gotkind;
	int temp;
	int done[max][max][2][3],m1[2][2],mt[2][2],mm[2][2];
	char GammaNps[256];

	gotN = 0;
	
	printf("enter N (0 < N < %d)  (or type q at any time to quit)\n",max);
	
	N = 0;
	while(gotN==0){
	  if (scanf("%d",&N)!=EOF) 
	    {if (N>0 && N <50) gotN = 1;}
	  if (gotN==0){
	    scanf("%c",&ch);
	    if (tolower(ch)=='q') return 0; 
	    if (ch!='\n') while (ch!='\n') scanf("%c",&ch);
	    printf("please enter q or an integer N, 0 < N < 50\n");
	  }
	}
	    

	if(N==0) return 0;
	if(N>0&&N<=max){
	        gotsbsp = 0;
		printf("upper or lower subscript? (u or l)\n");
		scanf("%c",&ch);

		while (gotsbsp==0){
		  if (tolower(ch)=='u') {sbsp=0; thescript='^'; gotsbsp=1;}
		  else if (tolower(ch)=='l'){sbsp=1; thescript='_'; gotsbsp=1;}
		  else if (tolower(ch)=='q') return 0; 
		  else if (ch!='\n') printf("please enter u or l\n");
		  while (ch!='\n') scanf("%c",&ch);
		  if (gotsbsp==0) scanf("%c",&ch);}


	        gotkind = 0;
		printf("Gamma-0 or Gamma-1? (0 or 1)\n");
		scanf("%c",&ch);

		while (gotkind==0){
		  if (tolower(ch)=='0') {kind=0; thetype='0'; gotkind=1;}
		  else if (tolower(ch)=='1') {kind=1; thetype='1'; gotkind=1;}
		  else if (tolower(ch)=='q') return 0; 
		  else if (ch!='\n') printf("please enter 0 or 1, or q to quit\n");		  
		  while (ch!='\n') scanf("%c",&ch);
		  if (gotkind==0)  scanf("%c",&ch);}


		printf("enter ps filename (including .ps part) \n");
		scanf("%s",GammaNps);

		gotName=0;
		while (gotName==0){
		  output=fopen(GammaNps,"r");		
		  if (fopen(GammaNps,"r")!=NULL){
		    fclose(output);
		    printf("%s already exists - ",GammaNps);
		    printf("are you sure you want to over write it?\n");
		    printf("enter n for no, and anything else to ");
		    printf("continue to write the file\n");
		    scanf("%c",&ch);
		    while (ch!='\n') scanf("%c",&ch);
		    scanf("%c",&ch);
		    if (ch!='n') gotName=1;
		    else {
		      printf("enter ps filename (including .ps part) \n");
		      scanf("%s",GammaNps);
		    }
		  }
		  else gotName=1;
		}

                printf("\n");
                printf("I'm now writing to the file %s\n",GammaNps);
                printf("when I've finished I'll ask you to tell me another N.\n");
                printf("\n");


		/*Output for the start of the PostScript file. */
		output = fopen(GammaNps,"w"); 
					 if (!output) {
			fprintf(stderr, "I couldn't open the output file!\n");
			return 1;
		}


		/* draw base line */

		for(i=0;i<max;++i){
			for(j=0;j<max;++j){
				for(k=0;k<2;++k){
					for(l=0;l<2;++l){
						done[i][j][k][l]=0;
					}
				}
			}
		}
		topmatter(N,sbsp);

		closed_line(-1*infinity,0,infinity*N,0);
		
		
		make_title(thescript,thetype,N);

		proj_eq(N, kind);



	/*	for(i=0;i<N;++i){
			fprintf(outputtxt,"next \n");
			for(j=0;j<N;++j){
				fprintf(outputtxt,"[%d, %d]",proj_st[i][j][0],proj_st[i][j][1]);
			}
		}

       */


		if (sbsp==0){
			matinit(m1,1,0,0,1);
			findfund(N,m1,done,0,sbsp);
		}
		if (sbsp==1){
			matinit(m1,1,-N/3,0,1);
			findfund(N,m1,done,N/3,sbsp);
		}

		temp=0;



		shorten(N,done);



		for(i=0;i<N;++i){


			for(j=0;j<N;++j){
				
				if(done[i][j][0][0]*done[i][j][1][1]-done[i][j][1][0]*done[i][j][0][1]==1){
					if(sbsp==0){
						transtri(done[i][j][0][0],done[i][j][0][1],done[i][j][1][0],done[i][j][1][1],1,0,0);

						temp=temp+1;
					}
					if(sbsp==1) {
						mm[0][0]=done[i][j][0][0];mm[0][1]=done[i][j][0][1];mm[1][0]=done[i][j][1][0];mm[1][1]=done[i][j][1][1];
						matmult(mt,S,mm);
						transtri(mt[0][0],mt[0][1],mt[1][0],mt[1][1],1,0,0);
					}


				}
			}
		}

/*		printf("\n index: %d \n",temp);  */

		fprintf(output, "\n\n showpage\n");
		fclose(output);
                printf("done!\n");
		}
		
		
		return 1;
}

int
main()
{

	int             i;


		for (i=0;i<5;++i){
		if(dopsfile()==0) return 0;
	}

			return 0;

}