// Intersection of Invaraiant Rings.
// Preliminary version (1.9.2000)
// Implementation by : Thomas Bayer, Institut fuer Informatik, Technische Universitaet Muenchen
// email : bayert@in.tum.de

// Intersection1(ideal inv1, ideal inv2, int ub)
// Intersection2(ideal alg1, ideal alg2, int ub)
// IntersectionOld(ideal inv1, ideal inv2);

// some algorithms are taken from 'rinvar.lib'
//

LIB "elim.lib";


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

int maxLoops = 50;
printlevel = 1;

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

proc Intersection1(ideal alg1, ideal alg2, int ub)
"USAGE:   Intersection1(alg1, alg2, ub); ideal alg1, alg2; int ub; 
PUROPSE: compute the intersection of the invariant rings generated by
         'alg1', 'alg2' respectively, where 'ub' is an upper bound
         for the generators of the intersection 'R'. 
RETURN:  ideal
ASSUME:  
NOTE:    the result is a generating set for the homogenous ideal 
         generated by the elements of R of positive degree.
         To obtain the generators of the invariant ring 'R' one
         has to apply the Reynolds operator to the generators. 
"
{
	int i, nrVars;
	intvec w, h1, w1, h2, w2;
	ideal I, Iold, J1, J2;

	nrVars = nvars(basering);
	def RIB = basering;
	for(i = 1; i <= nrVars; i++) { w[i] = 1;}

	// compute first gb to intersect I with R1

        degBound = ub;

	h1 = HilbertSeries(alg1, w);
	w1 = HilbertWeights(alg1, w);
     	def RIR1 = IdealAlgGB(alg1, h1, w1);

	setring RIR1;
	export(RIR1);
	ideal alg1, I1;
	alg1 = imap(RIB, alg1);
	map F1 = RIR1, maxideal(1);
	for(i = nrVars + 1; i <= nvars(RIR1); i++) { F1[i] = alg1[i - nrVars]; }

	// compute second gb to intersect I with R2

	setring(RIB);

	h2 = HilbertSeries(alg2, w);
	w2 = HilbertWeights(alg2, w);
	def RIR2 = IdealAlgGB(alg2, h2, w2);  

	setring RIR2;

	export(RIR2);
	ideal alg2, I2;
	alg2 = imap(RIB, alg2);
	map F2 = RIR2, maxideal(1);
	for(i = nrVars + 1; i <= nvars(RIR2); i++) { F2[i] = alg2[i - nrVars]; }
	
	setring(RIB);

	// prepare for the loop

	Iold = 0; 
 	I = jet(std(intersect(alg1, alg2)), ub); 

        i = 1;

	dbprint(dbPrt, "enter loop");

	while(!(SameQ(I, Iold)) && (i < maxLoops)) {
		 
		Iold = I;
		setring(RIR1); 
		
		I1 = IdealSubAlgIntersect(imap(RIB, I), imageid, F1, nrVars, h1, w1, ub);
                I1 = jet(I1, ub);
               
		setring(RIR2); 
	
		I2 = IdealSubAlgIntersect(imap(RIB, I), imageid, F2, nrVars, h2, w2, ub);  
	
                I2 = jet(I2, ub);
                setring(RIB);

	        J1 = jet(std(imap(RIR1, I1)), ub);
                J2 = jet(std(imap(RIR2, I2)), ub);
                
		I = intersect(J1, J2);
                I = std(jet(I, ub));
                I = jet(I, ub);  
		i = i + 1; 	
	}
        kill(RIR1);
        kill(RIR2);
	return(I);
}

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

proc Intersection2(ideal alg1, ideal alg2, int ub)
"USAGE:   Intersection2(alg1, alg2, ub); ideal alg1, alg2; int ub; 
PUROPSE: compute the intersection of the subrings of K[t1,...] 
         generated by 'alg1', 'alg2' respectively, where 'ub' is an 
         upper bound for the generators of the intersection 'R'. 
RETURN:  ideal
ASSUME:  
NOTE:    the result is a K-vectorspace basis of 'R<=ub'
"
{
	int d, i, j, nrVars, loop;
	intvec w, h1, w1, h2, w2;
	ideal I, B, SB;
        intmat dv;
        list idL;
        poly h;

	nrVars = nvars(basering);
	def RIB = basering;
	for(i = 1; i <= nrVars; i++) { w[i] = 1;}

        degBound = ub;

	// compute first gb to intersect I with R1

	
	h1 = HilbertSeries(alg1, w);
	w1 = HilbertWeights(alg1, w);
     	def RIR1 = IdealAlgGB(alg1, h1, w1);
	setring RIR1;
	export(RIR1);
	ideal alg1, I1;
	alg1 = imap(RIB, alg1);
	map F1 = RIR1, maxideal(1);
	for(i = nrVars + 1; i <= nvars(RIR1); i++) { F1[i] = alg1[i - nrVars]; }
        export(alg1); export(I1); export(F1);

	// compute second gb to intersect I with R2

	setring(RIB);

	h2 = HilbertSeries(alg2, w);
	w2 = HilbertWeights(alg2, w);
	def RIR2 = IdealAlgGB(alg2, h2, w2);  
	
	setring RIR2;

	export(RIR2);
	ideal alg2, I2;
	alg2 = imap(RIB, alg2);
	map F2 = RIR2, maxideal(1);
	for(i = nrVars + 1; i <= nvars(RIR2); i++) { F2[i] = alg2[i - nrVars]; } 
        export(alg2); export(I2); export(F2);
	setring(RIB);

	// prepare for the loop
        I = simplify(jet(intersect(alg1, alg2), ub), 2); 
        d = MinDeg(I);
        loop = 1;

        while(loop == 1) {
                idL = BasisElements1(RIR1, RIR2, RIB, I, ub);
                I = simplify(jet(idL[1], ub), 2);
                SB = idL[2];
                B = B, SB; B = simplify(B,2);
                for(i = 1; i <= ncols(I); i++) {
                        h = I[i];
                        if(ContainedQ(SB, h) == 1) { 
                                I[i] = 0;
                                I = I, NewElems(h);
                        }
                }
                I = simplify(I,2);
                d = MinDeg(I);
                if(d > ub || d == -1) { loop = 0;}
                else {I = std(I);}
	}

	kill(RIR1);
	kill(RIR2);
	return(B);
}

////////////////////////////////////////////////////////////////////
// no degree bound

proc IntersectionOld(ideal alg1, ideal alg2)
"USAGE:   Intersection1(alg1, alg2, ub); ideal alg1, alg2; int ub; 
PUROPSE: compute the intersection 'R' of the invariant rings of finite 
         groups generated by 'alg1', 'alg2'.
RETURN:  ideal
ASSUME:  'alg1', 'alg2' are generators of invariant rings of finite
         groups.
NOTE:    the result is a generating set for the homogenous ideal 
         generated by the elements of R of positive degree.
         To obtain the generators of the invariant ring 'R' one
         has to apply the Reynolds operator to the generators. 
"
{
	int i, nrVars, ub;
	intvec w, h1, w1, h2, w2;
	ideal I, Iold, J1, J2;

        ub = 20000;     // sufficient
	nrVars = nvars(basering);
	def RIB = basering;
	for(i = 1; i <= nrVars; i++) { w[i] = 1;}

	// compute first gb to intersect I with R1

	h1 = HilbertSeries(alg1, w);
	w1 = HilbertWeights(alg1, w);
     	def RIR1 = IdealAlgGB(alg1, h1, w1);

	setring RIR1;
	export(RIR1);
	ideal alg1, I1;
	alg1 = imap(RIB, alg1);
	map F1 = RIR1, maxideal(1);
	for(i = nrVars + 1; i <= nvars(RIR1); i++) { F1[i] = alg1[i - nrVars]; }

	// compute second gb to intersect I with R2

	setring(RIB);

	h2 = HilbertSeries(alg2, w);
	w2 = HilbertWeights(alg2, w);
	def RIR2 = IdealAlgGB(alg2, h2, w2);  

	setring RIR2;

	export(RIR2);
	ideal alg2, I2;
	alg2 = imap(RIB, alg2);
	map F2 = RIR2, maxideal(1);
	for(i = nrVars + 1; i <= nvars(RIR2); i++) { F2[i] = alg2[i - nrVars]; }
	
	setring(RIB);

	// prepare for the loop

	Iold = 0; 
 	I = std(intersect(alg1, alg2)); 

        i = 1;

	dbprint(dbPrt, "enter loop");

	while(!(SameQ(I, Iold)) && (i < maxLoops)) {
		 
		Iold = I;
		setring(RIR1); 
		
		I1 = IdealSubAlgIntersect(imap(RIB, I), imageid, F1, nrVars, h1, w1, ub);
               
		setring(RIR2); 
	
		I2 = IdealSubAlgIntersect(imap(RIB, I), imageid, F2, nrVars, h2, w2, ub);  
	;
                setring(RIB);

	        J1 = std(imap(RIR1, I1));
                J2 = std(imap(RIR2, I2));
                
		I = std(intersect(J1, J2));  // remove std if it is too slow 
		i = i + 1; 	
	}  
        kill(RIR1);
        kill(RIR2);
	return(I);
}

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

proc NewElems(poly f)
{
        ideal I;
        
        for(int i = 1; i <= nvars(basering); i++) { I[i] = var(i)*f; }
        return(I);
}

proc MinDeg(ideal I)
// Type : int
// Purpose : minimal integer of data
{
	int i; 
        intmat dv = mindeg(I);

	int min = dv[1, 1];

	for(i = 1; i <= size(dv); i = i + 1) {
		if(int(dv[1, i]) < min) { min = int(dv[1, i]); }
	}
	return(min);
}

proc BasisElements1(RIR1, RIR2, RIB, ideal J, int ub)
// ub .. upper bound
{
        setring RIB;
        ideal B, Iold, I, In, J1, J2;
       	int parttime, dbPrt;
	int d, dmin, i, j, nrVars, loop;
	intvec w, h1, w1, h2, w2;
        intmat dv;

	dbPrt = printlevel-voice+2;
	int runtime = rtimer; 

	nrVars = nvars(basering);
	for(i = 1; i <= nrVars; i++) { w[i] = 1;}

	// prepare for the loop

	Iold = 0; 
        I = simplify(jet(J, ub),2);
	parttime = rtimer;   
	dbprint(dbPrt - 1, "I = " + string(I)); 
	i = 1;
        d = MinDeg(I);
        if(d <= ub) { loop = 1;}
        else { loop = 0;}

	dbprint(dbPrt, "enter loop");

	while(loop == 1) {
		dbprint(dbPrt, "-- loop "+ string(i) + ", mindeg = " + string(d));
		dbprint(dbPrt - 1, " I = " + string(I));
		dmin = MinDeg(I); 
                dbprint(dbPrt, "size of I = " + string(size(I)) + " degrees = " + string(maxdeg(I)));
		dbprint(dbPrt, " dmin = " + string(dmin));

		setring(RIR1); 
		parttime = rtimer;
		I1 = IdealSubAlgIntersect(imap(RIB, I), imageid, F1, nrVars, h1, w1, ub);
		//itime[3] = itime[3] + rtimer - parttime;  
                if(MinDeg(I1) > ub) { loop = 0;}
             
		//dbprint(dbPrt - 1, "new I1 = " + string(I1));
		dbprint(dbPrt, "size of I1 = " + string(size(I1)) + " degrees = " + string(maxdeg(I1)));
                I1 = simplify(jet(I1, ub),2);
                dbprint(dbPrt, "size of I1 = " + string(size(I1)) + " degrees = " + string(maxdeg(I1)));
                if(loop == 1) {
		setring(RIR2); 
		parttime = rtimer;
		I2 = IdealSubAlgIntersect(imap(RIB, I), imageid, F2, nrVars, h2, w2, ub);  
		//itime[4] = itime[4] + rtimer - parttime;  
                I2 = simplify(jet(I2, ub),2);
                dbprint(dbPrt, "size of I2 = " + string(size(I2)) + " degrees = " + string(maxdeg(I2)));
                if(MinDeg(I2) > ub) { loop = 0;}
                }
		//dbprint(dbPrt - 1, "new I2 = " + string(I2));
		dbprint(dbPrt, "size of I2 = " + string(size(I2)) + " degrees = " + string(maxdeg(I2)));
		setring(RIB);
                if(loop == 1) {
		parttime = rtimer;

                J1 = jet(std(imap(RIR1, I1)), ub);
                J2 = jet(std(imap(RIR2, I2)), ub);
                
		I = simplify(jet(std(intersect(J1,J2)), ub), 2); 	
                
       		//itime[5] = itime[5] + rtimer - parttime;  
		i++;
                d = MinDeg(I);	
                dbprint(dbPrt, "minimal degree = " + string(d) + ", ub = " + string(ub));	
                dbprint(dbPrt, "size of I = " + string(size(I)) + " degrees = " + string(maxdeg(I)));
                if(d == dmin || d > ub) { loop = 0;}
                } // end loop == 0
                else { loop = 0;}
	} // end while

	dbprint(dbPrt, " result, I = " + string(I));
	//dbprint(dbPrt, " i = " + string(i) + ", maxLoops = " + string(maxLoops));
	kill(RIR1);
	kill(RIR2);
	//itime[6] = itime[6] + rtimer - runtime;   
        dbprint(dbPrt, "size of I = " + string(size(I)) + " degrees = " + string(maxdeg(I)));
        dbprint(dbPrt, " found " + string(size(jet(I,d))) + " elements of degree " + string(d));
        dbprint(dbPrt, "----------------------------------------------------");
	return(list(I, jet(I, dmin)));
}

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


proc IdealSubAlgIntersect(ideal I, ideal alg, map H, int nr, hs, w, int ub) 
{
	ideal In = I,alg;
        ideal J = std(In, w);
        ideal J1 = jet(nselect(J, 1, nr), ub, w);
	return(simplify(H(J1),2));
} 

proc SameQ(ideal I, ideal J)
{
	int i = 1;

	if(size(I) == size(J)) {
		while((i < size(J)) && (I[i] == J[i])) 
			{i++;};
	}
	return((I[i] == J[i]) && (size(I) == size(J)));
}
///////////////////////////////////////////////////////////////////////////////
// from 'rinvar.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);
}
///////////////////////////////////////////////////////////////////////////////

proc HilbertWeights(ideal I, wt)
"USAGE:     HilbertWeights(I, wt); ideal I, intvec wt
PUROPSE: compute the weights of the 'slack' varaibles needed for the computation
of the algebraic relations of the generators of I s.t. the Hilbert driven std
can be used.
RETURN: intvec 
ASSUME: basering = K[t_1,...,t_m], I is quasihomogenous w.r.t. wt
"
{
	int offset = size(wt);
	intvec wtn = wt;

	for(int i = 1; i <= size(I); i++) { wtn[offset + i] = deg(I[i], wt); }
	return(wtn);
}

///////////////////////////////////////////////////////////////////////////////
// The algorithms 'HilbertSeries' is taken from the 'rinvar.lib'.

proc HilbertSeries(ideal I, wt)
"USAGE:  HilbertSeries(I, wt); ideal I, intvec wt    
PUROPSE: compute the polynomial p of the Hilbert Series,represented by p/q, of 
the ring K[t_1,...,t_m,y_1,...,y_r]/I1 where 'wt' are the weights of the 
variables, computed, e.g., by 'HilbertWeights' 
RETURN:  poly
ASSUME:  I is homogenous w.r.t. wt and the ideal 'I1' is of the form
I[1] - y_1,...,I[r] - y_r.
"
{
	int i;
	intvec hs1;
	matrix coMx;
	poly f = 1;

	for(i = 1; i <= ncols(I); i++) { f = f * (1 - var(1)^deg(I[i], wt));}
	coMx = coeffs(f, var(1));
	for(i = 1; i <= deg(f) + 1; i++) { 
		hs1[i] = int(coMx[i, 1]);
	}
	hs1[size(hs1) + 1] = 0;
	return(hs1);
}

proc HilbertSeries1(wt)
"USAGE:  HilbertSeries(I, wt); ideal I, intvec wt    
PUROPSE: compute the polynomial p of the Hilbert Series represented by p/q of 
the ring K[t_1,...,t_m,y_1,...,y_r]/I where wt are the weights of the vars,
computed, e.g., by HilbertWeights 
RETURN:  poly
ASSUME:  I is homogenous
EXAMPLE:
"
{
	int i, j;
	intvec @hs1;
	matrix ma;
	poly f = 1;

	for(i = 1; i <= size(wt); i++) { f = f * (1 - var(1)^wt[i]);}
	ma = coef(f, var(1));
	j = ncols(ma);
	for(i = 0; i <= deg(f); i++) { 
		if(var(1)^i == ma[1, j]) {
			@hs1[i + 1] = int(ma[2, j]);
			j--;
		}
		else { @hs1[i + 1] = 0; }
	}
	@hs1[size(@hs1) + 1] = 0;
	return(@hs1);
}
///////////////////////////////////////////////////////////////////////////////

proc IdealAlgGB(ideal F, intvec hs, intvec w)
"USAGE:   IdealAlgGB(ideal I, F [, wt]);ideal I; F is a list/ideal, intvec wt.
PUROPSE: compute the Zariski closure of the image of the variety of I under 
the morphism F.
if I and F are quasihomogenous w.r.t. 'wt' then the Hilbert-driven 'std' is used.
RETURN:  ring [ideal imageid]  (variables have changed to Y(1..k))
	 imageid	ideal of the Zariski closure of F(X)
         imageid = ideal of the Zariski closure of F(X), X = variety of I.
NOTE:    This is the same function as 'ImageVariety' of the 'rinvar.lib'.
"
{
	int i, nrNewVars;
	string ringSTR1, order;
	
	def RARB = basering;
	nrNewVars = size(F);
	// create new ring for elimination, Y(1),...,Y(m) are slack variables

	poly mPoly = minpoly;
	order = "(dp(" + string(nvars(basering)) + "), dp);";
	ringSTR1 = "ring RAR1 = (" + charstr(basering) + "), (" + varstr(basering) + ", Y(1.." + string(nrNewVars) + ")), " + order;

	execute(ringSTR1); 
	minpoly = number(imap(RARB, mPoly));

	ideal J1, imageid, Fm;

	Fm = imap(RARB, F);

	// get the ideal of the graph of F : X -> Y and compute a standard basis

	for(i = 1; i <= nrNewVars; i++) { J1[i] = var(i + nvars(RARB)) - Fm[i];}  

       	J1 = std(J1, hs, w);	// Hilbert-driven algorithm 

	imageid = J1;

	export(imageid);
	return(RAR1);
}