/* partsnap.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 cusped 3-manifold obtained from it by the
 * modified twist construction with a single edge collapse.
 *
 * Copyright (c) 1999, 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>
#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 triangulation.*/
int noofnewtiles; /* Acts as size of each of next three arrays.*/
int *newfacebegin;/* Array: gives first edgeid associated with each tile.*/
int *newedges;/*Array: gives tileid adjacent across each edge of each tile.*/
              /* Since each tile has several edges, there are more edges than
               * tiles.*/

/* Variables for the adjacencies of the new tiles.*/
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. */

/* Variables for the chosen face and its inverse.*/
int face1,face2; /*Face1 is the chosen face and face2 is its inverse.*/

main(argc, argv) int argc; char *argv[];
{
            Readtilingforsnap();
            Computetwistdata();
            Writeforsnappea();
} /*end main() */
/******************* end main() ********************/

/* incr */
/* Description:
 * Return the next edge number of a tile in the clockwise direction.
 */
int incr(t,e)
    int t,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(t,e)
    int t,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.

 */

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

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

FILE *fp;


fprintf(stderr," Read which tiling file?\n");
gets(s);
if ((fp = fopen(s,"r")) == NULL)
{
  fprintf(stderr," cannot open file \n");
  exit();
}

/* 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(0);
    }
    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(0);
    }

  /* 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(0);
    }

    /* 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(0);
    }
    if((edgepair =
       (int *)malloc(nooftiles * sizeof(int))) ==
       NULL)
    {   fprintf(stderr,"allocation error - aborting\n");
        exit(0);
    }

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 Jim's twist construction.
 */

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

    /* standard input variable and functions */
    char s[256];
    extern char *gets();
    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 *Rcyclelength;     /* Array: gives the reduced 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. */
int *Rfacesize;       /*Array: gives the new size of each face.*/
int *Rfacebegin;	    /* Array: gives the first new edge of each face.*/ 

/* Allocate space for variables. */
    if((cycleofedge =
       (int *)malloc(noofedges * sizeof(int))) == NULL)
    {   fprintf(stderr,"allocation error 5 - aborting\n"); exit(0); }
    if((cyclelength =
       (int *)malloc(noofedges * sizeof(int))) == NULL)
    {   fprintf(stderr,"allocation error 5 - aborting\n"); exit(0); }
    if((Rcyclelength =
       (int *)malloc(noofedges * sizeof(int))) == NULL)
    {   fprintf(stderr,"allocation error 5 - aborting\n"); exit(0); }
    if((cyclemultiplier =
       (int *)malloc(noofedges * sizeof(int))) == NULL)
    {   fprintf(stderr,"allocation error 5 - aborting\n"); exit(0); }
        if((newcyclebegin =
       (int *)malloc(noofedges * sizeof(int))) == NULL)
    {   fprintf(stderr,"allocation error 5 - aborting\n"); exit(0); }
        if((tileincycle =
       (int *)malloc(noofedges * sizeof(int))) == NULL)
    {   fprintf(stderr,"allocation error 5 - aborting\n"); exit(0); }
        if((Rfacesize =
       (int *)malloc(nooftiles * sizeof(int))) == NULL)
    {   fprintf(stderr,"allocation error 5 - aborting\n"); exit(0); }
        if((Rfacebegin =
       (int *)malloc(nooftiles * sizeof(int))) == NULL)
    {   fprintf(stderr,"allocation error 5 - aborting\n"); exit(0); }
        if((edgeincycle =
       (int *)malloc(noofedges * sizeof(int))) == NULL)
    {   fprintf(stderr,"allocation error 5 - aborting\n"); exit(0); }
        if((oldedgebegin =
       (int *)malloc(noofedges * sizeof(int))) == NULL)
    {   fprintf(stderr,"allocation error 5 - aborting\n"); exit(0); }
        if((oldedgemultiplier =
       (int *)malloc(noofedges * sizeof(int))) == NULL)
    {   fprintf(stderr,"allocation error 5 - aborting\n"); exit(0); }
/* 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);
        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,"The tiles are numbered 0-%d. Enter the number\n",nooftiles-1);
    fprintf(stderr,"of the tile you wish to use for the construction.\n");
    gets(s);
    face1=atoi(s);
    face2=facepair[face1];
    fprintf(stderr,"\n");
    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 that do not involve the chosen tile\n");
    fprintf(stderr,"or its paired tile. 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);
        gets(s);
        cyclemultiplier[i] = atoi(s);
    } /*end fori*/

    /* Compute the array Rcyclelength.*/
    for(i=0;i < cycleno+1;i++){
      Rcyclelength[i]=0;
      for(j=0;j < cyclelength[i];j++){
        if((tileincycle[newcyclebegin[i]+j] != face1) &&
        (tileincycle[newcyclebegin[i]+j] != face2)){Rcyclelength[i]++;}
      } /*end forj*/
      fprintf(stderr,"Rcyclelength[%d] = %d.\n",i,Rcyclelength[i]);
    } /* 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] = Rcyclelength[0] * cyclemultiplier[cycleofedge[0]];
    for(i=1;i < noofedges; i++) {
        oldedgemultiplier[i] =
          Rcyclelength[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(0); }
        if((newtile1 =
       (int *)malloc(noofnewtiles * sizeof(int))) == NULL)
    {   fprintf(stderr,"allocation error t1 - aborting\n"); exit(0); }
        if((newtile2 =
       (int *)malloc(noofnewtiles * sizeof(int))) == NULL)
    {   fprintf(stderr,"allocation error t2 - aborting\n"); exit(0); }
        if((newtile3 =
       (int *)malloc(noofnewtiles * sizeof(int))) == NULL)
    {   fprintf(stderr,"allocation error t3 - aborting\n"); exit(0); }

    /* Compute the new size of each face. */
    for(t=0; t < nooftiles; t++){ /* each tile */
        Rfacesize[t]=0;
        for( e = 0; e < facesize[t] ; e++){ /* each edge of that tile */
            Rfacesize[t]=Rfacesize[t]+oldedgemultiplier[facebegin[t]+e];
        } /* end fore */
    fprintf(stderr,"Rfacesize[%d] = %d.\n",t,Rfacesize[t]);
    } /* end fort */

    /* Compute the first new edge of each face. */
    Rfacebegin[0]=0;
    for(t=1; t < nooftiles; t++){ /* each tile */
        Rfacebegin[t]=Rfacebegin[t-1]+Rfacesize[t-1];
    } /* end fort */

    /* 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;
                }
		    if ((t == face1) || (t == face2)) {
			newtile0[k+i]++;
                   if(newtile0[k+i]==Rfacebegin[t1]+Rfacesize[t1]) {
                      newtile0[k+i] = Rfacebegin[t1];} /*end ifn */
                } /* end ift */
                /* 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);
                    while(oldedgemultiplier[facebegin[t]+e2]==0){ 
                        e2 = decr(t,e2);}
                    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);
                    while(oldedgemultiplier[facebegin[t]+e2]==0){
                        e2 = incr(t,e2);}
                    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(t1,t2,e)
    int t1,t2,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 */


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

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

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

fprintf(stderr," Write which SnapPea input file?\n");
gets(s);
if ((fp = fopen(s,"w")) == NULL)
{
  fprintf(stderr," cannot open file \n");
  exit();
      }
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 ***********************/