/*
Copyright (C) 1998, 2000 Harald A. Helfgott

This program is free software; you can redistribute it and/or modify it 
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2, or (at your option)
any later version.

This library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

If you do not have a copy of the GNU General Public License, write to the Free
Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, 
USA. */
#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#ifdef INF
#include <gmp.h>
#endif
#include "ren.h"
#include "r250.h"

/* to compile this program run
   add mathlibrary
   gcc -DINF -I/mit/mathlibrary/include -L/mit/mathlibrary/lib -o ren ren.c r250.c -lgmp
   */

int mode;

#ifdef INF
void coup_set(couple *a, int q)
{
  if(q==0) {
    mpq_set_si(a->q,1,1);
    a->i=1;
  } else {
    mpq_set_si(a->q,q,1);
    a->i=0;
  }
}
#else
void coup_set(couple *a, num q)
/* sets e to zero, too */
{
  if(q==0) {
    a->q=1;
    a->i=1;
  } else {
    a->q = q;
    a->i = 0;
  }
}
#endif

num coup_spit(couple *a)
{
  if(a->i>0)
    return 0;
#ifdef INF
  return mpq_get_d(a->q);
#else
  return a->q;
#endif
}

void coup_e_set(couple *a, int i)
{
  a->i = i;
}

void coup_ass(couple *a, couple b)
{
#ifdef INF
  mpq_set(a->q,b.q);
#else
  a->q = b.q;
#endif
  a->i = b.i;
}

void coup_add(couple *a, couple b, couple c)
{
  if(b.i<c.i) {
    a->i=b.i;
#ifdef INF
    mpq_set(a->q,b.q);
#else
    a->q=b.q;
#endif
  } else
  if(c.i<b.i) {
    a->i=c.i;
#ifdef INF
    mpq_set(a->q,c.q);
#else
    a->q=c.q;
#endif
  } else
  if(c.i==b.i) {
    a->i=b.i;
#ifdef INF
    mpq_add(a->q,b.q,c.q);
#else
    a->q=b.q+c.q;
#endif
  }
}


void coup_sub(couple *a, couple b, couple c)
{
  if(b.i<c.i) {
    a->i=b.i;
#ifdef INF
    mpq_set(a->q,b.q);
#else
    a->q=b.q;
#endif
  } else
  if(c.i<b.i) {
    a->i=c.i;
#ifdef INF
    mpq_neg(a->q,c.q);
#else
    a->q = -c.q;
#endif
  } else
  if(c.i==b.i) {
    a->i=b.i;
#ifdef INF
    mpq_sub(a->q,b.q,c.q);
#else
    a->q=b.q-c.q;
#endif
  }
}

void coup_mul(couple *a, couple b, couple c)
{
#ifdef INF
  mpq_mul(a->q,b.q,c.q);
#else
  a->q = b.q*c.q;
#endif
  a->i = b.i+c.i;
}

void coup_div(couple *a, couple b, couple c)
{
#ifdef INF
  mpq_div(a->q,b.q,c.q);
#else
  if(c.q==0)
    fprintf(stderr,"Division by zero\n");
  a->q = b.q/c.q;
#endif
  a->i = b.i-c.i;
}

void coup_init(couple *a)
{
#ifdef INF
  mpq_init(a->q);
#else
  a->q = 0; 
#endif
  a->i = 0;
}

void coup_kill(couple *a)
{
#ifdef INF
  mpq_clear(a->q);
#endif
}

size_t coup_in(FILE *stream, int base, couple *op)
#ifdef INF
{
  mpz_t rop1,rop2;
  int c;
  int len;

  mpz_init(rop1);
  mpz_init(rop2);
  len=mpz_inp_str(rop1, stream, base);
  if(mpz_sgn(rop1)==0) {
    mpq_set_si(op->q,1,1);
    op->i = 1;
    c=getc(stream);
    if(c!='/') 
      ungetc(c,stream);
    else
      len+=mpz_inp_str(rop2, stream, base)+1;
    mpz_clear(rop1);
    mpz_clear(rop2);
    return len;
  } else {
    c=getc(stream);
    if(c!='/') {
      ungetc(c,stream);
      mpz_set_si(rop2,1);
      mpz_set(mpq_numref(op->q),rop1);
      mpz_set(mpq_denref(op->q),rop2);
      mpq_canonicalize(op->q);
      mpz_clear(rop1);
      mpz_clear(rop2);
      return len;
    } else {
      len++;
      len+=mpz_inp_str(rop2, stream, base);
      mpz_set(mpq_numref(op->q),rop1);
      mpz_set(mpq_denref(op->q),rop2);
      mpq_canonicalize(op->q);
      mpz_clear(rop1);
      mpz_clear(rop2);
      return len;
    }
    op->i=0;
  }
}
#else
{
  size_t res;
  float q;

  res=fscanf(stream,"%f",&q);
  coup_set(op,q);
  return res;
}
#endif

size_t coup_out(FILE *stream, int base, couple op)
#ifdef INF
{
  int len;
  len=mpz_out_str(stream,base,mpq_numref(op.q));
  if(mpz_cmp_si(mpq_denref(op.q),1)) {
    len += fprintf(stream,"/");
    len += mpz_out_str(stream,base,mpq_denref(op.q));
  }
  if(op.i>0) 
    len += fprintf(stream," e^%d",op.i);
  if(op.i<0) 
    len += fprintf(stream," e^(%d)",op.i);

  return len;
}
#else
     /* the program feels free not to use the base */
{
  size_t res;

  /* don't print number * e^x if x is positive because matlab doesn't
     like it...just print 0 instead */
  if (op.i > 0)
    res = fprintf(stream,"0");
  else {
    res = fprintf(stream,"%g",op.q);
    if(op.i<0) 
      res+=fprintf(stream," e^(%d)",op.i);
    if(op.i>0)
      res+=fprintf(stream," e^%d",op.i);
  }
  return res;
}
#endif

void matr_in(FILE *stream, int base, matr *mat)
     /* matrix assumed to be allocated and initialized */
{
  int i,j;

  for(i=0; i<mat->side; i++)
    for(j=0; j<mat->side; j++) 
      coup_in(stream,base,el(*mat,i,j));
}

void chatr_out(FILE *stream, chatr mat)
{
  int i,j,k;

  for(i=0; i<mat.side; i++) {
    for(j=0; j<mat.side; j++) 
      fprintf(stream,"%2d",*(chel(mat,i,j)));
    fprintf(stream,"\n");
  }
}

void matr_out(FILE *stream, int base, matr mat, int offset)
{
  int i,j,k;

  for(i=0; i<mat.side; i++) {
    for(j=0; j<mat.side; j++) {
      k = coup_out(stream,base,*(el(mat,i,j))); 
      if(k<offset)
	for(; k<offset; k++)
	  putc(' ',stream);
      else
	putc(' ',stream);
    }
    putc('\n',stream);
  }
}

void chatr_create(chatr *mat, int memside)
{
  int i,j;

  if(!(mat->a = (char *) calloc(memside*memside,sizeof(char))))
    crap_out();
  mat->memside = mat->side = memside;
  mat->i0 = mat->j0 = 0;
}

void chatr_kill(chatr *mat, int memside)
{
  free(mat->a);
}

void matr_create(matr *mat, int memside)
{
  int i,j;

  if(!(mat->a = (couple *) malloc(sizeof(couple)*memside*memside)))
    crap_out();
  mat->memside = mat->side = memside;
  mat->i0 = mat->j0 = 0;
  for(i=0; i<memside; i++) {
    for(j=0; j<memside; j++)
      coup_init(el(*mat,i,j));
  }
}

void matr_cpy(matr *mat, matr mat2)
{
  int i,len;

  if(mat->memside != mat2.memside) {
    fprintf(stderr,"Illegal use of matr_cpy %d %d",
	    mat->memside,mat2.memside);
    crap_out();
  }
  mat->side    = mat2.side;
  mat->i0      = mat2.i0;
  mat->j0      = mat2.j0;
#ifdef INF
  len = mat2.memside*mat2.memside;
  for(i=0; i<len; i++)
    coup_ass(mat->a+i,mat2.a[i]);
#else
  memcpy(mat->a,mat2.a,sizeof(couple)*mat2.memside*mat2.memside);
#endif
}

void matr_kill(matr *mat)
{
  int i,j;

  mat->i0=mat->j0=0;
  mat->side=mat->memside;
  for(i=0; i<mat->memside; i++)
    for(j=0; j<mat->memside; j++)
      coup_kill(el(*mat,i,j));

  free(mat->a);
}

void crap_out()
{
  fprintf(stderr,"PANIC! Out of memory, or some other stupid error.\n");
  abort();
}

void grrand(chatr *neu, matr weight)
{
  int n;
  int i,j;
  couple temp[10];

  for(i=0; i<10; i++)
    coup_init(temp+i);
  n = weight.side/2;
  neu->i0 = neu->j0 = ((neu->memside)/2)-n;
  neu->side = 2*n;
  for(i=0; i<2*n; i+=2)
    for(j=0; j<2*n; j+=2) {
      if((*(chel(*neu,i,j)) &&
	  *(chel(*neu,i+1,j+1))) ||
	 (*(chel(*neu,i+1,j)) &&
	  *(chel(*neu,i,j+1))))
	*(chel(*neu,i,j)) =
	  *(chel(*neu,i+1,j+1)) =
	  *(chel(*neu,i+1,j)) =
	  *(chel(*neu,i,j+1)) = 0;
      else {
	if(*chel(*neu,i,j)) {
	  *chel(*neu,i,j)     = 0;
	  *chel(*neu,i+1,j+1) = 1;
	} else
	  if(*chel(*neu,i+1,j+1)) {
	    *chel(*neu,i,j)     = 1;
	    *chel(*neu,i+1,j+1) = 0;
	  } 
	if(*chel(*neu,i+1,j)) {
	  *chel(*neu,i+1,j) = 0;
	  *chel(*neu,i,j+1) = 1;
	} else
	  if(*chel(*neu,i,j+1)) {
	    *chel(*neu,i+1,j) = 1;
	    *chel(*neu,i,j+1) = 0;
	  }
      }
    }
  for(i=0; i<2*n-1; i++) {
    if(i%2) {
      for(j=1; j<n; j++) {
	if(!*(chel(*neu,i-1,2*j-1)) 
	   && !*(chel(*neu,i  ,2*j-1)) 
	   && !*(chel(*neu,i-1,2*j  )) 
	   && !*(chel(*neu,i  ,2*j  ))) {
	  if(dr250()<(coup_spit(el(weight,i,  2*j-1))*
		      coup_spit(el(weight,i+1,2*j  )))/
	            ((coup_spit(el(weight,i,  2*j-1))*
		      coup_spit(el(weight,i+1,2*j  ))) +
		     (coup_spit(el(weight,i,  2*j  ))*
		      coup_spit(el(weight,i+1,2*j-1)))))
	    *(chel(*neu,i, 2*j-1)) = 
	    *(chel(*neu,i+1, 2*j)) = 1;
	  else
	    *(chel(*neu,i,   2*j)) = 
	    *(chel(*neu,i+1, 2*j-1)) = 1;
	}
      }
    } else {
      for(j=0; j<n; j++) {      
	if(((i==0) &&
	    !*(chel(*neu,0,2*j)) &&
	    !*(chel(*neu,0,2*j+1))) ||
	   ((i>0) &&
	    !*(chel(*neu,i-1,2*j)) &&
	    !*(chel(*neu,i-1,2*j+1)) &&
	    !*(chel(*neu,i,2*j)) &&
	    !*(chel(*neu,i,2*j+1)))) {
	  coup_mul(temp,*(el(weight,i,  2*j)),
		   *(el(weight,i+1,  2*j+1)));
	  coup_mul(temp+1,*(el(weight,i,  2*j+1)),
		   *(el(weight,i+1,  2*j)));
	  coup_add(temp+2,temp[0],temp[1]);
	  coup_div(temp+3,temp[0],temp[2]);
	  if(dr250()<coup_spit(temp+3))
	    *(chel(*neu,i,  2*j)) = 
	      *(chel(*neu,i+1,2*j+1)) = 1;
	  else
	    *(chel(*neu,i,  2*j+1)) = 
	      *(chel(*neu,i+1,2*j))   = 1;
	}
      }
    }
  }
  for(i=0; i<10; i++)
    coup_kill(temp+i);
}

void grow(matr *neu, matr weight, matr alt)
/* Pre: neu->memside, weight.memside, weight.side even */
/* neu is taken to be allocated and initialized */
/* Post: neu->side = weight.side = alt.side+2*/
     /* must work even if neu->a == weight.a */
{
  int n;
  int i,j;
  couple temp[10];

  for(i=0; i<10; i++)
    coup_init(temp+i);

  n = weight.side/2;
  neu->i0 = neu->j0 = ((neu->memside)/2)-n;
  neu->side = 2*n;

  for(i=0; i<2*n; i+=2)
    for(j=0; j<2*n; j+=2) {
      if(i>0 && j>0)
	coup_ass(temp,*(el(alt,i-1,j-1)));
      else
	coup_set(temp,0);
      if(i>0 && j<(2*n-2))
	coup_ass(temp+1,*(el(alt,i-1,j)));
      else
	coup_set(temp+1,0);
      if(i<(2*n-2) && j>0)
	coup_ass(temp+2,*(el(alt,i,j-1)));
      else
	coup_set(temp+2,0);
      if(i<(2*n-2) && j<(2*n-2))
	coup_ass(temp+3,*(el(alt,i,j)));
      else
	coup_set(temp+3,0);
      coup_set(temp+4,1);
      coup_sub(temp+4,temp[4],temp[0]);       
      coup_sub(temp+4,temp[4],temp[1]);
      coup_sub(temp+4,temp[4],temp[2]);       
      coup_sub(temp+4,temp[4],temp[3]);
      coup_mul(temp+5,*(el(weight,i,j)),*(el(weight,i+1,j+1)));
      coup_mul(temp+6,*(el(weight,i+1,j)),*(el(weight,i,j+1)));
      coup_add(temp+7,temp[5],temp[6]);
      coup_div(temp+8,temp[5],temp[7]);
      coup_div(temp+9,temp[6],temp[7]);
      coup_mul(temp+8,temp[8],temp[4]);
      coup_mul(temp+9,temp[9],temp[4]);
      coup_add(el(*neu,i,j),temp[3],temp[8]);
      coup_add(el(*neu,i+1,j+1),temp[0],temp[8]);
      coup_add(el(*neu,i,j+1),temp[2],temp[9]);
      coup_add(el(*neu,i+1,j),temp[1],temp[9]);
    }

  for(i=0; i<10; i++)
    coup_kill(temp+i);
}

void newrenew(matr *neu, couple *fact)
/* Pre: neu->memside, neu->side positive and even */
/* neu and fact are taken to be allocated and initialized */
{
  int i,j;
  couple temp[7];

  for(i=0; i<7; i++)
    coup_init(temp+i);
  coup_set(fact,1);
  for(i=0; i<neu->side; i+=2)
    for(j=0; j<neu->side; j+=2) {
      coup_mul(temp  ,*(el(*neu,i,j)),  *(el(*neu,i+1,j+1)));
      coup_mul(temp+1,*(el(*neu,i+1,j)),*(el(*neu,i,j+1)));
      coup_add(temp+2,temp[0],temp[1]);
      coup_mul(fact,*fact,temp[2]);
      coup_div(temp+3,*(el(*neu,i+1,j+1)),temp[2]);
      coup_div(temp+4,*(el(*neu,i+1,j)),temp[2]);
      coup_div(temp+5,*(el(*neu,i,j+1)),temp[2]);
      coup_div(temp+6,*(el(*neu,i,j)),temp[2]);
      coup_ass(el(*neu,i  ,j  ),temp[3]);
      coup_ass(el(*neu,i  ,j+1),temp[4]);
      coup_ass(el(*neu,i+1,j  ),temp[5]);
      coup_ass(el(*neu,i+1,j+1),temp[6]);
    }
  for(i=0; i<7; i++)
    coup_kill(temp+i);
}

matr goodedgprob(matr weight, couple *fact, chatr *match)
     /* pre: weight.side even and positive */
     /*      weight can be modified (and will!)*/
{
  matr a1, a2;
  int n;
  int i;
  couple temp[5];
  int howfar = 0;

  for(i=0; i<5; i++)
    coup_init(temp+i);

  n = weight.side/2;
  coup_set(fact,1);
  fprintf(stderr, "Doing step 1...                                   end mark V\n");
  for(i=0; i<n; i++) {    
    while(i * 60 / (n-1) > howfar) {
      fprintf(stderr, ".");
      howfar ++;
    }
    newrenew(&weight,&(temp[0]));
    weight.side -= 2;
    weight.i0++; weight.j0++;
    /* matr_out(stdout,10,weight,10); */
    coup_mul(fact, *fact, temp[0]);
  }
  if(mode&6) {
    matr_create(&a1,2*n);
    chatr_create(match,2*n);
    a1.i0=a1.j0=weight.i0=weight.j0=match->i0=match->j0=n-1;
    a1.side = weight.side = match->side = 2;
    a2=a1;
    newrenew(&weight,&(temp[0]));
    coup_mul(temp,   (*(el(weight,0,0))), (*(el(weight,1,1)))); 
    coup_mul(temp+1, (*(el(weight,0,1))), (*(el(weight,1,0)))); 
    coup_add(temp+2, temp[0], temp[1]);
    coup_div(temp+3, temp[0], temp[2]);
    coup_div(temp+4, temp[1], temp[2]);
    coup_ass(el(a1,0,0),temp[3]);
    coup_ass(el(a1,1,1),temp[3]);
    coup_ass(el(a1,0,1),temp[4]);
    coup_ass(el(a1,1,0),temp[4]);
    if(dr250()<coup_spit(temp+3)) {
      *(chel(*match,0,0))=1;
      *(chel(*match,1,1))=1;
    } else {
      *(chel(*match,0,1))=1;
      *(chel(*match,1,0))=1;
    }
    fprintf(stderr, "\nDoing step 2...                                   end mark V\n");
    howfar = 0;
    for(i=n-2; i>=0; i--) {
      while((n-2-i) * 60 / (n-2) > howfar) {
	fprintf(stderr, ".");
	howfar ++;
      }
      weight.side+=2; weight.i0--; weight.j0--;
      newrenew(&weight,&(temp[0]));
      if(mode&4)
	grrand(match,weight);
      if(mode&2) {
	grow(&a2,weight,a1);
	a1 = a2; 
      }
    }
  }
  
  for(i=0; i<5; i++)
    coup_kill(temp+i);
  fprintf(stderr, "\ndone.\n");
  return a1;
}

main(int argc, char *argv[])
{
  int n;
  matr mat,m2;
  chatr match;
  couple fact;

  if(argc<2)
    mode=4;
  else 
    mode=atoi(argv[1]);
  r250_init(time(NULL));
  fscanf(stdin,"%d",&n);
  matr_create(&mat,2*n);
  matr_in(stdin,10,&mat);
  coup_init(&fact);
  m2 = goodedgprob(mat,&fact,&match);
  if(mode&1) {
    coup_out(stdout,10,fact);
    fprintf(stdout,"\n");
  }
  if(mode&2)
    matr_out(stdout,10,m2,10);
  if(mode&4)
    chatr_out(stdout,match);
}