(* ---------------- Subring Intersection ---------------- *)
(* Version       : 1.0
Last modified : 12.1.98

Written by : 

Thomas Bayer, Institut fuer Informatik, Technische Universitaet Muenchen
www:    http://wwwmayr.informatik.tu-muenchen.de/personen/bayert/
email : bayert@in.tum.de *)
		


(*
	 vars = variable list of the polynomial ring 
	  
	 InvRIngIntersection[base1_List, base2_List, vars_List, p_Integer:0] 
	  base1, base2 : Basis for the invariant rings, 
  p = characteristic of the groun field
			Stopps after MAXUB iterations

  InvRingIntersection[base1_List, base2_List, vars_List, p_Integer:0, 
        ub_Integer] 
	  base1, base2 : Basis for the invariant rings, 
  p = characteristic of the groun field
		 Stopps after ub iterations
				
		SubAlgIntersection[base1_List, base2_List, vars_List, p_Integer, 
        ub_Integer] 	
			 base1, base2 : Basis for the subrings, 
  p = characteristic of the groun field
		  Stopps after ub iterations
*)

(* ----------------------------------------------------------------------------------------- *)
							
INVDEBUG = Table[True, {i, 10}];
DBSUBRING = 1;
MAXUB = 10;

Verbose[on] := Block[{},INVDEBUG = Table[True, {i, 10}];];
Verbose[off] := Block[{},INVDEBUG = Table[False, {i, 10}];];
Verbose[int] := Block[{},
		INVDEBUG = Table[False, {i, 10}];
		INVDEBUG[[DBSUBRING]] = True;
];

InvRingIntersection[base1_List, base2_List, vars_List, p_Integer:0] := 
  InvRingIntersection[base1, base2, vars, p,MAXUB] 

InvRingIntersection[base1_List, base2_List, vars_List, p_Integer, 
    ub_Integer] := 
  Module[{vars1, vars2, n, gb1, gb2, I1, I2,int, gens1, gens2,exit, rr1, rr2, 
      u, v, count},
		If[INVDEBUG[[DBSUBRING]], Print["Intersect ", base1, " and ", base2]];
		I1 = base1;
		I2 = base2;
		rr1 = Table[u[i] -> base1[[i]], {i, 1, Length[base1]}];
		rr2 = Table[v[i] -> base2[[i]], {i, 1, Length[base2]}];
		vars1 = Join[vars, Table[u[i], {i, 1, Length[base1]}]];
		vars2 = Join[vars, Table[v[i], {i, 1, Length[base2]}]];
		gb1 = GroebnerBasis[Table[base1[[i]] - u[i], {i, 1, Length[base1]}], vars1, 
        Modulus->p];
		gb2 = GroebnerBasis[Table[base2[[i]] - v[i], {i, 1, Length[base2]}], vars2, 
        Modulus->p];
		oldInt = {}; count = 0;
		exit = False; i = 0;
		While[Not[exit],
					count++; If[count == ub, exit = True]; 
					int = IdealIntersection[I1, I2, vars, p];
					If[INVDEBUG[[DBSUBRING]],
									 Print["-- Loop # ", count, " -------------"];
									Print[" Intersection = ", int];
									PrintDeg[int, vars];
					];
					If[SameQ[int, oldInt], 
									exit = True,
				
									I1 = GroebnerBasis[Join[int, gb1], vars1, Modulus->p];
									I1 = 
          Select[I1, SameQ[Intersection[Variables[#], vars], {}] &] /. rr1;
									I2 = GroebnerBasis[Join[int, gb2], vars2, Modulus->p];
									I2 = 
          Select[I2, SameQ[Intersection[Variables[#], vars], {}] &] /. rr2;
									If[INVDEBUG[[DBSUBRING + 1]],
									    Print[" new I1 = ", Expand[I1]];
											PrintDeg[Expand[I1], vars];
											If[INVDEBUG[[DBSUBRING + 3]], 
														
            Do[ Print[AlgebraMember[base1, I1[[k]], vars, p]], {k, 1, 
                Length[I1]}];
											];
									    Print[" new I2 = ", Expand[I2]];
											PrintDeg[Expand[I2], vars];
											If[INVDEBUG[[DBSUBRING + 4]], 
														
            Do[ Print[AlgebraMember[base2, I2[[k]], vars, p]], {k, 1, 
                Length[I2]}];
											];
									];
									oldInt = int;
				];
	 	];
		int
]

SubAlgIntersection[base1_List, base2_List, vars_List, p_Integer, 
    ub_Integer] := 
  Module[{B, vars1, vars2, n, gb1, gb2, I1, I2,int,newelems, exit, rr1, rr2, 
      u, v, count, stable, st1, st2, Itemp},
		If[INVDEBUG[[DBSUBRING]], Print["Intersect ", base1, " and ", base2]];
		n = Length[var]; st1 = st2 = {}; t21 = 0;
		I1 = base1;
		I2 = base2;
		rr1 = Table[u[i] -> base1[[i]], {i, 1, Length[base1]}];
		rr2 = Table[v[i] -> base2[[i]], {i, 1, Length[base2]}];
		vars1 = Join[vars, Table[u[i], {i, 1, Length[base1]}]];
		vars2 = Join[vars, Table[v[i], {i, 1, Length[base2]}]];
		gb1 = GroebnerBasis[Table[base1[[i]] - u[i], {i, 1, Length[base1]}], vars1, 
        Modulus->p];
		gb2 = GroebnerBasis[Table[base2[[i]] - v[i], {i, 1, Length[base2]}], vars2, 
        Modulus->p];
		oldInt = {}; count = 0; B = {};
		exit = False; i = 0;
		While[Not[exit],
					count++; If[count == ub, exit = True]; 
					int = IdealIntersection[I1, I2, vars, p];
					If[INVDEBUG[[DBSUBRING]],
									 Print["-- Loop # ", count, " -------------"];
									Print[" Intersection = ", int];
					];
					If[SameQ[int, oldInt], 
									exit = True,
				
					stable = Intersection[int, oldInt];
					If[Not[SameQ[stable, {}]], 
          Print[" Stable : ", stable, " exit = ", exit];
					d = Min[Table[Deg[stable[[i]], vars], {i, 1, Length[stable]}]];
					minI = Min[Table[Deg[I[[i]], vars], {i, 1, Length[I]}]];
					If[d<=minI, 
					stable = 
              Table[If[Deg[stable[[i]], vars] == d, stable[[i]], 0, 0], {i, 
                  1, Length[stable]}];
					Itemp = Complement[int, stable];
					Print[Table[Deg[stable[[i]], vars], {i, 1, Length[stable]}]," : ",  
              Table[Deg[Itemp[[i]], vars], {i, 1, Length[Itemp]}]];
									If[
              Min[Table[Deg[stable[[i]], vars], {i, 1, Length[stable]}]] < 
                Min[Table[Deg[Itemp[[i]], vars], {i, 1, Length[Itemp]}]],
										Print["Remove stable, ", stable];
									B = Union[B, stable];
									newelems  = 
                Flatten[
                  Table[vars[[i]] stable[[j]], {i, 1, Length[vars]}, {j, 1, 
                      Length[stable]}]];
									int = 
                GroebnerBasis[Union[Complement[int, stable], newelems], vars, 
                  Modulus->p];
								Print[" new int ",int];
						], Print[" Stable element not of minimal degree !!"]]; (* 
            end d<=minI *) 
					]; (* end stable*)
									I1 = GroebnerBasis[Join[int, gb1], vars1, Modulus->p];
									I1 = 
          Select[I1, SameQ[Intersection[Variables[#], vars], {}] &] /. rr1;
									I2 = GroebnerBasis[Join[int, gb2], vars2, Modulus->p];
									I2 = 
          Select[I2, SameQ[Intersection[Variables[#], vars], {}] &] /. rr2;
									If[INVDEBUG[[DBSUBRING]],
									    Print[" new I1 = ", Expand[I1]];
									    Print[" new I2 = ", Expand[I2]];
									];
									oldInt = int;
				];
	 	];
		B
]



AlgContQ[ poly_,algBase_List, vars_List, p_Integer:0] := 
  Module[{vars1, n, gb1,I1,int,  rr1, },
		n = Length[vars];
		rr1 = Table[u[i] -> algBase[[i]], {i, 1, Length[algBase]}];
		vars1 = Join[vars, Table[u[i], {i, 1, Length[algBase]}]];
		gb1 = GroebnerBasis[Table[algBase[[i]] - u[i], {i, 1, Length[algBase]}], 
        vars1, Modulus->p];
		PolynomialReduce[poly, gb1, vars1]
]
AlgebraMember[base_List, f_, var_List, p_Integer:0]:= Module[{gb, rr},
		n = Length[var];
		rr = Table[u[i] -> base[[i]], {i, 1, Length[base]}];
		vars = Join[var, Table[u[i], {i, 1, Length[base]}]];
		gb = GroebnerBasis[Table[base[[i]] - u[i], {i, 1, Length[base]}], vars1, 
        Modulus->p];
		PolynomialReduce[f, gb, vars][[2]]
	]
PrintDeg[base_List, vars_List] := Block[{},
			Print[" Degrees : ", Table[Deg[base[[i]], vars], {i, 1, Length[base]}]];
		]

IdealIntersection[base1_List, base2_List, var_List, p_Integer:0] := 
  Module[{int,t},
		int = GroebnerBasis[Join[t base1, (1 - t) base2], Join[{t}, var], 
        Modulus->p];
		Select[int, SameQ[Intersection[Variables[#], {t}], {}] &]
	]
		
NewElems[poly_, gb_List, vars_List, vv_List] := Module[{basis},
		    basis = GroebnerBasis[Join[{poly}, gb], vars];
		    Select[basis, SameQ[Intersection[Variables[#], vv], {}] &]
	]

ExtractPolynomials[B_List, gb_List, vars_List, vv_List] :=
  Complement[
    Select[B, 
      SameQ[Intersection[Variables[PolynomialReduce[#, gb, vars][[2]]], 
            vv], {}] &], {0}]

ElimBase[gb_List, vars_List, u_Symbol] := Module[{vv, base},
		vv = Join[vars, Table[u[i], {i, 1, Length[gb]}]];
		base = GroebnerBasis[Table[gb[[i]] - u[i], {i, 1, Length[gb]}], vv];
		{base, vv}
			]

ElimRules[gb_List, u_Symbol] := Table[u[i] -> gb[[i]], {i, 1, Length[gb]}]




(* Polynomials Tools *)
(* Version       : 0.1
Last modified : 2.12.96

Written by : Thomas Bayer, RISC, J. Kepler universtiy, Linz, Austria
email : Thomas.Bayer@risc.uni-linz.ac.at      *)

(* ----------------------------------------------------------------------- *)



PolyToFun[poly_, var_List] := Function[var, Evaluate[poly]]
PolyToFun[poly_] := PolyToFun[poly, Variables[poly]]

DegVec[mono_Times, vars_List:Variables[mono]]  := 
  Table[Exponent[mono, vars[[i]]], {i, Length[vars]}]
DegVec[mono_Symbol, vars_List:Variables[mono]] := 
  Table[Exponent[mono, vars[[i]]], {i, Length[vars]}]
DegVec[mono_Power, vars_List:Variables[mono]] := 
  Table[Exponent[mono, vars[[i]]], {i, Length[vars]}]
DegVec[mono_Integer, vars_List:Variables[mono]] := 
  Table[0, {i, Length[vars]}]
DegVec[mono_Rational, vars_List:Variables[mono]] := 
  Table[0, {i, Length[vars]}]
DegVec[mono_Real, vars_List:Variables[mono]] := Table[0, {i, Length[vars]}]
DegVec[mono_Complex, vars_List:Variables[mono]] := 
  Table[0, {i, Length[vars]}]
(*Deg[poly_] := Deg[poly, Variables[poly], PTotal]
Deg[poly_, vars_List] := Deg[poly, vars, PTotal]
Deg[poly_, Ordering_Symbol] := Deg[poly, Variables[poly], Ordering] *)

Deg[poly_, var_List:{}, Ordering_:PTotal] := Module[{pList, vars},
    If[SameQ[var, {}], vars = Variables[poly], vars = var];
    pList = Expand[poly];
    If[Not[SameQ[pList[[0]], Plus]], 
	pList = {DegVec[pList, vars]}, 
    	pList[[0]] = List; 
    	pList = Map[DegVec[#, vars] &, pList];
    ];
    Apply[Plus, Sort[pList, Ordering][[1]]]
]

HeadTerm[poly_Plus, vars_List, PTotal] := Module[{pList, deg},
    deg = Deg[poly, vars];
    pList = Expand[poly];
    pList[[0]] = List;
    Apply[Plus, Select[pList, Deg[#, vars] == deg &]]
]

HeadTerm[poly_Plus, vars_List, Ordering_:PLex] := Module[{pList},
    pList = Expand[poly];
    pList[[0]] = List; 
    Sort[pList, Ordering[#1, #2, vars] &][[1]]
]

HeadTerm[poly_Times, vars_List, Ordering_:PTotal] := Module[{pList},
    pList = Expand[poly];
    If[ SameQ[pList[[0]], Times], Return[poly]];
    pList[[0]] = List; 
    Sort[pList, Ordering[#1, #2, vars] &][[1]]
]

HeadTerm[poly_Power, vars_List, Ordering_:PTotal] := Module[{pList},
    pList = Expand[poly];
    If[ SameQ[pList[[0]], Power], Return[poly]];
    pList[[0]] = List; 
    Sort[pList, Ordering[#1, #2, vars] &][[1]]
]

HeadTerm[poly_Symbol, vars_List, Ordering_:PTotal]   := poly
HeadTerm[poly_Integer, vars_List, Ordering_:PTotal]  := poly
HeadTerm[poly_Rational, vars_List, Ordering_:PTotal] := poly
HeadTerm[poly_Real, vars_List, Ordering_:PTotal]     := poly
HeadTerm[poly_Complex, vars_List, Ordering_:PTotal]  := poly

HeadCoeff[poly_, vars_List, Ordering_:PTotal] := 
  Block[{hTerm = HeadTerm[poly, vars, Ordering], dv},
    dv = DegVec[hTerm, vars];
    hTerm / Product[vars[[i]]^dv[[i]], {i, Length[dv]}]
]

HeadMono[poly_, vars_List, Ordering_:PTotal] := 
	HeadTerm[poly, vars, Ordering] / HeadCoeff[poly, vars, Ordering]
					

(* True : l1 >= l2, False : l1 < l2	*)

PLex[l1_List, l2_List] := Block[{i = 1, exit = False, len = Length[l1]}, 
    While[!exit,
   	If[l1[[i]] > l2[[i]], Return[True]];
	If[l1[[i]] < l2[[i]], Return[False]];
        i++;
	If[i > len, Return[True]];
    ]
]

PTotal[l1_List, l2_List] := Apply[Plus, l1] >= Apply[Plus, l2]

PLex[m1_, m2_, vars_] := PLex[DegVec[m1, vars], DegVec[m2, vars]]
PTotal[m1_, m2_, vars_] := PTotal[DegVec[m1, vars], DegVec[m2, vars]]

LM := HeadMono
HM := HeadMono
LT := HeadTerm
HT := HeadTerm