/* cellsnap.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 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 *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 = 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);
    } /*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 = 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;
                 /* 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 ********************/



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