/* pairsnap.c */
/* Description:
 * Reads a tiling of a 2-sphere (with face-pairing data)
 * from a file and writes the associated input file for SnapPea
 * for a closed 3-manifold obtained from it by the 
 * twist construction. 
 *
 * Copyright (c) 1998, 2000  James W. Cannon, William J. Floyd and
 * Walter R. Parry
 * 
 * This work was supported in part by National Science Foundation 
 * grants DMS-9803868 and DMS-9971783.
 *
 * This is free software and may be freely copied, modified, and
 * redistributed as noted below.
 *      The copyright notice must remain intact.
 *      Modifications must include a notice giving the reason
 *      for the modification and including name and date.
 *      This is provided "as is" and comes with no warranty or
 *      guarantee of fitness. Comments, bugs, and suggestions may be sent
 *      to floyd@math.vt.edu.
 *
 */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

void Readtilingforsnap();
void Computetwistdata();
void Writeforsnappea();
int edgefromtiletotile(int t1,int t2,int e);
int incr(int t, int e);
int decr(int t, int e);







#define MAXVALENCE 50
        /* The assumed maximal number of adjacencies to a tile. */
        
/* Variables for the old tiling.*/
int noofedges; /* Acts as size of *edges. */
int nooftiles; /* Acts as size of each of next three arrays.*/
int *facesize; /* Array: gives no of edges of each tile.*/
int *facebegin;/* Array: gives first edgeid associated with each tile.*/
int *edges;    /* Array: gives tileid adjacent across 
                * each edge of each tile.*/
               /* Since each tile has several edges, there are more edges 
                * than tiles.*/

/* Variable for the face-pairings.*/
int *facepair; /* Array: gives face(tile) pairing for each tile.*/
int *edgepair; /* Array: gives edge of facepair[i] identified with
                * edge 0 of face i.*/

/* Variables for the adjacencies of the new tiles.*/               
int noofnewtiles; /* Acts as size of each of next four arrays.*/
int *newtile0; /* Array: the new tile across from face 0 of a new tile. */
int *newtile1; /* Array: the new tile across from face 1 of a new tile. */
int *newtile2; /* Array: the new tile across from face 2 of a new tile. */
int *newtile3; /* Array: the new tile across from face 3 of a new tile. */




int main()
{
            Readtilingforsnap();
            Computetwistdata();
            Writeforsnappea(); 
} /*end main() */
/******************* end main() ********************/

/* incr */
/* Description:
 * Return the next edge number of a tile in the clockwise direction.
 */
int incr(int t,int e)
{   /* loop variables */
    int i,j,k;

    i = e+1;
    if( i == facesize[t]
    ){
        i = 0;
    }/*endif */

    return i;

} /*end incr */

/******************* end incr ********************/

/* decr */
/* Description:
 * Return the previous edge number of a tile in the clockwise direction.
 */
int decr(int t,int e)
{   /* loop variables */
    int i,j,k;

    i = e-1;
    if( i == -1
    ){
        i = facesize[t] - 1;
    }/*endif*/

    return i;
} /* end decr */

/******************* end decr ********************/


/* Readtilingforsnap */
/* Description:

 * This function reads in a tiling of a 2-sphere. The format is
 * similar to the format for subdivide.c and for tilepack.c, except
 * now there is no type information. 

 * The text portions of the format are 1-string comments that must
 * appear in the file.
 * The comments are read and discarded by fscanf as a single string.
 * The content of the comment lines is only to help those writing input
 * files.
 * Here is the format:

   Number_of_tiles_including_the_four_standard_ends:
   <number>
   Size_of_each_of_these_tiles:
   <numbers (no tile-number id is to be included)>
   Tile-ids_--_Adjacent_tiles_in_correct_order:
   <id number ... number>
   <id number ... number>
   ...
   Tile-id_--_image_tile_id_and_edge_no_of_edge_0:
   <id number id number edge number>
   <id number id number edge number>
   ... 
   (Any number of comments or extra data not to be used by the program>

  * The program does a fair amount of processing of these rules in order
  * to have indexes into the information.

 */

void Readtilingforsnap()
{   /* loop variables */
    int i,j,k,l,m;

    /* standard input variable and functions */
    char s[256];
    extern char *fgets();
    extern int atoi();

FILE *fp; 


fprintf(stderr," Read which tiling file?\n");
fgets(s,sizeof(s),stdin);
s[strlen(s)-1]=s[strlen(s)];
if ((fp = fopen(s,"r")) == NULL)
{
  fprintf(stderr," cannot open file \n");
  exit(EXIT_FAILURE);
}

/* Read nooftiles.*/
fscanf(fp," %s",s);/* Skip comment string.*/
fscanf(fp," %d",&nooftiles);

/* Read the sizes of the tiles( *facesize).*/
fscanf(fp," %s",s);/* Skip comment string. */
    if((facesize =
       (int *)malloc(nooftiles * sizeof(int))) ==
       NULL)
    {   fprintf(stderr,"allocation error - aborting\n");
        exit(EXIT_FAILURE);
    }
    for( i = 0; i < nooftiles ; i++){
        fscanf(fp," %d",&l);
        facesize[i] = l;
    }/*endfor */

/* Calculate the array *facebegin.*/
  /* Allocate space for the array. */

    if((facebegin =
       (int *)malloc(nooftiles * sizeof(int))) ==
       NULL)
    {   fprintf(stderr,"allocation error - aborting\n");
        exit(EXIT_FAILURE);
    }

  /* Calculate the array. */

    /* Calculate facebegin.*/
    facebegin[0] = 0;
    for( i = 0; i < nooftiles -1  ; i++){
        facebegin[i+1] = facebegin[i]+facesize[i];
    }/*endfor */
    noofedges = facebegin[nooftiles - 1] + facesize[nooftiles - 1];

/* For each tile read the id and the adjacent tiles.*/
fscanf(fp," %s",s);/* Skip comment string.*/

    /* Allocate space for edges.*/

    if((edges =
       (int *)malloc(noofedges * sizeof(int))) ==
       NULL)
    {   fprintf(stderr,"allocation error - aborting\n");
        exit(EXIT_FAILURE);
    }

    /* Read the edges. */
    for( i = 0; i < nooftiles ; i++){
        /* Read and discard id.*/
        fscanf(fp," %d",&l);

        /* Read the tiles adjacent to the given tile.*/
        j = facesize[i];
        for( k = 0; k < j ; k++){
        fscanf(fp," %d",&l);
        m = facebegin[i]+k;
        edges[m] = l;
        }/*endfor */
    }/*endfor */
    
    /* Allocate space for the face-pairing data. */

    if((facepair =
       (int *)malloc(nooftiles * sizeof(int))) ==
       NULL)
    {   fprintf(stderr,"allocation error - aborting\n");
        exit(EXIT_FAILURE);
    }                                                        
    if((edgepair =
       (int *)malloc(nooftiles * sizeof(int))) ==
       NULL)
    {   fprintf(stderr,"allocation error - aborting\n");
        exit(EXIT_FAILURE);
    }                     

fscanf(fp," %s",s);/* Skip comment string.*/       
                            
    /* Read the face-pariing data. */
    for( i = 0; i < nooftiles ; i++){
        /* Read and discard id.*/
        fscanf(fp," %d",&l);
        fscanf(fp," %d",&l); /* Read face-pairing */
        facepair[i] = l;
        fscanf(fp," %d",&l); /* Read the image of edge 0 under 
                             * the face-pairing. */
        edgepair[i] = l;
    } /* endfor i */
   
fclose(fp);

} /*end Readtilingforvertex */

/******************* end Readtilingforvertex ********************/

/* Computetwistdata */
/* Description:
 * From the data of a tiling of a 2-sphere and face-pairing data, this
 * programs prompts the user for a multiplier for each edge cycle, and
 * then computes the face-pairings for a triangulation of the resulting
   * closed 3-manifold obtained from the twist construction.
 */

void Computetwistdata()
{   /* loop variables */
    int i,j,k,l,m;

    /* standard input variable and functions */
    char s[256];
    extern char *fgets();
    extern int atoi();

    /* tile and edge loop variables */
    int t,e,t1,e1,t2,e2,e3,v;
    
    /* Variables to describe the edge cycles.*/
int cycleno;          /* The number of edge cycles. */
int *cycleofedge;     /* Array: gives the edge cycle of each edge. */
int *cyclelength;     /* Array: gives the length of each cycle. */
int *cyclemultiplier; /* Array: gives the multiplier of each cycle. */
int *newcyclebegin;   /* Array: gives the beginning in the next two
                       * arrays of each cycle. */
int *tileincycle;     /* Array: gives the tile numbers of 
                       * the edges in a cycle.*/
int *edgeincycle;     /* Array: gives the edge numbers of 
                       * the edges in a cycle.*/
    /* Variables to describe the new tiles.*/ 
int *oldedgemultiplier; /* Array: gives the multiplier for each old edge.*/
int *oldedgebegin;      /* Array: gives the first new tile number for an 
                         * old edge. */
                         
    
/* Allocate space for variables. */
    if((cycleofedge =
       (int *)malloc(noofedges * sizeof(int))) == NULL)
    {   fprintf(stderr,"allocation error 5 - aborting\n"); exit(EXIT_FAILURE); }
    if((cyclelength =
       (int *)malloc(noofedges * sizeof(int))) == NULL)
    {   fprintf(stderr,"allocation error 5 - aborting\n"); exit(EXIT_FAILURE); }
    if((cyclemultiplier =
       (int *)malloc(noofedges * sizeof(int))) == NULL)
    {   fprintf(stderr,"allocation error 5 - aborting\n"); exit(EXIT_FAILURE); }
        if((newcyclebegin =
       (int *)malloc(noofedges * sizeof(int))) == NULL)
    {   fprintf(stderr,"allocation error 5 - aborting\n"); exit(EXIT_FAILURE); }
        if((tileincycle =
       (int *)malloc(noofedges * sizeof(int))) == NULL)
    {   fprintf(stderr,"allocation error 5 - aborting\n"); exit(EXIT_FAILURE); }
        if((edgeincycle =
       (int *)malloc(noofedges * sizeof(int))) == NULL)
    {   fprintf(stderr,"allocation error 5 - aborting\n"); exit(EXIT_FAILURE); }
        if((oldedgebegin =
       (int *)malloc(noofedges * sizeof(int))) == NULL)
    {   fprintf(stderr,"allocation error 5 - aborting\n"); exit(EXIT_FAILURE); }
        if((oldedgemultiplier =
       (int *)malloc(noofedges * sizeof(int))) == NULL)
    {   fprintf(stderr,"allocation error 5 - aborting\n"); exit(EXIT_FAILURE); }    
/* Initialize the edgecycle array. */
    for( i = 0; i < noofedges ; i++){
        cycleofedge[i] = -1;
    }/*endfor */

    cycleno = -1;
    
/* Process each edge provided
 * that the edge has not already been processed.*/
    for( t = 0; t < nooftiles ; t++){ /* each tile */
        for( e = 0; e < facesize[t] ; e++){ /* each edge of that tile */
            /* Edge not already used?  i == -1*/
            i = cycleofedge[facebegin[t]+e];
            if( i == -1
            ){
            	cycleno++;
            	j=0;
            	if (cycleno == 0
            	) { newcyclebegin[cycleno] = 0;
                k=0;
            	}
            	else
            	{
            		newcyclebegin[cycleno] = 
             		newcyclebegin[cycleno - 1] + cyclelength[cycleno - 1];
            		k = newcyclebegin[cycleno];            	
            	} /* end if cycleno */
            	v=-1;
            	t2=t;
            	e2=e;
            	while(v == -1)
            	{
                cycleofedge[facebegin[t2]+e2] = cycleno;
                tileincycle[k] = t2;
                edgeincycle[k] = e2;
            	    t1 = edges[facebegin[t2]+e2];
            	    e1 = edgefromtiletotile(t1,t2,e2);
            	    cycleofedge[facebegin[t1]+e1] = cycleno;
            	    t2 = facepair[t1];
            	    e2 = edgepair[t1] - e1;
            	    if (e2 < 0) {
            		    e2 = e2 + facesize[t2];
            	    } /* if e2 */
             	    k++;
            	    v = cycleofedge[facebegin[t2]+e2];        	
            	} /*endwhile v*/
            	cyclelength[cycleno] = k - newcyclebegin[cycleno];
            }/*endif i */
        }/*endfor e = edge no*/
    }/*endfor t = tile no*/
    
    /* Prompt the user for the multipliers. */
    fprintf(stderr,"The number of cycles is %d\n",cycleno+1);
    fprintf(stderr,"\n");
    for(i=0;i < cycleno+1; i++) {
    	fprintf(stderr,"The (tile,edge) pairs in cycle %d are\n",i+1);
    	for(j=0;j < cyclelength[i]; j++) {
    		fprintf(stderr,"(%d,%d) ",
    		tileincycle[newcyclebegin[i]+j],
    		edgeincycle[newcyclebegin[i]+j]);
    	} /* end forj*/
        fprintf(stderr,"\n");
    } /*end fori*/
    fprintf(stderr,"For each cycle, each edge in the cycle will be\n");
    fprintf(stderr,"subdivided into a multiple times the numer of edges\n");
    fprintf(stderr,"in that cycle. Enter the multiplier to use\n"); 
    fprintf(stderr,"for each cycle.\n");
    for(i=0;i < cycleno+1; i++) {
        fprintf(stderr,"Enter the multiplier for cycle %d\n",i+1);
        fgets(s,sizeof(s),stdin);
        s[strlen(s)-1]=s[strlen(s)];
    	cyclemultiplier[i] = atoi(s);
    } /*end fori*/
    
    /* Find the number of edges each old edge is subdivided into and the index
     * number of the first new tile meeting each old edge.*/   
    oldedgebegin[0] = 0;
    oldedgemultiplier[0] = cyclelength[0] * cyclemultiplier[cycleofedge[0]];
    for(i=1;i < noofedges; i++) {
    	oldedgemultiplier[i] = 
          cyclelength[cycleofedge[i]] * cyclemultiplier[cycleofedge[i]];
    	oldedgebegin[i] = oldedgebegin[i-1] + oldedgemultiplier[i-1];
    } /*end fori*/
    noofnewtiles=oldedgebegin[noofedges-1]+oldedgemultiplier[noofedges-1];
    /*Allocate space for the tile adjacency information.*/
        if((newtile0 =
       (int *)malloc(noofnewtiles * sizeof(int))) == NULL)
    {   fprintf(stderr,"allocation error t0 - aborting\n"); exit(EXIT_FAILURE); }
        if((newtile1 =
       (int *)malloc(noofnewtiles * sizeof(int))) == NULL)
    {   fprintf(stderr,"allocation error t1 - aborting\n"); exit(EXIT_FAILURE); }
        if((newtile2 =
       (int *)malloc(noofnewtiles * sizeof(int))) == NULL)
    {   fprintf(stderr,"allocation error t2 - aborting\n"); exit(EXIT_FAILURE); }
        if((newtile3 =
       (int *)malloc(noofnewtiles * sizeof(int))) == NULL)
    {   fprintf(stderr,"allocation error t3 - aborting\n"); exit(EXIT_FAILURE); }  
    
    /* Compute the adjacencies for the new tiles. The new tiles correspond
     * to subintervals of the old edges, and are processed in order of the
     * old edges. */
    for( t = 0; t < nooftiles ; t++){ /* each tile */
        for( e = 0; e < facesize[t] ; e++){ /* each edge of that tile */
            j = facebegin[t] + e;
            k = oldedgebegin[j];
            for(i=0; i < oldedgemultiplier[j]; i++) { /* each subinterval*/
                /* Compute newtile0, the image tile of the face the 
                 * triangle corrsponding to the current tile is glued to.*/
                t1 = facepair[t]; /* the face t is glued to*/
                e1 = edgepair[t] - e; /* e1 is the image of edge e*/
                if ( e1 < 0) {
                	e1 = e1 + facesize[t1];
                } /* modular arithmetic*/
                if (i == oldedgemultiplier[j]-1) {
                	e3=decr(t1,e1);
                	l = facebegin[t1]+e3;
                	newtile0[k+i] = oldedgebegin[l]+oldedgemultiplier[l]-1;
                } else {
                	l = facebegin[t1]+e1;
                	newtile0[k+i] = oldedgebegin[l] + oldedgemultiplier[l]-2-i;
                }
                /* Compute newtile1, the tile of the face across from the
                 * subinterval.*/
                t2 = edges[facebegin[t]+e]; /* the tile across t from e*/
                e2 = edgefromtiletotile(t2,t,e); /*the edge of that tile*/
                l = facebegin[t2] +e2; /* the edge number of that edge*/
                newtile1[k+i] = oldedgebegin[l]+ oldedgemultiplier[l] -1-i;                              
                 /* Compute newtile2, the tile of the face to the left.*/
                if (i == 0) {
                	e2 = decr(t,e);
                	l = facebegin[t]+e2;
                	newtile2[k+i] = oldedgebegin[l]+ oldedgemultiplier[l]-1;
                } else {
                	newtile2[k+i] = k+i-1; /* decrease by 1 if not at far left*/
                }
                /* Compute newtile3, the tile of the face to the right.*/
                if (i == oldedgemultiplier[j]-1) {
                	e2 = incr(t,e);
                	l = facebegin[t]+e2;
                	newtile3[k+i] = oldedgebegin[l];
                } else {
                	newtile3[k+i] = k+i+1; /* increase by 1 if not at far right*/
                }                              
            } /*end for i*/
        } /*end for e*/
    } /*end for t*/
} /*end Computetwistdata */

/******************* end Computetwistdata ********************/


/*Edgefromtiletotile */
/* Description:
 * Return the edge number of tile t1 that is adjacent to tile t2
 * corresponding to edge e of tile t2.
 */
int edgefromtiletotile(int t1,int t2,int e)
{   /* loop variables */
    int i,j,k;
    /* edge variables to check whether t1 and t2 are incident in more
     * than one adjacent face. */
    int incre,decrj,eplus,jminus;
    i = facebegin[t1];
    j = 0;
    while( edges[i + j] != t2
    ){
        j++;
    }/*endwhile */
    incre = incr(t2,e);
    eplus=0;
    while (
    edges[facebegin[t2] + incre] == t1)
    {  incre = incr(t2,incre);
         eplus++;
      }
    decrj = decr(t1,j);
    jminus=0;
    while (
    edges[facebegin[t1] + decrj] == t2)
    {
        decrj = decr(t1,decrj);
        jminus++;
    }
    if (eplus == jminus) {
        k = j;
    }
    if (eplus > jminus) {
        k = j;
        for (i=0; i < (eplus - jminus); i++) {
            k = incr(t1,k);
        }
    }
    if (eplus < jminus) {
        k = j;
        for (i=0; i < (jminus - eplus); i++) {
            k = decr(t1,k);
        }
    }
    return k;

} /*end edgefromtiletotile */

/******************* end edgefromtiletotile ********************/



/* Writeforsnappea */
/* Description:
 * Write to file the face-pairings for the resulting triangulation
 * in suitable format for input to SnapPea 2.5. */

/* Here is the format.
 * % Triangulation
 *
 * name of 3-manifold
 * no_solution 0.0
 * unknown_orientability
 * CS_unknown
 *
 * 0 0
 *
 * number of tetrahedra
 *
 * numbers of tetrahedra across from each face of tetrahedra i
 * 0132 0132 0132 0132
 * -1 -1 -1 -1
 * 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 * 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 * 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 * 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 * 0.0 0.0    
 * 
 * repeat above 9 lines for each tetrahedron */


void Writeforsnappea()
{   /* loop variables */
    int i,j,k,l,m;
    int linesize;

    /* standard input variable and functions */
    char s[256];
    extern char *fgets();
    extern int atoi();

FILE *fp;  /* move this to the variable declaration section */

fprintf(stderr," Write which SnapPea input file?\n");
fgets(s,sizeof(s),stdin);
s[strlen(s)-1]=s[strlen(s)];
if ((fp = fopen(s,"w")) == NULL)
{
  fprintf(stderr," cannot open file \n");
  exit(EXIT_FAILURE);
      }
fprintf(fp,"%% Triangulation\n");
fprintf(fp,"\n");
fprintf(fp,"%s\n",s);
fprintf(fp,"no_solution 0.0\n");
fprintf(fp,"unknown_orientability\n");
fprintf(fp,"CS_unknown\n");
fprintf(fp,"\n");
fprintf(fp,"0 0\n");
fprintf(fp,"\n");
fprintf(fp,"%d\n",noofnewtiles);
fprintf(fp,"\n");
for( i = 0; i < noofnewtiles ; i++){
	 fprintf(fp,"%d %d %d %d\n",
	 newtile0[i],newtile1[i],newtile2[i],newtile3[i]);
         fprintf(fp,"0132 0132 0132 0132\n");
	 fprintf(fp,"-1 -1 -1 -1\n");
	 fprintf(fp,"0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n");
	 fprintf(fp,"0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n");
	 fprintf(fp,"0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n");
	 fprintf(fp,"0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n");
	 fprintf(fp,"0.0 0.0\n");
	 fprintf(fp,"\n");   
	}/*endfor i */
	
fclose(fp);
} /*end Writeforsnappea */

/******************* end Writeforsnappea ***********************/