/* bitwist.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
 * bitwist construction.
 *
 * Copyright (c) 1998, 2000, 2005 James W. Cannon, William J. Floyd,
 * and Walter R. Parry
 *
 * This work was supported in part by National Science Foundation
 * grants DMS-9803868, DMS-9971783, and DMS-0203902.
 *
 * 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 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 *faceacrossedge;    /* Array: gives tileid adjacent across
                * each edge of each tile.*/
               /* Since each tile has several edges, there are more edges
                * than tiles.*/
int *edgeacrossedge;    /* 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_and_edges_in_correct_order:
   <id tile_# edge_# tile_# edge_# ... edge_#>
   <id tile_# edge_# tile_# edge_# ... edge_#>
   ...
   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((faceacrossedge =(int *)malloc(noofedges * sizeof(int))) == NULL)
{   fprintf(stderr,"allocation error - aborting\n");
    exit(EXIT_FAILURE);
}
if((edgeacrossedge = (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;
        faceacrossedge[m] = l;
        fscanf(fp," %d",&l);
        m = facebegin[i]+k;
        edgeacrossedge[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 *cycletwist;  /* Array: gives the twist parameter 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 edgeshift; /* Gives the shift in new edges for a face pairing.*/
    int *oldedgemultiplier; /* Array: gives the multiplier for
                             * each old edge.*/
    int *oldedgetwist; /* Array: gives the twist parameter for
                             * each old edge.*/
    int *oldedgebegin; /* Array: gives the first new tile number for an
                        * old edge. */
    int *oldedgesticker; /* Array: gives 0 if there is no sticker before
                          * the old edge and 2 if there is a sticker. */
    int *newfacebegin;   /* Array: gives first new edge number of a face.*/
    int *newfaceend;     /* array: gives last new edge number of a face.*/

/* 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((cycletwist = (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);
}
if((oldedgetwist = (int *)malloc(noofedges * sizeof(int))) == NULL)
{   fprintf(stderr,"allocation error 5 - aborting\n");
    exit(EXIT_FAILURE);
}
if((oldedgesticker = (int *)malloc(noofedges * sizeof(int))) == NULL)
{   fprintf(stderr,"allocation error 5 - aborting\n");
    exit(EXIT_FAILURE);
}
if((newfacebegin = (int *)malloc(nooftiles * sizeof(int))) == NULL)
{   fprintf(stderr,"allocation error 5 - aborting\n");
    exit(EXIT_FAILURE);
}
if((newfaceend = (int *)malloc(nooftiles * 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 = faceacrossedge[facebegin[t2]+e2];
                e1 = edgeacrossedge[facebegin[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);
    fprintf(stderr,"For each cycle, each edge in the cycle will be\n");
    fprintf(stderr,"twisted in the positive direction or in the\n");
    fprintf(stderr,"negative direction. Enter 1 if it should be\n");
    fprintf(stderr,"twisted in the positive direction and -1 if it\n");
    fprintf(stderr,"should be twisted in the negative direction.\n");
}/*endfor i*/
for(i=0;i < cycleno+1; i++) {
    fprintf(stderr,"Enter the twist parameter for cycle %d\n",i+1);
    fgets(s,sizeof(s),stdin);
    s[strlen(s)-1]=s[strlen(s)];
    cycletwist[i] = atoi(s);
} /*end fori*/

/* Give the twist parameter of each old edge. */
for( i = 0; i < noofedges ; i++){
    oldedgetwist[i] = cycletwist[cycleofedge[i]];
}/*endfor */

/* Compute whether there is a sticker before each old edge. */
for( t = 0; t < nooftiles ; t++){ /* each tile */
    for( e = 0; e < facesize[t] ; e++){ /* each edge of that tile */
        i = facebegin[t]+e;
        if((oldedgetwist[i]==1) && 
               (oldedgetwist[facebegin[t]+decr(t,e)]==-1)){
fprintf(stderr,"sticker before edge %d\n",i);
            oldedgesticker[i] = 2;
        } else {
            oldedgesticker[i] = 0;
        }/*endif*/
    }/*endfor e*/
}/*enddor t*/

/* 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] = oldedgesticker[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]+
    oldedgesticker[i];
} /*end fori*/

/* Find the first and last new edge number of each tile. */
for(i=0;i < nooftiles; i++) {
    j = facebegin[i];
    newfacebegin[i] = oldedgebegin[j] - oldedgesticker[j];
    j = facebegin[i] + facesize[i] - 1;
    newfaceend[i] = oldedgebegin[j] + oldedgemultiplier[j] -1;
fprintf(stderr,"t = %d,newfacebegin = %d,newfaceend = %d\n",
i,newfacebegin[i], newfaceend[i]);
}/*endfor*/

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);
}

for( t = 0; t < nooftiles ; t++){ /* each tile */
    /* Compute newtile0, the image tile of the face the
    * triangle corrsponding to the current tile is glued to.*/
    j = facebegin[t];
    k = oldedgebegin[j];
    t1 = facepair[t]; /* the face t is glued to*/
    e1 = edgepair[t]; /* e1 is the image of edge 0*/
    l = facebegin[t1] + e1;
    edgeshift = oldedgebegin[l] + oldedgemultiplier[l] -1 -oldedgetwist[l];
    if ( edgeshift > newfaceend[t1]) {
        edgeshift = newfacebegin[t1];
    } /* modular arithmetic*/
    if ( edgeshift < newfacebegin[t1]) {
        edgeshift = newfaceend[t1];
    } /* modular arithmetic*/
    for(i=newfacebegin[t]; i < newfaceend[t]+1; i++){ /*each new edge*/
        newtile0[i] = edgeshift - i + newfacebegin[t];
        if (newtile0[i] < newfacebegin[t1]) {
        newtile0[i]=newtile0[i]-newfacebegin[t1]+newfaceend[t1]+1;
    } /* modular arithmetic*/
    }/*endfor i*/
} /*end for t*/

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];
        /*Compute newtile1 for stickers.*/
        if(oldedgesticker[j]==2){
fprintf(stderr,"t=%d,e=%d,j=%d,k=%d\n",t,e,j,k);
            newtile1[k-2] = k-1;
            newtile1[k-1] = k-2;
        }/*endif*/
        for(i=0; i < oldedgemultiplier[j]; i++) { /* each subinterval*/
            /*endif oldedgesticker*/
            /* Compute newtile1, the tile of the face across from the
             * subinterval.*/
            t2 = faceacrossedge[facebegin[t]+e]; /*the tile across t from e*/
            e2 = edgeacrossedge[facebegin[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;
        } /*end for i*/
    } /*end for e*/
} /*end for t*/

/* 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(i=newfacebegin[t]; i < newfaceend[t]+1; i++){ /*each new edge*/
        /* Compute newtile2, the tile of the face to the left.*/
        newtile2[i] = i-1;
        if(i==newfacebegin[t]) {newtile2[i] =newfaceend[t];}
        newtile3[i] = i+1;
        /* Compute newtile3, the tile of the face to the right.*/
        if(i==newfaceend[t]){ newtile3[i] =newfacebegin[t];}
    }/*endfor i*/
}/* endfor t*/

} /*end Computetwistdata */

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

/* 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 ***********************/