Package $primary

export InitialIdeal;
export PrimaryPoincare;
export TgCone;
export IntersectionMultiplicity;

Define About()
  PrintLn "    package for: HilbertPoincare series for primary ideals,";
  PrintLn "              Tangent Cone, Initial Ideal.";
  PrintLn "    Authors : A.M.Bigatti, L.Robbiano";
  PrintLn "    Date    : 2000 - CoCoA-4";
  PrintLn "    Date    : 2014 - CoCoA-5";
EndDefine;

------[   help & examples   ]--------

Define Man()
  PrintLn "*** CoCoA-5 ***";
  PrintLn "Some functions in this package are now official functions.  See";
  PrintLn "    ? InitialIdeal";
  PrintLn "    ? PrimaryPoincare";
  PrintLn "    ? TgCone";
  PrintLn "    ? IntersectionMultiplicity";
  PrintLn "*** CoCoA-5 ***";
  PrintLn;
  PrintLn "*** CoCoA-4 manual for this package ***";  
  PrintLn "Suggested alias for this package:";
  PrintLn "    Alias P := $primary;";
  PrintLn "";
  PrintLn "SYNTAX";
  PrintLn "    PrimaryPoincare(I: Ideal, Q: Primary Ideal): TAGGED(\"PSeries\")";
  PrintLn "    TgCone(I: Ideal): Ideal";
  PrintLn "    InitialIdeal(I: Ideal, L: List): Ideal";
  PrintLn "    P.E(N: Int, HV: HVector): Int";
  PrintLn "";
  PrintLn "DESCRIPTION";
  PrintLn "Let P be a polynomial ring, M the maximal ideal generated by the";
  PrintLn "indeterminates. This package deals with the problem of computing the";
  PrintLn "Hilbert-Poincare' series of (P/I) / (Q+I)/I, where (Q+I)/I is a";
  PrintLn "primary ideal for M/I.";
  PrintLn "";
  PrintLn "    The main function is";
  PrintLn "            PrimaryPoincare(I, Q);";
  PrintLn "    which computes the Hilbert-Poincare' series of (Q+I)/I,";
  PrintLn "    where (Q+I)/I is a primary ideal for M/I.";
  PrintLn "";
  PrintLn "    Some related functions are also given, namely";
  PrintLn "            TgCone(I)";
  PrintLn "    which computes the initial ideal of I with respect to M;";
  PrintLn "            InitialIdeal(I, L)";
  PrintLn "    which computes the initial ideal of I with respect to the ideal";
  PrintLn "    generated by a subset L of the indeterminates;";
  PrintLn "            E(N, HV)";
  PrintLn "    which computes the (generalized) coefficients of the Hilbert";
  PrintLn "    Polynomial with respect to a particular binomial basis";
  PrintLn "";
  PrintLn ">EXAMPLES<";
  PrintLn "";
  PrintLn "    Alias P := $primary; -- proposed alias";
  PrintLn "";
  PrintLn "    Use S ::= QQ[x,y,z];";
  PrintLn "    I := ideal(x^3-y*z, y^2-x*z, z^2-x^2*y);";
  PrintLn "    TgCone(I);";
  PrintLn "    InitialIdeal(I, [y,z]);";
  PrintLn "    Q := ideal(y, z);";
  PrintLn "    PS := PrimaryPoincare(I, Q); PS;";
  PrintLn "    $hp.PSerToHilbert(PS); -- the Hilbert Function associated to PS";
  PrintLn "";
  PrintLn "    Use S ::= ZZ/(32003)[x,y,z,w];";
  PrintLn "    I := ***ideal(x^5 - yz, y^4 - xz^2, xy^3 - zw, x^2z - yw,";
  PrintLn "    y^2z^2 - w^3, y^3z - x^2w^2, x^3w - z^2, xyw^2 - z^3,";
  PrintLn "    x^3y^2 - w^2, xz^4 - y^2w^3, yz^5 - xw^5, y^3w^5 - z^7,";
  PrintLn "    x^2w^7 - z^8, z^9 - yw^8)***;";
  PrintLn "    TgCone(I);";
  PrintLn "    Q  := ideal(x, y, z);";
  PrintLn "    PS := PrimaryPoincare(I, Q); PS;";
  PrintLn "    $hp.PSerToHilbert(PS);";
  PrintLn "";
  PrintLn "    Use S ::= ZZ/(32003)[x,y,z];";
  PrintLn "    I  := ideal(0);";
  PrintLn "    Q  := ideal(x, y, z^2);";
  PrintLn "    PS := PrimaryPoincare(I, Q); PS;";
  PrintLn "    $hp.PSerToHilbert(PS);";
  PrintLn "    HV := $hp.PSerHVector(PS); -- the H-vector associated to PS";
  PrintLn "";
  PrintLn "    Use S ::= ZZ/(32003)[x,y,z,w];";
  PrintLn "    I  := ***ideal(-yz + xw, z^3 - yw^2, -xz^2 + y^2w, -y^3 + x^2z)***;";
  PrintLn "    Q  := ideal(x, y, z^2, w^3);";
  PrintLn "    PS := PrimaryPoincare(I, Q); PS;";
  PrintLn "    $hp.PSerToHilbert(PS);";
  PrintLn "    HV := $hp.PSerHVector(PS);";
  PrintLn "    P.E(0, HV);";
  PrintLn "    [ P.E(J,HV) | J In 0..($hp.PSerDim(PS)-2) ];";
  PrintLn "    [ P.E(J,HV) | J In 0..(len(HV)-1) ];";
  PrintLn "";
  PrintLn "    Use S ::= ZZ/(32003)[x,y,z,w];";
  PrintLn "    I  := ***ideal(x^3-y^7, x^2y - xw^3-z^6)***;";
  PrintLn "    Q  := ideal(x, y, z, w);";
  PrintLn "    PS := PrimaryPoincare(I, Q); PS;";
  PrintLn "    $hp.PSerToHilbert(PS);";
  PrintLn "    HV := $hp.PSerHVector(PS);";
  PrintLn "    [ P.E(J,HV) | J In 0..($hp.PSerDim(PS)-2) ];";
  PrintLn "    [ P.E(J,HV) | J In 0..(len(HV)-1) ];";
  PrintLn "";
  PrintLn "    Use S ::= ZZ/(32003)[x,y,z];";
  PrintLn "    I  := ideal(z^3);";
  PrintLn "    Q  := ***ideal(x^2, y^2, xz, yz)***;";
  PrintLn "    PS := PrimaryPoincare(I, Q); PS;";
  PrintLn "    $hp.PSerToHilbert(PS);";
  PrintLn "    HV := $hp.PSerHVector(PS);";
  PrintLn "    [ P.E(J,HV) | J In 0..($hp.PSerDim(PS)-2) ];";
  PrintLn "    [ P.E(J,HV) | J In 0..(len(HV)-1) ];";
EndDefine;

------[   Main functions   ]--------

Define InitialIdeal(I, IndetList)
  P := RingOf(I);
  N := NumIndets(P);
  If GradingDim(P)>1 Or GradingDim(P)=0
    Then error("Not implemented for GradingDim = ", sprint(GradingDim(P)));
  EndIf;
  R := RingOf(WeightsMatrix(P));
  M := ElimMat(ConcatHor(WeightsMatrix(P), RowMat(R,[1])), [N+1]);
  NewP := NewPolyRing(CoeffRing(P), SymbolRange("x",1,N+1), M, 0);
  phiP_NewP := PolyAlgebraHom(P, NewP, first(indets(NewP),N));
  phiNewP_P := PolyAlgebraHom(NewP, P, concat(indets(P),[0]));
  t := (indets(NewP))[N+1];
  IL := apply(phiP_NewP, IndetList);
  HGens := [$.RelHomogenized(F,IL,t) | F In apply(phiP_NewP, gens(I))];
  InGB := interreduced([ InitialForm(F,t) | F In GBasis(ideal(HGens)) ]);
  Return ideal(apply(phiNewP_P, InGB));
EndDefine;


Define TgCone(I)
  P := RingOf(I);
  N := NumIndets(P);
  If GradingDim(P)>1
    Then error("Not implemented for GradingDim = ", sprint(GradingDim(P)));
  EndIf;
  R := RingOf(WeightsMatrix(P));
  If GradingDim(P)>0 Then
    M := HomogElimMat(ConcatHor(WeightsMatrix(P), RowMat(R,[1])), [N+1]);
  Else
    M := ElimMat(N+1, [N+1]);
  EndIf;
  NewP := NewPolyRing(CoeffRing(P), SymbolRange("x",1,N+1), M, GradingDim(P));
  phiP_NewP := PolyAlgebraHom(P, NewP, first(indets(NewP), NumIndets(P)));
  phiNewP_P := PolyAlgebraHom(NewP, P, concat(indets(P),[0]));
  t := (indets(NewP))[N+1];
  HGens := [ homog(F,t) | F In apply(phiP_NewP, gens(I)) ];
  InGB  := interreduced([ InitialForm(F,t) | F In GBasis(ideal(HGens)) ]);
  Return ideal(apply(phiNewP_P, InGB));
EndDefine;


Define PrimaryPoincare(I, Q)
  R := RingOf(I); -- check if = RingOf(Q)
  K := CoeffRing(R);
  NInds := NumIndets(R);
  LGensQ := len(gens(Q));
  W2 := BlockMat2x2(NewMat(1,NInds, 0), NewMat(1,LGensQ, 1),
	      	    NewMat(1,NInds, 1), NewMat(1,LGensQ, 0));
  W1 := ConcatVer(NewMat(1, NInds+LGensQ, 1), W2);
  xt := concat(SymbolRange("x",1,NInds), SymbolRange("t",1,LGensQ));
  R1 := NewPolyRing(K, xt, CompleteToOrd(W1), 1);
  R2 := NewPolyRing(K, xt, CompleteToOrd(W2), 2);
  X := first(indets(R1), NInds);
  T := last(indets(R1), LGensQ);
  phi := PolyAlgebraHom(R, R1, X);
  J := ideal(R1, apply(phi, gens(I))) + ideal(R1, T - apply(phi, gens(Q)));
  IIJ  := InitialIdeal(InitialIdeal(J, T), X);
  IIJinR2 := ideal(apply(PolyAlgebraHom(R1,R2,indets(R2)), gens(IIJ)));
  PS := HilbertSeries(R2/IIJinR2);
  Qt := RingQQt(1);
  Qty := RingQQt(2);
  N  := $hp.HSNumToPoly(Qty, PS)/((1-indet(Qty,2))^NInds);
  Return $hp.PSeries(PolyAlgebraHom(Qty,Qt,[indet(Qt,1),one(Qt)])(N), LGensQ);
EndDefine;


Define E(J, HV)
  Return sum([ binomial(K,J)*HV[K+1] | K In J..(len(HV)-1) ]);
EndDefine;


Define IntersectionMultiplicity(PrPt, F, G)
// need more checks!!!
  If PrPt = [0,0,0] Then error("Not a projective point: "+sprint(PrPt)); EndIf;
  QQxy := NewPolyRing(RingQ(), ["x","y"]);
  x := indets(QQxy);
  If PrPt[1]<>0 Then
    P := (1/PrPt[1])*[PrPt[2],PrPt[3]];
    phi := PolyAlgebraHom(RingOf(F), QQxy, [1, x[1]+P[1], x[2]+P[2]]);
    Return IntMultAff00(phi(F), phi(G));
  EndIf;
  If PrPt[2]<>0 Then
    P := (1/PrPt[2])*[PrPt[1],PrPt[3]];
    phi := PolyAlgebraHom(RingOf(F), QQxy, [x[1]+P[1], 1, x[2]+P[2]]);
    Return IntMultAff00(phi(F), phi(G));
  EndIf;
  If PrPt[3]<>0 Then
    P := (1/PrPt[3])*[PrPt[1],PrPt[2]];
    phi := PolyAlgebraHom(RingOf(F), QQxy, [x[1]+P[1], x[2]+P[2], 1]);
    Return IntMultAff00(phi(F), phi(G));
  EndIf;
EndDefine; -- IntersectionMultiplicity


------[   Auxiliary functions   ]--------

-- Eval0 sostituisce 0 alla indeterminata X
Define Eval0(F, X)  Return NR(F, [X]); EndDefine;

Define IntMultAff00(F,G)
  x := indets(RingOf(F));
  F0 := Eval0(F, x[2]);
  G0 := Eval0(G, x[2]);
  If F0=0 And G0=0 Then
    Return "infty";
  Elif Eval0(F0, x[1])<>0 Or Eval0(G0, x[1])<>0 Then
    Return 0;
  Elif G0=0 Then -- G divisible by x[2]
    Return 1 + IntMultAff00(x[2], F0/x[1]) + IntMultAff00(G/x[2], F);
  Elif F0=0 Then -- F divisible by x[2]
    Return 1 + IntMultAff00(x[2], G0/x[1]) + IntMultAff00(F/x[2], G);
  Elif deg(G0) <= deg(F0) Then
    Return IntMultAff00(G, F - G*(F0-NR(F0,[G0]))/G0);
  Else
    Return IntMultAff00(F, G - F*(G0-NR(G0,[F0]))/F0);
  EndIf;
EndDefine; -- IntMultAff00


Define InitialForm(F, t)
  D := deg(LT(F), t);
  Return sum([M In monomials(F) | deg(M,t)=D])/t^D;
EndDefine;


Define RelHomogenized(F, L, t)
  If F=0 Then Return F; EndIf;
  Ms  := monomials(F);
  IIL := [IndetIndex(X) | X In L];
  RelDeg := Func(M) Lg:=log(M); Return sum([Lg[X]|X In IIL]); EndFunc;
  RelDegs := [ RelDeg(M) | M In monomials(F)];
  D := max(RelDegs);
  Return sum([ Ms[I]*t^(D-RelDegs[I]) | I In 1..len(Ms)]);
EndDefine;


Define ElimOrdMat()
EndDefine; -- ElimOrdMat


EndPackage; -- Package $primary
