//////////////////////////////////////////////////////////////////////
// Computation of Invaraiants up to a given degree
// Preliminary version (8.11.2000)
// Thomas Bayer, Institut fuer Informatik, Technische Universitaet Muenchen
// www:    http://wwwmayr.informatik.tu-muenchen.de/personen/bayert/
// email : bayert@in.tum.de


//////////////////////////////////////////////////////////////////////
/*
PRODECURES:
        AlgContainedQ(f,alg, index)     test if f is contained in K[alg] 
        GroupEquations(G)               compute equations of the matrix group G
        HomogenizeAction(action, index) homogenize 'action' w.r.t. var(index) 
        HomQ(action)                    test if 'action' is homogeneous
        Invariants(G, action[,opt])     invariants of G w.r.t. action 
        Relations(I)                    ideal for algebraic relations

        ContainedQ(data, f)             test if f is contained in 'data'

*/
//////////////////////////////////////////////////////////////////////

LIB "elim.lib";

int MINGEN = 1;
int MAXINVARS = 20;     // maximal 20 invariants for elimination
printlevel = 1;

//////////////////////////////////////////////////////////////////////

proc Invariants(ideal G, ideal action, list #)
"USAGE:   Invariants(G, action [, degrees]); ideal G, action; int/list degrees.
PUROPSE: 
ASSUME:  
RETURN:  
EXAMPLE: example 
"
{
	int d, dn, degAction, i, j, dimV, nrVars, gVars, nrInv, loop, invQ, sizeJ, homQ, ub;
        intmat dv;
        list degList;
        poly mp;
        string ringSTR1, ringSTR2, addSTR;
	
        // ringSTR1 for homogenization, ringSTR2 for elimination

        int dbPrt = printlevel-voice+2;

        // initialize the list of degrees

        if(size(#) == 1) {
                ub = #[1]; 
                for(i = 1; i <= ub + 1; i = i + 1) { degList[i] = i;}
        }
        else {
                degList = #;
                ub = #[size(#)];
                degList[size(degList) + 1] = ub + 1;
        }

        // prepare a new ring 'RIR' and homogenize the action and the group
        // the ring 'RIS' will be used for elimination

	nrVars = nvars(basering);
        gVars = nrVars - ncols(action); 
        dimV = ncols(action);
        homQ = HomQ(action);
   
        mp = minpoly;
        ringSTR1 = "ring RIR = (" + charstr(basering) + "),("; 
        ringSTR2 = "ring RIS = (" + charstr(basering) + "),(";
        if(homQ == 1) { ringSTR1 = ringSTR1 + varstr(basering) + ", ";}
        else {
                for(i = 1; i <= gVars; i++) {
                        ringSTR1 = ringSTR1 + string(var(i)) + ", ";
                }      
                ringSTR1 = ringSTR1 + "s(" + string(gVars + 1) + "), "; 
                gVars = gVars + 1;
                for(i = gVars; i <= nrVars; i++) {
                        ringSTR1 = ringSTR1 + string(var(i)) + ", ";  
                }  
                nrVars = nrVars + 1; 
        }
        ringSTR1 = ringSTR1 + "X), dp;";   
        
        def RIB1 = basering;
	def RIB = basering; export(RIB); 
      

        execute(ringSTR1);
        if(npars(basering) > 0) { minpoly = number(imap(RIB, mp)); }


        for(i = gVars + 1; i <= nrVars; i++) {
                ringSTR2 = ringSTR2 + string(var(i)) + ", ";
        }   
        ringSTR2 = ringSTR2 + "Y(1.." + string(MAXINVARS) + ")), (dp(" + string(nrVars - gVars);
        ringSTR2 = ringSTR2 + "), dp(" + string(MAXINVARS) + "));";        
  
        execute(ringSTR2); export(RIS);
        ideal alg;
        poly g;

        setring(RIR);
        ideal action, alg, B, Bn, G, Gh, I, In, J, Jn, newJ;
        matrix coMx;
        poly g;

        action = imap(RIB, action);
        G = imap(RIB, G);
        if(homQ == 0) {
                action = HomogenizeAction(action, gVars); 
                G = G, var(gVars) -  1; 
        }
        Gh = homog(std(G), var(nrVars + 1));  

        //print(" G = " + string(Gh));
        //print(" action = " + string(action));

        degAction = deg(action[1]);

        map F;
        map Fd = RIR, maxideal(1); 
        map Fn = RIR, maxideal(1);
      
        Fd[nrVars + 1] = 1;      // dehomogenize;
        Fn[nrVars + 1] = 0;      // delete
        for(i = 1; i <= gVars; i++) { Fd[i] = 0;F[i] = var(i);} 
                    
        // 
	// compute first gb to intersect I with R1

        d = degList[1]; 

        I = action^d, Gh;
        nrInv = 0;
        j = 1;
        loop = 1;

        while(loop == 1) {
                dbprint(dbPrt, "-- loop, d = " + string(d) + ", ub = " + string(ub) + " degBound = " + string(degAction * d) + " ---------");
                In = I, Gh;
                I = 0;                  // save memory
                degBound = degAction * d;   // correct ?    
                dbprint(dbPrt, " std of size = " + string(size(In)));
                Jn = std(In);
                In = 0;                 // save memory
                dbprint(dbPrt, " nselect , size = " + string(size(Jn)));
                J = nselect(Jn, 1, gVars);  
                Jn = 0;                         // save memory    
                dbprint(dbPrt, " size of J = " + string(ncols(J))); 
                dbprint(dbPrt, string(maxdeg(J))); 
                dbprint(dbPrt - 2, " J = " + string(J));
                dbprint(dbPrt, " -- test elements -- ");
                sizeJ = ncols(J);

                if(size(J) > 0) { 

                Bn = 0;
                // test all elements of the groebner basis of J 
                for(i = 1; i <= sizeJ; i++) {
                        g = J[i];  
                        invQ = 0;
                        if(deg(g) == d * degAction) {
                                invQ = 1;
                                g = Fd(g);      // dehomogenize
                                if(MINGEN == 1){
                                        setring(RIS);
                                        g = imap(RIR,g);  
                                        if(!AlgContainedQ(g, alg, dimV)) {
                                                nrInv = nrInv + 1;
                                                alg = alg, g - Y(nrInv);
                                                alg = std(alg);
                                                setring(RIR);
                                                B = B, g;
                                                addSTR = " added ";
                                        }
                                        else { setring(RIR); addSTR = "";}
                                }
                                else {
                                        if(!ContainedQ(g, B)) {B = B, g; }
                                }
                        } 
                        if(invQ == 1) {
            dbprint(dbPrt, string(i) + ": deg = " + string(deg(g)) + ", " + string(g) + " is invariant, " + addSTR);
                        }
                        else {dbprint(dbPrt, string(i) + ": deg = " + string(deg(g)) + ", " + string(g));}
                }
                dbprint(dbPrt, " invariants = " + string(nrInv));
                B = simplify(B, 2);
                } // end if(size(Jn) > 0)
                j = j + 1;
                d = degList[j];
                if(d <= ub) { I = action^d; }
                else { loop = 0;} 
	}
        setring(RIB1);
        ideal B = simplify(imap(RIR,B), 2);
        kill(RIR);
        kill(RIS);
        kill(RIB);
        degBound = 0;
	return(B);
}

////////////////////////////////////////////////////////////////////

proc HomogenizeAction(ideal action, int index)
"USAGE:    HomogenizeAction(action, index); ideal action, int index 
PUROPSE: homogenize each polynomial of 'action' w.r.t. var(index)
         s.t. all polynomials of 'action' are of the same degree 
ASSUME:  
RETURN:  
EXAMPLE: example 
"
{
        ideal hAction;
        int d, i, j;
        matrix coMx;
        poly f, p;

        d = maxdeg1(action);

        p = 1;
        for(i = 1; i <= nvars(basering); i++) { p = p*var(i);}
        for(i = 1; i <= size(action); i++) {
                coMx = coef(action[i], p);
                f = 0;      
                for(j = 1; j <= ncols(coMx); j++) {
                        f = f + coMx[1,j]*coMx[2,j]*var(index)^(d - deg(coMx[1,j]));
                }
                hAction[i] = f;
        }
        return(hAction);
}

proc HomQ(action)
{
        int d, i, homQ;
        intvec w = qhweight(action);

        homQ = 1;
        for(i = 1; i <= nvars(basering); i++) { if(w[i] != 1) { homQ = 0;}}

        if(homQ == 1) {
                d = maxdeg1(action);
                for(i = 1; i <= size(action); i++) { if(deg(action[i]) != d) { homQ = 0;}}
        }
        return(homQ);
}
      
////////////////////////////////////////////////////////////////////
// Auxillary stuff


proc AlgContainedQ(poly f, ideal algG, int d)
{
        poly f1 = reduce(f, algG); 
        ideal J = nselect(ideal(f1),1,d);
        return(size(J != 0));
}

////////////////////////////////////////////////////////////////////

proc GroupEquations(list G)
"USAGE:   GroupEquations(G); list G
PUROPSE: compute the ideal of the finite group G
ASSUME:  G is a list of n x n matrices, basering contains at least 
         n^2 variables
RETURN:  ideal
EXAMPLE: example 
"
{
        int i, j, k, n; 
        ideal I, J, action;

        n = ncols(G[1]);        // number of variables for G

        // M[i][j] -> var((i - 1) * n + j)

        for(i = 1; i <= n; i++) {
                action[i] = 0;
                for(j = 1; j <= n; j++) {
                        I = I, var((i - 1) * n + j) - G[1][i, j];
                        action[i] = action[i] + var((i - 1) * n + j) * var(n^2 + j);
                }
        }
        print(" action = " + string(action));
        for(k = 2; k <= size(G); k++) {
                J = 0;
                for(i = 1; i <= n; i++) {
                        for(j = 1; j <= n; j++) {
                                J = J, var((i - 1) * n + j) - G[k][i, j];
                        }
                }
                I = std(intersect(I, J));
        }
        return(I);
}

////////////////////////////////////////////////////////////////////

proc Relations(I, int d)
"USAGE:   Relations(I); list/ideal I, int d
PUROPSE: ideal for computing relations
ASSUME:  
RETURN:  ideal
EXAMPLE: example 
"
{
        int i; 
        ideal J;

        for(i = 1; i <= size(I); i++) {
                J[i] = I[i] - var(d + i);
        }
        return(J);
}

proc ReduceAndTest(I,alg,int d)
"USAGE:   Relations(I); list/ideal I, int d
PUROPSE: ideal for computing relations
ASSUME:  
RETURN:  ideal
EXAMPLE: example 
"
{
        int i; 
        ideal J,J1;
        poly f;

        for(i = 1; i <= size(I); i++) {
                f = reduce(I[i], alg);
                J = J, select(ideal(f),1,d);
        }
        return(simplify(J,2));
}

       
////////////////////////////////////////////////////////////////////

proc Get(B, list index)
// forget elements of B which do not contain any of the variables var(lb..ub)
{
        int i;
        ideal Bn;
        matrix ma;
        poly p = 1;

        for(i = 1; i <= size(index); i++) { p = p*var(index[i]);} print(p);

        for(i = 1; i <= size(B); i++) {
                ma = coef(B[i],p);
                if(ma[1,1] != 1) { Bn = Bn, B[i]; }
        }
        return(simplify(Bn, 2));
} 


/////////////////////////////////////////////////////////////////

proc Lead(ideal I)
{
        for(int i = 1; i <= ncols(I); i++) {print(string(i) + ": " + string(lead(I[i])));};
} 

/////////////////////////////////////////////////////////////////
// from 'zeroset.lib'

proc ContainedQ(data, f, list #)
"USAGE:    ContainedQ(data, f [, opt]); list data; f is of any type, int opt
PUROPSE:  test if 'f' is an element of 'data'.  
RETURN:   int
          0 if 'f' not contained in 'data'
          1 if 'f' contained in 'data'
OPTIONS:  opt = 0 : use '==' for comparing f with elements from data
          opt = 1 : use 'SameQ' for comparing f with elements from data
"
{
	int opt, i, found;
	if(size(#) > 0) { opt = #[1];}
	else { opt = 0; }
	i = 1;
	found = 0;

	while((!found) && (i <= size(data))) {
		if(opt == 0) {
			if(f == data[i]) { found = 1;}
			else {i = i + 1;}
		}
		else {
			if(SameQ(f, data[i])) { found = 1;}
			else {i = i + 1;}
		}
	}
	return(found);
}