/************Disclaimer*************************************/
/*This work was supported by the University of Minnesota   */
/*Undergraduate Research Opportunities Program		   */
/*Under the direction of Professor Victor Reiner.	   */
/*Programmed by Michael Barany.  No liability, yada, yada, */
/*yada... Thanks for your support!			   */
/***********************************************************/


/**********Known bugs/issues*************/
/* Tuttes with coefficients >2M	        */
/* 				        */
/*					*/
/****************************************/

#include <stdio.h>
#include <malloc.h>
#include <stdlib.h>


/* basic matrix and row-reduce prototypes */
int **imatrix(int d,int n,int start);
void getcols(int **m, FILE *usermatrix, int d, int n,int start);
void nrerror(char errortext[]);
void freeimatrix(int **m,int d, int n,int start);
void printimatrix(int **m, int d, int n,int start);
void copyimatrix(int **m, int **newmat, int d, int n,int start);
void swaprows(int **m, int row1, int row2, int n);
void swapcols(int **m, int col1, int col2, int d);
void setuppivot(int **m, int pivcol, int pivrow, int d, int n);
void elim(int **m, int pivcol, int pivrow, int d, int n);
int iabs(int a);
void mod(int **m, int d, int n, int p);
void rowreduce(int **m, int d, int n, int p);

/*polynomial prototypes                                 */
/*A_(i,j)(x^i)(y^j)                                     */
void initpoly(int **m,int d,int n);
void printpoly(int **m,int d,int n); 

int *ivector(int low, int high);
void printvector(int *v,int n);

/*Tutte prototypes*/
int justloopsnbridges(int **m, int d, int n);
int done(int *record,int n);
int lastreczero(int *record,int n);
int lastrecmod(int *record, int n);
void sendtotutte(int **m, int **tutte,int d,int n);
void deletein(int **m, int d, int *n, int edge);
void contractin(int **m, int *d, int *n, int edge);
void d1totutte(int **m, int **tutte, int n);
int rebuild(int **m1,int **m, int d,int n,int *d1,int *n1,int *record,\
int p);

int  main(int argc, char *argv[])
{
	FILE *usermatrix;
	int **m,**m1,**tutte;
	int d,d1,n,n1,th,i,j,p,notfinished,nextedge;
	int *record;

	if (argc!=5){
		printf("Usage: ./tutte.out d n th matrix.dat\n");
	}
	else if (!(usermatrix=fopen(argv[4], "r"))){
		printf("Error opening file %s\n", argv[4]);
	}
	else{
		/*get the user's matrix*/
		d=atoi(argv[1]);
		n=atoi(argv[2]);
		p=atoi(argv[3]);
		d1=d;
		n1=n;
		m1=imatrix(d1,n1,1);
		getcols(m1,usermatrix,d1,n1,1);
		printimatrix(m1,d1,n1,1);
		/*store a reference copy*/
		m=imatrix(d,n,1);
		copyimatrix(m1,m,d1,n1,1);

		/*initialize tutte*/
		tutte=imatrix(d,n,0);
		initpoly(tutte,d,n);

		/*create del/mod recording vector*/
		record=ivector(1,n);

		notfinished=1; /*indicates not finished*/
		while(notfinished){ /*while not finished*/
			rowreduce(m1,d1,n1,p);
			nextedge=justloopsnbridges(m1,d1,n1);
			while(nextedge>0){ /*while not at base case*/
				i=lastreczero(record,n);
				*(record+i)=1;
				contractin(m1,&d1,&n1,nextedge);
                		rowreduce(m1,d1,n1,p);   
				nextedge=justloopsnbridges(m1,d1,n1);
        		}
			if(nextedge==-1) d1totutte(m1,tutte,n1);
			else sendtotutte(m1,tutte,d1,n1);
			notfinished=rebuild(m1,m,d,n,&d1,&n1,record,p);
		}
		printpoly(tutte,d,n);
		return 0;
	}
}

int **imatrix(int d,int n, int start)
{
	int **m;
	int i;
	/* allocate pointers to rows */
	m=(int **) malloc((unsigned) (d-start+1)*sizeof(int*));
	if (!m) nrerror("allocation failure 1 in imatrix()");
	m -= start;
	/* allocate rows and set pointers to them */
	for(i=start;i<=d;i++)
		{
		m[i]=(int *)malloc((unsigned) (n-start+1)*sizeof(int));
		if (!m[i]) nrerror("allocation failure 2 in imatrix1()");
		m[i] -= start;
		}
	return m;
}

int *ivector(int low, int high){
	int *x;
	x=(int *) malloc((unsigned) (high-low+1)*sizeof(int));
	if (!x) nrerror("allocation failure in ivector()");
	return (x-low);
}

void initpoly(int **m,int d,int n){
	int i,j;
	for (i=0;i<=d;i++){
		for (j=0;j<=n;j++){
			*(*(m+i)+j)=0;
		}
	}
}

void printimatrix(int **m, int d, int n, int start)
{
        int i, j;
        for(i=start;i<=d;i++)
        {
                printf("\n[ ");
                for(j=start;j<=n;j++)
                {
                        printf("%d  ", m[i][j]);
                }
                printf("]");
        }
        printf("\n");
}

void printpoly(int **m, int d, int n)
{
        int i, j, count=0;
        for(i=0;i<=d;i++) {
                for(j=0;j<=n;j++) {
			if(m[i][j]){
		  if (count==0) {
                        printf("%d",m[i][j]);
                        if(i) printf(" x^%d",i);
                        if(j) printf(" y^%d",j);
		  }
		  else {
			printf(" + %d",m[i][j]);
			if(i) printf(" x^%d",i);
			if(j) printf(" y^%d",j);
		  }
		  count++;
               		}
		 }
        }
	printf("\n");
}

void nrerror(char errortext[])
{
	void exit();

	fprintf(stderr, "Numerical Recipes run-time error...\n");
	fprintf(stderr, "%s\n", errortext);
	fprintf(stderr, "...now exiting to system...\n");
	exit(1);
}

void getcols(int **m, FILE *usermatrix, int d, int n, int start)
{
	int i, j;
	int temp;
	for (j=start; j<=n; j++)
	{
		for (i=start; i<=d; i++)
		{
			fscanf(usermatrix,"%d", &temp);
			*(*(m+i)+j)=temp;
		}
	}
}

void freeimatrix(int **m, int d, int n,int start)
{
	int i;
	for(i=d; i>=start; i--) free((char*) (m[i]+start));
	free((char*) (m+start));
}


void copyimatrix(int **m, int **newmat, int d, int n, int start) {
	int i, j;

        for (i=start; i<=d; i++)
        {
                for (j=start; j<=n; j++)
                {
                        *(*(newmat+i)+j)=*(*(m+i)+j);
                }
        }	
}

void swaprows(int **m, int row1, int row2, int n) {
	int j;
	int temp;

	for (j=1; j<=n;j++)
	{
		temp=*(*(m+row1)+j);
		*(*(m+row1)+j)=*(*(m+row2)+j);
		*(*(m+row2)+j)=temp;
	}
}

void swapcols(int **m, int col1, int col2, int d) {
	int i;
        int temp;
         
        for (i=1; i<=d;i++)
        {        
                temp=*(*(m+i)+col1);    
                *(*(m+i)+col1)=*(*(m+i)+col2);
                *(*(m+i)+col2)=temp;
        }
} 

void setuppivot(int **m, int pivcol, int pivrow, int d, int n) {
	int temp,i;
	if (pivrow<d){
		temp=pivrow+1;
		/*put a non-zero entry in the pivot position*/
		while ((*(*(m+pivrow)+pivcol))==0 && temp<=d) {
			swaprows(m,pivrow,temp, n);
			temp=temp+1;
		}
	}
}

void elim(int **m, int pivcol, int pivrow, int d, int n) {
	int i,j;
	int temp1,temp2;
	for(i=pivrow+1;i<=d;i++) {
		temp1=(*(*(m+i)+pivcol));
		temp2=(*(*(m+pivrow)+pivcol));
		for(j=1; j<=n; j++){
	*(*(m+i)+j)=temp2*(*(*(m+i)+j))-temp1*(*(*(m+pivrow)+j));
		}
	}
        for(i=1;i<=pivrow-1;i++) {  
                temp1=(*(*(m+i)+pivcol));
		temp2=(*(*(m+pivrow)+pivcol));
                for(j=1; j<=n; j++){        
        *(*(m+i)+j)=temp2*(*(*(m+i)+j))-temp1*(*(*(m+pivrow)+j));
                }
        }
}

void rowreduce(int **m, int d, int n,int p) {
	int i=1,j=1;
	while (i<=d && j<=n){
		setuppivot(m,j,i,d,n);
		if(*(*(m+i)+j)){
			elim(m,j,i,d,n);
			mod(m,d,n,p);
			i++;
			j++;
		}
		else {
			j++;
		}
	}
}

int iabs(int a){
	if (a>=0) return a;
	else return -1*a;
}

int justloopsnbridges(int **m, int d, int n){
	/*returns 0 if it's just loops and bridges, otherwise returns 
		first other edge*/
	int i=1,j=1,count=0,indicator=0;
	if(d==1 && n>=1) return -1;
	else{
		while(i<=d && indicator==0){
			count=0;
			j=1;
			while(j<=n && count<=1){
				if (*(*(m+i)+j)) count++;
				j++;
			}
			if (count>1) indicator++;
			else i++;
		}
		if (indicator) {
			j=1;
			count=0;
			while (count==0){
				if(*(*(m+i)+j)) count++;
				else j++;
			}
			return j;
		}
		else return 0;
	}
}

int done(int *record, int n) {
	/*returns 0 if record is entirely nulls and deletions*/
	int i=1, count=0;
	while(i<=n && count==0){
		count += *(record+i) % 2;
		i++;
	}
	return count;
}

int lastreczero(int *record,int n){
	int i=0, count=0;
	while(count==0){
		i++;
		if(i==n+1) count=1;
		else count += *(record+i);
	}
	i--;
	return i;
}

int lastrecmod(int *record, int n){
	int i=0, count=0;
	while(i<=n && count==0){
		i++;
		count += *(record+i) % 2;
	}
	return i;
}

void sendtotutte(int **m, int **tutte,int d,int n) {
	int i=1,j=1, dummy;
	int a=0,b=0;
	/* takes a forest with loops and converts to a tutte term*/

	if(d==0) b+=n;
	else if(n==0) a++;
	else {
		while(j<=n) {
			if(*(*(m+i)+j)==0) i++;
			else {
				j++;
				a++;
				i=1;
			}
			if(i==d+1) {
				j++;
				b++;
				i=1;
			}
		}
	}
	*(*(tutte+a)+b)+=1;
}

void d1totutte(int **m, int **tutte, int n){
	int j,count=0;
	for(j=1;j<=n;j++){
		if(*(*(m+1)+j)) count++;
	}
	if(count==0) *(*(tutte)+n)+=1;
	else{
		*(*(tutte+1)+n-count)+=1;
		for(j=1;j<=count-1;j++) *(*(tutte)+n-j)+=1;
	}
}

void printvector(int *v, int n) {
	int i;
	printf("record[ ");
	for(i=1;i<=n;i++) printf("%d ",*(v+i));
	printf("]\n");
}

void deletein(int **m, int d, int *n, int edge){
        int i,j;  
        for (i=1;i<=d;i++){
                for (j=edge+1;j<=*n;j++) *(*(m+i)+j-1)=*(*(m+i)+j);
        }
	(*n)--;
}


void contractin(int **m, int *d, int *n, int edge){
        int i,j,row,notzero;
        /*find the row corresponding to edge*/
        i=1;
        j=edge;
        notzero=0;
        while(notzero==0){
                notzero=*(*(m+i)+j);
                i++;
        } /*this will give 1 larger than the desired row*/
        row=i-1;
        for(i=1;i<row;i++){
                for(j=edge+1;j<=*n;j++) *(*(m+i)+j-1)=*(*(m+i)+j);
        }
        for(i=row+1;i<=*d;i++){
                for(j=edge+1;j<=*n;j++) *(*(m+i-1)+j-1)=*(*(m+i)+j);
        }
	(*d)--;
	(*n)--;
}

int rebuild(int **m1,int **m, int d,int n,int *d1,int *n1,int *record,\
int p) {
	int i,j,nextedge;
	*d1=d;
	*n1=n;
	copyimatrix(m,m1,d,n,1);
	j=lastrecmod(record,n);
	for(i=n;i>j;i--){
		rowreduce(m1,*d1,*n1,p);
		nextedge=justloopsnbridges(m1,*d1,*n1);
		if(*(record+i)==2) deletein(m1,*d1,n1,nextedge);
		if(*(record+i)==1) contractin(m1,d1,n1,nextedge);
	}
	rowreduce(m1,*d1,*n1,p);
	nextedge=justloopsnbridges(m1,*d1,*n1);
	deletein(m1,*d1,n1,nextedge);
	nextedge=done(record,n); /*actually, just reusing the var*/
	for(i=1;i<j;i++) *(record+i)=0;
	*(record+j)=2;
	return nextedge;
}

void mod(int **m, int d, int n, int p) {
        int i,j,temp;
        for(i=1;i<=d;i++){
                for(j=1;j<=n;j++){
                        temp=*(*(m+i)+j);
                        if(temp<0 && temp>-200000){
                                *(*(m+i)+j)-=p*temp;
                        }
                        else if(temp<-200000){
                                printf("you've hit a bug the programmer\
 didn't anticipate.  Please e-mail the graph you tried to run to \
rosevilledebate@yahoo.com\n");
                                exit(1);
                        }
                        *(*(m+i)+j)=temp % p;
                }
        }
}