/* Given a Grobner basis it returns the "PBW" basis
 * With -d option, it prints monomials with degree < argument (default=MAXDEG).
 * With -b option, it prints monomials with lexicog. order < argument (default=all monomials).
 * */



/* Example of input:
%Rack (for my program base)
abcde
%relations up to degree 7
a*a, b*b, c*c, d*d, e*e,
	d*c + c*a + b*d + a*b,
	e*a + c*b + b*e + a*c,
	e*b + d*e + b*a + a*d,
	e*c + d*a + c*d + a*e,
	e*d + d*b + c*e + b*c,
	b*a*b*a + a*b*a*b,
	c*a*c - b*c*a - b*a*b + a*b*c,
	d*b*d - c*d*b - c*b*c + b*c*d,
	c*b*a - b*c*b - b*a*c + a*c*a,
	d*b*a + d*a*d - b*a*b - a*d*b,
	d*a*c - c*b*d - c*a*b - a*d*a - a*c*d,
	d*a*b + c*d*b + c*b*c - b*d*a - b*c*d - a*b*a,
	c*b*c*a - c*a*b*c + b*a*c*b - a*c*a*b,
	c*b*c*b + b*c*b*c,
	d*a*d*a + a*d*a*d,
	d*b*c*a - d*a*d*b + c*d*b*c + b*c*b*d + b*c*a*b,
	d*b*c*b + c*a*b*a + b*a*d*b + b*a*b*c + a*d*b*c + a*b*a*d,
	d*b*c*d + c*a*d*b + c*a*b*c + b*d*b*c + a*b*d*b,
	c*a*b*c*a - b*a*b*c*b + a*c*a*b*a + a*b*a*b*c,
	c*a*b*c*b + c*a*b*a*c - b*c*a*b*a - a*b*a*c*a,
	c*a*b*a*b + b*a*b*c*b - a*c*a*b*a - a*b*a*b*c,
	c*a*b*a*c*b - b*a*c*a*b*a - a*b*a*c*a*b,
	c*a*b*a*c*a*b - b*c*a*b*a*c*a,
*/

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

#define MAXSIZE 50 // max dim for V
#define MAXGROBNER 800 // max number of grobner relations
#define MAXDEG 40 // max degree in grobner relations

int cardx; // size of X (basis of V)
char rack[MAXSIZE];  // the set X
char *leading[MAXGROBNER]; // leading terms of grobner basis
int grobnernr; // number of grobner relations
int pol[MAXDEG+2]; // Hilbert pol.
int md; // max degree for producebasis
char hilbert[MAXDEG*MAXSIZE*10]; // Hilbert pol for maple
int maxdeg=MAXDEG; // max degree, can be smaller than MAXDEG when -d is used in input
char belmt[MAXDEG]; // bound element. The result will be those elements below it, if it is not \0

isinstr(char c, char *str) {
	int i=0;
	while ((str[i] != c) && str[i])
		i++;
	if (str[i])
		return (i+1);
	else
		return (0);
}



readgrobner(FILE *in) {
	int samerel;
	int j,p,k;
	char c;
	grobnernr=0;
	while (1) {
		while (!isinstr(c=getc(in), rack) && c != EOF) {
			if (c=='%')
				while (getc(in) != '\n');
		}
		if (c==EOF) {
			fclose(in);
			return(0);
		}
		leading[grobnernr]=calloc(maxdeg, sizeof(char));
		leading[grobnernr][0]=c;
		j=1;
		samerel=1;
		while (samerel) {
			while(!isinstr(c=getc(in), rack) && !isinstr(c, "^+-,;"));
			if (isinstr(c, rack)) {
				leading[grobnernr][j++]=c;
			}
			else {
				switch (c) {
					case '^':
						fscanf(in, "%d", &p);
						for (k=1; k<p; k++)
							leading[grobnernr][j+k-1]=leading[grobnernr][j-1];
						j += p-1;
						break;
					case '+':
					case '-':
						samerel=0;
						while (!isinstr(c=getc(in),",;") && (c != EOF));
						break;
					case ',':
					case ';':
						samerel=0;
						break;
				}
			}
		}
//		leading[grobnernr][j+1]='\0';
//	printf("Read rel %d with leading term %s\n", grobnernr, leading[grobnernr]);
		grobnernr++;
	}
}


printleadings() {
	int i;
	for (i=0; i<grobnernr; i++)
		printf("Rel %d: %s\n", i, leading[i]);
}


length(char *a) {
	int i=0;
	while(a[i++]);
	return(i-1);
}

addinttostr(char *str, int i) {
	int j, k;
	while (*str++);
	str--;
	for (j=10, k=1; i/j; j*=10, k++);
	*(str+k) = '\0';
	for (; k>0; k--) {
		*(str+k-1) = i % 10 + '0';
		i /= 10;
	}
}


addstrtostr(char *str, char *add) {
	while (*str++);
	str--;
	while (*str++ = *add++);
}


deletechars(char *str, int i) {
	while (*str++);
	*(str-i-1) = '\0';
}

substr(char *a, char *b) {
	int i, j, la, lb;
	la=length(a);
	lb=length(b);
	for (i=0; i<=lb-la; i++) {
		for (j=0; j<la; j++) {
			if (a[j] != b[i+j])
				break;
		}
		if (j==la)
			return(i+1);
	}
	return(0);
}

iszero(char *mon, int deg) {
	int i,l;
	for (i=0; i<grobnernr; i++) {
		l=length(leading[i]);
		if (l>deg)
			continue;
		else {
			if (substr(leading[i], mon))
				return (i+1);
		}
	}
	return(0);
}



producebasis(char *mon, int deg) {
	int i,j;
	char c;
	if (deg>=maxdeg)
		return(0);
	for (i=cardx-1; i>=0; i--) {
		if (!*belmt && i>=isinstr(*mon, rack))
			continue;
		mon[deg]=rack[i];
		if (!iszero(mon, deg+1) && (!*belmt || strcmp(mon, belmt)<0)) {
			printf("(%2d)  %s\n", deg+1, mon);
			pol[deg+1]++;
			md=(deg>=md)?deg+1:md;
//	scanf("%c", &c);
			mon[deg+2]='\0';
			producebasis(mon, deg+1);
		}
		if (deg==0) {
			addstrtostr(hilbert, "(");
			for (j=0; j<md; j++) {
				printf("%d t^%d + ", pol[j], j);
				addinttostr(hilbert, pol[j]);
				addstrtostr(hilbert, "*t^");
				addinttostr(hilbert, j);
				addstrtostr(hilbert, " + ");
			}
			printf("%d t^%d\n", pol[md], md);
			addinttostr(hilbert, pol[j]);
			addstrtostr(hilbert, "*t^");
			addinttostr(hilbert, j);
			addstrtostr(hilbert, ")*");
			for (j=1; j<=md; j++)
				pol[j]=0;
			pol[0]=1;
			md=0;
		}
	}
	mon[deg]='\0';
}




main (int argc, char **argv) {
	int i,j;
	char c, mon[maxdeg];
	FILE *in;
	*belmt='\0';
	if (argc<2) {
		printf("Usage: %s [-d deg] [-b elmt] basis-file\n-d for max degree used (default=40)\n-b for printing all elements before elmt\nWarning: -d goes always before -b\n", argv[0]);
		abort();
	}
	if (strstr(argv[1], "-d")) {
		if (!strcmp(argv[1], "-d")) {
			if (!sscanf(argv[2], "%d", &maxdeg)) {
				printf("When using -d you must provide a number\n");
				abort();
			}
			argv += 2;
		}
		else {
			if (!sscanf(argv[1]+2, "%d", &maxdeg)) {
				printf("When using -d you must provide a number\n");
				abort();
			}
			argv += 1;
		}
	}
	if (strstr(argv[1], "-b")) {
		if (!strcmp(argv[1], "-b")) {
			strcpy(belmt, argv[2]);
			argv += 2;
		}
		else {
			strcpy(belmt, argv[1]+2);
			argv += 1;
		}
/*sac*/	printf("Bound: %s\n", belmt);
	}
	in = fopen(argv[1],"r");
	if (in==NULL) {
		printf("Can't open basis file %s\n", argv[1]);
		abort();
	}
	c=getc(in);
	if (c=='%') {
		while ((c=getc(in)) != '\n');
		c=getc(in);
	}
	rack[i=0]=c;
	while ((rack[++i]=getc(in)) != '\n');
	rack[i]='\0';
	cardx=i;
//	printf("Read card: size=%d, rack={%s}\n", cardx, rack);
	readgrobner(in);
//	printleadings();
	mon[0]=rack[cardx-1]; mon[1]='\0';
	pol[0]=1; hilbert[0]='\0';
	for (i=1; i<=maxdeg+1; i++)
		pol[i]=0;
	producebasis(mon, 0);
	deletechars(hilbert, 1);
	printf("Hilbert polynomial:\n%s\n", hilbert);
}