/* In the standard form, reads from a file a rack and a cocycle with
 * integers values, then a monomial and finally a list of derivatives
 * to apply to it. Then, it applies them.
 * With the parameter -s reads any line after the cocycle as a monomial
 * and searches for each of them for a derivation in the same degree which
 * gives nonzero on it.
 * With the parameter -sb reads one monomial and one derivation. The derivation
 * can have a degree < degree of monomial. It applies the derivation to the
 * monomial and then searches for a way to complete it for it to give a nonzero
 * result.
 * With -ss acts like -s except that it reads several lines in the form
 * n mon, where n is an integer, and looks for a nonzero derivation of the element
 * which is the linear combination of monomials with coefficients the 'n' in each line.
 * With the parameter -d reads one or more monomials and one or more derivations,
 * divided by a line beginning with "%D". Then takes derivatives for each monomial
 * and each derivation.
 * An example of the syntax of the input file (for standard form) is:
 * ----------------
%Rack
abcdefghij
%Braid:
aefgbcdhij
ebhiafgcdj
fhcjeagbid
gijdefahbc
bacdehifgj
cbadhfjeig
dbcaijghef
acbdfeghji
adcbgfejih
abdcegfihj
%Cocycle
-1 -1 -1 -1 -1 -1 -1 -1 -1 -1
-1 -1 -1 -1 -1 -1 -1 -1 -1 -1
-1 -1 -1 -1 -1 -1 -1 -1 -1 -1
-1 -1 -1 -1 -1 -1 -1 -1 -1 -1
-1 -1 -1 -1 -1 -1 -1 -1 -1 -1
-1 -1 -1 -1 -1 -1 -1 -1 -1 -1
-1 -1 -1 -1 -1 -1 -1 -1 -1 -1
-1 -1 -1 -1 -1 -1 -1 -1 -1 -1
-1 -1 -1 -1 -1 -1 -1 -1 -1 -1
-1 -1 -1 -1 -1 -1 -1 -1 -1 -1
%Monomial
abacabacdabacabacdabadcabacdefegefeghihj
%Derivatives
jihifjhgifhedijeidfeigdjihdihgidchfhcbea
  * ----------------------------
  * This computes derivatives of the longest monomial in the
  * Nichols algebra generated by transpos. in S_5. */


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

#define MAXRACK 35 	/* Max size of rack allowed */
#define MAXDEG 50	/* Max degree allowed */
#define MAXTERM 4000000 /* Max number of terms in each derivation */
#define MAXMON 10000 /* Max number of monomials to derivate with option -sb */
#define MAXDER 20  /* Max number of derivations with option -d */

static char rack[MAXRACK];		/* the elements of the rack */
static int cardx;			/* cardinality of the rack */
static int braid[MAXRACK*MAXRACK];	/* braiding */
static int cocycle[MAXRACK*MAXRACK];/* 2-cocycle */
static int mon[MAXDEG];			/* monomial to derivate */
static int der[MAXDEG];			/* derivations to apply */
static int dertoim[MAXDEG];			/* derivations to apply for printing matters */
static int *smon[MAXMON];			/* monomials to derivate with option -sb */
static int *sder[MAXDER];			/* derivations to apply with option -d */
static int degmon;					/* degree of mon[] */
static int degder;					/* number of derivations */
static int *toder[MAXTERM];			/* terms to derivate */
static int *result[MAXTERM];		/* results of sucessive derivatives */
static int numterms;				/* number of terms (<=MAXTERM) */
static int numtermsres;				/* number of terms resulting */
static int srch=0;					/* do we want it to search ? */
static int count[MAXRACK];			/* counts number of times each element in the rack appears */
static float percentcompleted=0;	/* computes how much is completed with option -s */
static int whentoprint=6;			/* the higher, the more info printed with option -s */
static int nummonomials=0;			/* number of monomials to look for with option -sb */
static int numderivs=0;				/* number of derivations to compute with option -d */

int compterms(void const*, void const*);
int cmpcount(void const*, void const*);
int charinstr(char*, char);

readrackcocmonder (FILE *in) {  // reads rack, cocycle, monomial and derivations to apply
	int siga,i,j,n;
	char c;
	siga=1; i=0;
	while (siga) { // read rack set
		c=getc(in);
		if (c=='%') {
			while (getc(in) != '\n');
			continue;
		}
		if (c=='\n')
			siga=0;
		else
			rack[i++]=c;
	}
	rack[i]='\0';
	cardx=i;
	i=0;
	while (i<cardx*cardx) { // read rack braid
		c=getc(in);
		if (c=='%') {
			while (getc(in) != '\n');
			continue;
		}
		else
			if (charinstr(rack, c)) {
				braid[i++]=charinstr(rack, c)-1;
			}
	}
	getc(in); while (getc(in) != '\n'); // skip a line with a comment
	for (i=0; i<cardx*cardx; i++) { // read cocycle
		fscanf(in, "%d", &n);
		cocycle[i]=n;
	}
	getc(in);
	siga=1; i=0;
	while (siga) { // read monomial
		c=getc(in);
		if (c==EOF)
			return(1);
		if (c=='%') {
			if (getc(in)=='D' && srch==3) {
				while (getc(in) != '\n');
				siga=0;
			}
			else {
				while (getc(in) != '\n');
				continue;
			}
		}
		if (c=='\n') {
			if (srch==1 || srch==3 || srch==4) {
				smon[nummonomials]=calloc(MAXDEG, sizeof(int));
				for (j=0; j<i; j++)
					smon[nummonomials][j+1]=mon[j];
				smon[nummonomials++][0]=i;
				i=0;
			}
			else
				siga=0;
		}
		else if ( (srch != 4) && (charinstr(rack,c))) {
			if (i==MAXDEG-1) {
				printf("I'm sorry, degree of monomial > max allowed degree.\n");
				exit(1);
			}
			else
				mon[i++]=charinstr(rack,c)-1;
		}
		else if ( (srch==4) ) {
			if (i==MAXDEG-1) {
				printf("I'm sorry, degree of monomial > max allowed degree.\n");
				exit(1);
			}
			else if (i==0) {
				ungetc(c, in);
				mon[i++]=atoif(in);
//	printf("\nread mon[%d]=%d ", nummonomials, mon[0]);
			}
			else {
				mon[i++]=charinstr(rack,c)-1;
//	printf("%c", rack[mon[i-1]]); fflush(stdout);
			}
		}
	}
	degmon=i;
	siga=1; i=0;
	if (srch==1 || srch==4)
		return(1);
	while (siga) { // read derivations
		c=getc(in);
		if (c=='%') {
			while (getc(in) != '\n');
			continue;
		}
		if (c=='\n') {
			if (srch==3) {
				sder[numderivs]=calloc(MAXDEG, sizeof(int));
				for (j=0; j<i; j++)
					sder[numderivs][j+1]=der[j];
				sder[numderivs++][0]=i;
				i=0;
			}
			else
				siga=0;
		}
		else if (charinstr(rack,c)) {
			if (i==degmon && srch==0) {
				printf("Degree of derivations > degree of monomial = %d. Result.... 0\n", degmon);
				exit(0);
			}
			else
				der[i++]=charinstr(rack, c)-1;
		}
		if (c==EOF)
			siga=0;
	}
	degder=i;
}


charinstr(char *s, char c) {
	int i=0;
	while (*(s+i) && *(s+i)!=c)
		i++;
	if (*(s+i))
		return(i+1);
	else
		return(0);
}


beginstr(char *a, char *b) {
	while (*a && *b && *a==*b) {
		a++; b++;
	}
	if (*a)
		return(0);
	else
		return(1);
}

equalstr(char *a, char *b) {
	while(*a && *b && *a==*b) {
		a++; b++;
	}
	if (*a || *b)
		return(0);
	else
		return(1);
}



deriveterm(int i, int d, int l) {
	int j, k, doi;
	for (j=1; j<=l; j++) {
		if (toder[i][j]==d) {
			result[numtermsres]=calloc(l+1, sizeof(int));
			for (k=1; k<j; k++)
				result[numtermsres][k]=toder[i][k];
			result[numtermsres][0]=toder[i][0];
			for (k=j+1; k<=l; k++) {
				doi=toder[i][k];
				result[numtermsres][k-1]=braid[d*cardx+doi];
				result[numtermsres][0] *= cocycle[d*cardx+doi];
			}
			numtermsres++;
		}
	}
}


compterms(void const *a, void const *b) {
	int i=1;
	int **aa=(int **)a;
	int **bb=(int **)b;
	while (i<=degmon) {
		if ((*aa)[i]-(*bb)[i])
			return((*aa)[i]-(*bb)[i]);
		i++;
	}
	return(0);
}

collectterms(int ldegmon) {
	int i, j, k, s;
	numterms=0;
	degmon=ldegmon;
	qsort(result, numtermsres, sizeof(int*), compterms);
	for (i=0; i<numtermsres;) {
		s=result[i][0];
		for (j=i+1; j<numtermsres && !compterms(result+i, result+j); s+=result[j++][0]);
		if (s) {
			toder[numterms]=calloc(degmon+1, sizeof(int));
			toder[numterms][0]=s;
			for (k=1; k<=degmon; k++)
				toder[numterms][k]=result[i][k];
			numterms++;
		}
		i=j;
	}
	for (i=0; i<numtermsres; i++)
		free(result[i]);
}



derive(int printlevel) { // needed variables: degder, numterms, degmon, der, toder
	int i,j,cpdegder=degder;
	while (degder) {
		if (printlevel==1) {
			printf("(%c, %4d term%s)", rack[*der], numterms, (numterms==1)?" ":"s"); fflush(stdout);
			if (!((cpdegder-degder+1)%5))
				printf("\n");
		}
		numtermsres=0;
		for (i=0; i<numterms; i++) {
			deriveterm(i, *der, degmon);
			free(toder[i]);
		}
		for (i=0; i<degder; i++)
			der[i]=der[i+1];
		collectterms(degmon);
		degder--;
		degmon--;
	}
	qsort(toder, numterms, sizeof(int*), compterms);
	if (!numterms) {
		if (printlevel==1)
			printf("\n\nThis gave 0. Sorry...\n\n");
		else
			printf("0 \n");
	}
	else {
		if (printlevel==1)
			printf("\n");
		for (i=0; i<numterms; i++) {
			if (printlevel==1)
				printf("\n");
			printf("%+d * ", toder[i][0]);
			for (j=1; j<=degmon; j++) {
				printf("%c", rack[toder[i][j]]);
			}
			printf(" ");
		}
		if (printlevel==1)
			printf("\n\n");
		else
			printf("\n");
	}
}

cmpcount(void const *a, void const *b) {
	int *aa=(int *)a;
	int *bb=(int *)b;
//	return(count[*aa]-count[*bb]); // orders according to count
//	return(count[*bb]%2-count[*aa]%2); // puts first odds
	if (count[*aa]&1) {
		if (count[*bb]&1)
			return(count[*aa]-count[*bb]);
		else
			return (-1);
	}
	else {
		if (count[*bb]&1)
			return(1);
		else
			return(count[*aa]-count[*bb]);
	}
}

isitnonzero(int locdegder, int **loctoder, int locnumterms, int locdegmon) {
	int i,j,k;
	int order[cardx];
	int loccount[cardx];
	int *cptoder[locnumterms];
	if (!locdegmon && locnumterms) {
		return(1);
	}
//	printf("Will run isitnonzero with %d numterms, %d degmon\n", numterms, degmon);
	for (i=0; i<locnumterms; i++) {
		cptoder[i]=calloc(degmon+1, sizeof(int));
		for (j=0; j<=degmon; j++)
			cptoder[i][j]=loctoder[i][j];
	}
	for (i=0; i<cardx; i++)
		order[i]=i;
	countelems(locnumterms, loctoder, locdegmon);
	qsort(order, cardx, sizeof(int), cmpcount);
	for (i=0; i<cardx; i++)
		loccount[i]=count[i];
//	printf("Order: "); for(i=0; i<cardx; i++) printf("(%c,%d) ",rack[order[i]], count[order[i]]); printf("\n");
	for (i=0; i<cardx; i++) {
		if (locdegder<10) {
			percentcompleted += 1./power(cardx, locdegder+1);
		}
		if (!loccount[order[i]])
			continue;
		for (j=0; j<locnumterms; j++) {
			toder[j]=calloc(locdegmon+1, sizeof(int));
			for (k=0; k<=locdegmon; k++)
				toder[j][k]=cptoder[j][k];
		}
//	printf("Will take der %c with der order %d... ", rack[order[i]], locdegder);
		der[locdegder]=order[i];
		dertoim[locdegder]=order[i];
//		dertoim[locdegder]=i; // this is in case one wants to have a better idea of the remaining time
		numtermsres=0;
		for (j=0; j<locnumterms; j++) {
			deriveterm(j, order[i], locdegmon);
			free(toder[j]);
		}
//	printf(" got %d terms... ", numtermsres);
		collectterms(locdegmon);
//	printf(" and %d terms after collect\n", numterms);
		if (numterms) {
			int *passtoder[numterms];
			int loclocnumterms=numterms;
//	printf("run isitnonzero, degmon=%d, ", locdegmon-1);
//	for (k=0; k<=locdegder; k++) {printf("%c", rack[der[k]]);} printf("\n");
//	printf("Terms: "); for(k=0; k<numterms; k++) {printf("  %d ", toder[k][0]); for(j=1; j<=locdegmon-1; j++) printf("%c", rack[toder[k][j]]);} printf("\n");
			for (j=0; j<numterms; j++) {
				passtoder[j]=calloc(locdegmon, sizeof(int));
				for (k=0; k<locdegmon; k++) {
					passtoder[j][k]=toder[j][k];
				}
				free(toder[j]);
			}
			if (isitnonzero(locdegder+1, passtoder, numterms, locdegmon-1)) {
				for (j=0; j<loclocnumterms; j++)
					free(passtoder[j]);
				return(1);
			}
			else {
				for (j=0; j<loclocnumterms; j++)
					free(passtoder[j]);
//				return(0);
			}
		}
	}
	for (i=0; i<locnumterms; i++)
		free(cptoder[i]);
	if (locdegder <= whentoprint) {
		if (locdegder==whentoprint) {
			printf("\ncompleted: %f. Der: ", percentcompleted*(cardx-1)/cardx);
			for (i=0; i<locdegder; i++)
				printf("%c", rack[dertoim[i]]);
			printf(".  ");
		}
		printf(" %d ", locdegder); fflush(stdout);
	}
	return(0);
}


power(int a, int b) {
	int r=1;
	for (; b>0; b--)
		r *= a;
	return(r);
}


atoif(FILE *in) {
	int r=0, s=1;
	char c;
	while ( (((c=getc(in))<'0') || (c>'9')) && (c != '-') && (c != '+') );
	if (c=='-')
		s=-1;
	if ((c >= '0') && (c <= '9'))
		r=c-'0';
	while ( (c=getc(in)) != ' ' )
		r = r*10+c-'0';
	return(r*s);
}


countelems(int locnumterms, int **loctoder, int locdegmon) {
	int i,j;
	for (i=0; i<cardx; i++)
		count[i]=0;
	for (i=0; i<locnumterms; i++) {
		for (j=1; j<=locdegmon; j++) {
			count[loctoder[i][j]]++;
		}
	}
}



main (int argc, char **argv) {
	FILE *in;
	int i,j,k;
	if (beginstr("-s", argv[1])) {
		srch=1;
		if (equalstr("-sb", argv[1]))
			srch=2;
		if (equalstr("-ss", argv[1]))
			srch=4;
		argv++;
	}
	if (beginstr("-d", argv[1])) {
		srch=3;
		argv++;
	}
	in = fopen(argv[1],"r");
	if (!in) {
		printf("Couldn't open file %s...\n", argv[1]);
		exit(1);
	}
	readrackcocmonder(in);
//	printf("Read rack with %d elements. Monomial of deg %d and der of deg %d\n", cardx, degmon, degder);
	fclose(in);
	if (srch) {
//		whentoprint=(whentoprint>degmon/2)?whentoprint:degmon/2;
		if (srch==1) {
			for (j=0; j<nummonomials; j++) {
				int *starttoder[0];
				percentcompleted=0;
				degmon=smon[j][0];
				starttoder[0]=calloc(degmon+1, sizeof(int));
				starttoder[0][0]=1;
				printf("Let's deal with: ");
				for (i=1; i<=degmon; i++) {
					starttoder[0][i]=mon[i-1]=smon[j][i];
					printf("%c", rack[smon[j][i]]);
				}
				printf("\n");
				if (isitnonzero(0, starttoder, 1, degmon)) {
					printf("\nThe monomial   ");
					for (i=1; i<=smon[j][0]; i++)
						printf("%c", rack[smon[j][i]]);
					printf("    gave nonzero with\n");
					for (i=0; i<smon[j][0]; i++)
						printf("%c", rack[der[i]]);
					printf("\n\n");
				}
				else {
					printf("\n\nSorry.... the monomial   ");
					for (i=1; i<=smon[j][0]; i++)
						printf("%c", rack[smon[j][i]]);
					printf("    is zero.\n\n");
				}
			}
		}
		else if (srch==2) {
			int cpmon[degmon+1], cpder[degmon];
			printf("Let's deal with: ");
			for (i=1; i<=degmon; i++) {
				printf("%c", rack[mon[i-1]]);
			}
			printf("\n\n");
			toder[0]=calloc(degmon+1, sizeof(int));   // put mon into der[0] to begin process
			toder[0][0]=1;
			for (i=1; i<=degmon; i++) {
				toder[0][i]=cpmon[i]=mon[i-1];
				cpder[i]=der[i-1];
			}
			numterms=1;
			cpmon[0]=degmon;
			cpder[0]=degder;
			derive(1);
			percentcompleted=0;
			if (isitnonzero(0, toder, numterms, degmon)) {
				printf("\nThe monomial   ");
				for (i=1; i<=cpmon[0]; i++)
					printf("%c", rack[cpmon[i]]);
				printf("    gave nonzero with\n");
				for (i=1; i<=cpder[0]; i++)
					printf("%c", rack[cpder[i]]);
				for (i=0; i<cpmon[0]-cpder[0]; i++)
					printf("%c", rack[der[i]]);
				printf("\n\n");
			}
			else {
				printf("\n\nSorry.... it gives zero with those derivations \n\n  ");
			}
		}
		else if (srch==3) {
			for (j=0; j<nummonomials; j++) {
				printf("\nMonomial: ");
				for (i=1; i<=smon[j][0]; i++)
					printf("%c", rack[smon[j][i]]);
				printf(":\n");
				for (k=0; k<numderivs; k++) {
					printf("\tDerivative ");
					for (i=1; i<=sder[k][0]; i++)
						printf("%c", rack[sder[k][i]]);
					printf(":\n\t\t");
					toder[0]=calloc(smon[j][0]+1, sizeof(int));   // put mon into der[0] to begin process
					toder[0][0]=1;
					for (i=1; i<=smon[j][0]; i++)
						toder[0][i]=smon[j][i];
					for (i=1; i<=sder[k][0]; i++) {
						der[i-1]=sder[k][i];
					}
					numterms=1;
					degmon=smon[j][0];
					degder=sder[k][0];
					derive(3);
				}	
			}
		}	
		else if (srch==4) {
			int cpdegmon;
//	printf("entering in search mode = 4...\n"); degmon=smon[0][0]-1;
			degmon=smon[0][0]-1;
			cpdegmon=degmon;  // save it to print the derivatives (isitnonzero changes its value)
			printf("I'll consider the combination\n");
			for (i=0; i<nummonomials; i++) {
				printf("%d ",smon[i][1]);
				for(j=1; j<=degmon; j++)
					printf("%c",rack[smon[i][j+1]]);
				printf("\n");
			}
			for (i=0; i<nummonomials; i++) { // put smon into toder to begin process
				if (smon[i][0]-1 != degmon) {
					printf("\nToo bad. Monomial # %d has different degree from first monomial!\n\n", i);
					abort();
				}
				toder[i]=calloc(degmon+1, sizeof(int));
				for (j=0; j<=degmon; j++)
					toder[i][j]=smon[i][j+1];
			}
			if (isitnonzero(0, toder, nummonomials, degmon)) {
				printf("\nThe element gave nonzero with\n");
				for (j=0; j<cpdegmon; j++)
					printf("%c", rack[der[j]]);
				printf("\n\n");
			}
			else {
				printf("\n\n...well, I guess this is good: it gave 0.\n\n");
			}
		}
	}
	else {
		toder[0]=calloc(degmon+1, sizeof(int));   // put mon into der[0] to begin process
		toder[0][0]=1;
		for (i=1; i<=degmon; i++)
			toder[0][i]=mon[i-1];
		numterms=1;
		derive(1);
		if (numterms) {
			printf("Number of terms: %d\n", numterms);
			countelems(numterms, toder, degmon);
			printf("Number of elems in these terms by element:\n");
			for (i=0; i<cardx; i++)
				printf("(%c : %d) ", rack[i], count[i]);
			printf("\n\n");
		}
	}
}