Package $hp -- Hilbert-Poincare Series

export AffHilbertFn;
export AffHilbertSeries;
export EvalHilbertFn;
export HVector;
export HilbertFn;
export HilbertPoly;
export HilbertSeries;
export RegularityIndex;
export dim;
export multiplicity;
//export HilbertSeriesMultiDeg;

Define About()
  PrintLn "    Author:  A.M.Bigatti";
  PrintLn "    Date:      17 Nov 2005";
  PrintLn "    Added:   non-homogeneous cases";
  PrintLn "    Date:      26 May 2003";
  PrintLn "    Fixed:   PoincareShifts(...)";
  PrintLn "    Date:      13 May 2003";
  PrintLn "    Added:   EvalHilbertFn for univariate non-standard PS";
  PrintLn "    Date:      2 Aug 2002";
  PrintLn "    Added:   Poincare for modules with shifts";
  PrintLn "    Date:      7 Feb 2000";
  PrintLn "    BugFix:  ZMOD Hilbert Function values";
  PrintLn "    BugFix:  TruncBin for BinExp (fixed 0-diml HF accordingly)";
  PrintLn "    BugFix:  DenIsStandard used SP.PPExp";
  PrintLn "    BugFix:  Zero dimensional EvalHilbertFn";
  PrintLn "    Date:      10 Apr 2001";
  PrintLn "    BugFix:  New Tags handling";
EndDefine; -- About

--------------------------------------------------------------

-- TAG "PSeries"; --  [P]oincare [Series]
--     "HSDen";   --  [H]ilbert [S]eries [Den]ominator
--     "Hilbert"; --  [Hilbert] Function
-- TYPE PSeries := [SpPoly,HSDen];      --> TAGGED "PSeries"
--      Hilbert := [INT,..,INT,SpPoly]  --> TAGGED "Hilbert"
--      HSDen   := LIST(SpPP);          --> TAGGED "HSDen"


Define AffHilbertSeries(M)
  If type(M) = RING Then
    If IsPolyRing(M) Then W := WeightsMatrix(M);
    Elif IsQuotientRing(M) Then  W := WeightsMatrix(BaseRing(M));
    Else error("PolyRing or QuotientRing expected");
    EndIf;
  Else W := WeightsMatrix(RingOf(M));
  EndIf;
  If NumRows(W)>1 Then
    error("AffHilbertSeries: grading must be a positive ZZ-grading: "+Sprint(W));
  EndIf;
  If MakeSet(GetRow(W,1)) <> [1] Then
    error("AffHilbertSeries: first row of ordering matrix must have only 1");
  EndIf;
  If type(M) = RING Then
    If IsQuotientRing(M) Then PS := PoincareQuotient(LT(DefiningIdeal(M)));
    Elif IsPolyRing(M) Then PS := PoincareQuotient(ideal(M,[]));
    EndIf;
  Else
    PS := HilbertSeries(LT(M));
  EndIf;
  NewDenFactors := untagged(HSDen(PS));
  append(ref NewDenFactors.multiplicities, 1);
  append(ref NewDenFactors.factors, 1-indet(RingQQt(1),1));
  Return $.Make(untagged(PS).num, NewDenFactors);
EndDefine; -- Poincare


Define IsStdGraded(R) -- ring
  If GradingDim(R) <> 1 Then Return False; EndIf;
  Return WeightsMatrix(R) = NewMat(1, NumIndets(R),  1);
EndDefine; -- IsStdGraded

  
Define PoincareQuotient(I)
  Return tagged(HilbertSeriesQuot(I), "PSeries");
EndDefine; -- PoincareQuotient


Define HilbertSeries(X, Opt ShiftsList)
  ERRShifts := "ShiftsList only for MODULE";
  If type(X) = RING Then
    If IsPolyRing(X) Then Return PoincareQuotient(ideal(X,[]));
    Elif IsQuotientRing(X) Then
      If IsDefined(ShiftsList) Then error(ERRShifts); EndIf;
      Return PoincareQuotient(DefiningIdeal(X));
    EndIf;
  Elif type(X) = IDEAL Then
    If IsDefined(ShiftsList) Then error(ERRShifts); EndIf;
    R := RingOf(X);
    F0 := HilbertNumQuot(ideal(R,[]));
    FX := HilbertNumQuot(X);
    If RingOf(F0) = RingOf(FX) Then
      F := F0-FX;
    Else
      F := F0 - apply(PolyAlgebraHom(RingOf(FX), RingOf(F0), indets(RingOf(F0))), FX);
    EndIf;
    Return PSeries(F, DenMake(R));
  Elif type(X) = MODULE Then
    If not(IsDefined(ShiftsList)) Then ShiftsList := []; EndIf;
    Return $.PoincareShifts(X, ShiftsList);
  EndIf;
  error("RING, IDEAL, or MODULE expected");
EndDefine; -- poincare


Define PoincareShifts(M, ShiftsList)  // M: MODULE
  If type(M) <> MODULE Then error("MODULE expected"); EndIf;
  If not(IsHomog(M)) Then
    PrintLn "WARNING! HilbertPoincare input not homogeneous: computing LT...";
  EndIf;
  NC := NumCompts(M);
  If ShiftsList = [] Then ShiftsList := NewList(NC,0); EndIf;
  If len(ShiftsList)<>NC Then
    error("Poincare: Shifts List has wrong length");
  EndIf;
  PS0 := $.QuotUPoincareShifts(SubmoduleRows(ModuleOf(M),ZeroMat(BaseRing(M),0,NC)), ShiftsList);
  PS := $.QuotUPoincareShifts(LT(M), ShiftsList);
  PN := untagged(PS0).num - untagged(PS).num;
  PSM := $.PSeries(PN, untagged(PS).DenFactors);
  Return PSM;
EndDefine; -- UPoincareShifts


Define QuotUPoincareShifts(M, ShiftsList)
  If type(M)=IDEAL Then NC := 1; Else NC := NumCompts(M); EndIf;
  If len(ShiftsList)<>NC Then
    Return error("Poincare: Shifts List has wrong length");
  EndIf;
  If [ N In ShiftsList | N<0 ] <> [] Then
    Return error("Poincare: Shifts must be non-negative");
  EndIf;
  If type(M) = IDEAL Then
    LT_Ideals := [ LT(M) ];
  Else
    LT_Ideals := [ ideal(RingOf(M), [V[I] | V In gens(LT(M))]) | I In 1..NC ];
  EndIf;
  P := RingOf(LT_Ideals[1]);
  PSs  := [ PoincareQuotient(LT_Ideals[I]) | I In 1..NC ];
  Den := untagged(HSDen(first(PSs)));
  QQt := RingQQt(1);
  t := indet(QQt, 1);
  PN := sum([t^(ShiftsList[I])*untagged(PSs[I]).num | I In 1..NC]);
  Return PSeries(PN, Den);
EndDefine; -- QuotUPoincareShifts


Define HilbertFn(...)
  If len(ARGV)=1 Then
    Return $.PSerToHilbert(HilbertSeries(ARGV[1]));
  Else
    Return $.PSerToHilbert(HilbertSeries(ARGV[1]),ARGV[2]);
  EndIf;
EndDefine; -- Hilbert

  
Define AffHilbertFn(...)
  If len(ARGV)=1 Then
    Return $.PSerToHilbert(AffHilbertSeries(ARGV[1]));
  Else
    Return $.PSerToHilbert(AffHilbertSeries(ARGV[1]),ARGV[2]);
  EndIf;
EndDefine; -- Hilbert

	
Define HilbertPoly(Q)
  QQt := RingQQt(1);
  Return $.PSerToHilbertPoly(HilbertSeries(Q));
EndDefine;


Define HVector(Q)  Return $.PSerHVector(HilbertSeries(Q)); EndDefine;

Define multiplicity(M)
  If type(M) = RING Then
    If IsPolyRing(M) Then Return 1;
    Elif IsQuotientRing(M) Then
      P := BaseRing(M);
      If not(IsStdGraded(P)) Then
	R := NewPolyRing(CoeffRing(P), SymbolRange("dim_indet",1,NumIndets(P)));
	phi := PolyAlgebraHom(P, R, indets(R));
	Return multiplicity(R/ideal(apply(phi, gens(DefiningIdeal(M)))));
      EndIf;
      If IsHomog(DefiningIdeal(M)) Then Return MultHom(M); Else Return MultAff(M); EndIf;
    Else MultErr := "Not defined for this type of ring"; error(MultErr);
    EndIf;
  EndIf;
  If IsHomog(M) Then Return MultHom(M); Else Return MultAff(M); EndIf;
EndDefine;


Define MultHom(M) Return $.PSerMultiplicity(HilbertSeries(M)); EndDefine;
Define MultAff(M) Return $.PSerMultiplicity(AffHilbertSeries(M)); EndDefine;

Define dim(M)
  If type(M) = RING Then
    If IsPolyRing(M) Then Return NumIndets(M);
    Elif IsQuotientRing(M) Then
      P := BaseRing(M);
      I := DefiningIdeal(M);
      If not(IsStdGraded(P)) Then
	R := NewPolyRing(CoeffRing(P), SymbolRange("dim_indet",1,NumIndets(P)));
	phi := PolyAlgebraHom(P, R, indets(R));
	Return dim(R/ideal(apply(phi, gens(I))));
      EndIf;
      If IsHomog(I) Then Return DimHom(M); Else Return DimAff(M);
      EndIf;
    Else DimErr := "Not defined for this type of ring"; error(DimErr);
    EndIf;
  EndIf;
  If IsHomog(M) Then Return DimHom(M); Else Return DimAff(M); EndIf;
EndDefine;

Define DimHom(M) Return $.PSerDim(HilbertSeries(M)); EndDefine;
Define DimAff(M) Return $.PSerDim(AffHilbertSeries(M))-1; EndDefine;

Define RegularityIndex(X)
  If tag(X) = "PSeries" Then Return $.PSerRegularityIndex(X); EndIf;
  If tag(X) = "Hilbert" Then Return $.HFRegularityIndex(X); EndIf;
  RegularityIndexERR := Sprint(type(tagged(0,"PSeries"))) +
	" or " +  sprint(type(tagged(0,"Hilbert"))) + " expected";
  error(RegularityIndexERR);
EndDefine; -- RegularityIndex

-------------------------------

Define PSerHVector(PSer)
  QQt := RingQQt(1);
  Return ToHVec(untagged(HSSimplified(PSer)).num);
EndDefine;

Define PSerDim(PSer)
  Return sum(untagged(HSDen(HSSimplified(PSer))).multiplicities);
EndDefine;

Define PSerMultiplicity(PSer)
  Return AsINT(sum($.PSerHVector(PSer)));
EndDefine;

-------------------------------

Define HSDen(PSer)
  Return tagged(untagged(PSer).DenFactors,"HSDen");
EndDefine; -- Den

Define Make(N, HSDen)
  P := RingOf(N);
  If P<>RingQQt(NumIndets(P)) Then
    QQt := RingQQt(NumIndets(P));
    N := PolyAlgebraHom(P, QQt, indets(QQt))(N);
  EndIf;
  Return tagged(Record[num:=N, DenFactors:=HSDen], "PSeries");
EndDefine; -- Make

Define IsStandard(PSer)
  Return $.DenIsStandard(untagged(PSer).DenFactors);
EndDefine; -- IsStandard

Define DenMakeStandard(N)
//  Return NewList(N,[1]);
  QQt := RingQQt(1);
  Return Record[multiplicities := [N],
		factors := [1-indet(QQt,1)],
		RemainingFactor := one(QQt)];
EndDefine; -- DenMakeStandard

Define DenMake(P)
  If IsStdGraded(P) Then Return DenMakeStandard(NumIndets(P)); EndIf;
//  Return [ wdeg(x) | x In indets(P)];
  QQt := RingQQt(GradingDim(P));
  Return Record[multiplicities := NewList(NumIndets(P),1),
		factors := [1 - LogToTerm(QQt, wdeg(x)) | x In indets(P)],
		RemainingFactor := one(QQt)];
EndDefine; -- DenMake

Define HSDenToPoly(Kx, PSer)
//  Return product([1-LogToTerm(Kx, PP) | PP In untagged(HSDen(PSer))]);
  Factors := untagged(HSDen(PSer)).factors;
  Exp := untagged(HSDen(PSer)).multiplicities;
  QQt := RingOf(Factors[1]);
  phi := PolyAlgebraHom(QQt, Kx, indets(Kx));
  Return product([(phi(Factors[i]))^Exp[i] | i In 1..len(Factors)]);
EndDefine; -- HSDenToPoly

Define HSNumToPoly(Kx, PSer)
  QQt := RingOf(untagged(PSer).num);
  Return PolyAlgebraHom(QQt, Kx, indets(Kx))(untagged(PSer).num);
EndDefine; -- HSNumToPoly

Define ToRatFun(FrFld, PSer)
  N := HSNumToPoly(BaseRing(FrFld), PSer);
  D := HSDenToPoly(BaseRing(FrFld),PSer);
  phi := EmbeddingHom(FrFld);
  Return phi(N)/phi(D);
EndDefine; -- ToRatFun

-------------------------------

Define PSeries(P, DenExp_Or_Den)
  If type(DenExp_Or_Den)=INT Then
    DenExp := DenExp_Or_Den;
    If DenExp<0 Then error("PSeries: Expected non-negative INT"); EndIf;
    Den := $.DenMakeStandard(DenExp);
  Else
    Den := DenExp_Or_Den;
  EndIf;
  Return $.Make(P, Den);
EndDefine;

-------------------------------

Define HSSimplified(PSer)
  QQt := RingQQt(1);
  If not($.IsStandard(PSer)) Then
    error("HSSimplified: Operator not available for non-standard HPSeries");
  EndIf;
  If IsZero(untagged(PSer).num) Then Return
    $.Make(zero(QQt), $.DenMakeStandard(0));
  EndIf;
  t := indet(QQt, 1);
  HPNum    := untagged(PSer).num;
//  HPDenExp := len(untagged(HSDen(PSer)));
  HPDenExp := sum(untagged(PSer).DenFactors.multiplicities);
  While IsDivisible(HPNum, 1-t) Do
    HPNum    := HPNum/(1-t);
    HPDenExp := HPDenExp - 1;
  EndWhile;
  Return $.Make(HPNum, $.DenMakeStandard(HPDenExp));
EndDefine;

-------------------------------

Define DenIsStandard(HSDen)
  fact := 1-indet(RingOf(HSDen.RemainingFactor),1); // 1-t
  Foreach f In HSDen.factors Do
    If f<>fact Then Return false; EndIf;
  EndForeach;
  Return True;
EndDefine; -- DenIsStandard

--------------------------------
--<<   Hilbert Functions    >>--
--------------------------------

Define PSerToHilbert(...)
  TopLevel ERR;
  If len(ARGV) = 2 Then
    Return $.EvalHilbertFn(ARGV[1],ARGV[2]);
  Elif
    len(ARGV) = 1 Then
    Return $.PSerToHilbertFn(ARGV[1]);
  Else error(ERR.BAD_PARAMS_NUM, ": 1 or 2 expected (Hilbert)");
  EndIf;
EndDefine;

-------------------------------

Define AuxHilbertPoly(HSSimplifiedPS)
  SPS := untagged(HSSimplifiedPS);
  If IsZero(SPS.num) Then Return 0; EndIf;
  HV  := ToHVec(SPS.num);
  DIM := sum(SPS.DenFactors.multiplicities);
  Return sum([HV[I]*Bin(indet(RingQQt(1),1)+DIM-I,DIM-1) | I In 1..len(HV)]);
EndDefine; -- AuxHilbertPoly


Define PSerToHilbertPoly(PS)
  QQt := RingQQt(1);
  SPS := $.HSSimplified(PS);
  If untagged(HSDen(SPS)).factors = [] Then Return 0; EndIf;
  Return $.AuxHilbertPoly(SPS);
EndDefine; -- PSerToHilbertPoly


Define ToHVec(SpP)
  If IsZero(SpP) Then Return []; EndIf;
  t := indet(RingOf(SpP),1);
  Skeleton := [t^i | i In 0..deg(SpP)];
  HV := coefficients(SpP, Skeleton);
  Return [AsINT(n) | n In HV];
EndDefine;

Define PSerToHilbertFn(PS)
  QQt := RingQQt(1);
  If IsZero(untagged(PS).num) Then Return tagged([ [], 0],"Hilbert"); EndIf;
  SPS := untagged(HSSimplified(PS));
//  DIM := len(untagged(HSDen(SPS)));
  DIM := sum(SPS.DenFactors.multiplicities);
  HV := ToHVec(SPS.num);
  If DIM=0 Then Return tagged([HV,0],"Hilbert"); EndIf;
  REG := len(HV)-DIM;
  If REG > 0 Then
    TrBins := concat( NewList(REG-1,0),
                   [$.TruncBin(I,DIM-1) | I In (DIM-1)..(DIM+REG)]);
    HF1 :=[sum([HV[I]*TrBins[N-I+REG+1] | I In 1..(N+1)]) | N In 0..(REG-1)];
  Else
    HF1 := [];
  EndIf;
  Return tagged([HF1, PSerToHilbertPoly(SPS)],"Hilbert");
EndDefine; -- PSerToHilbertFn


Define PSerRegularityIndex(PS)
  QQt := RingQQt(1);
  If IsZero(untagged(PS).num) Then Return 0; EndIf;
  SPS := untagged(HSSimplified(PS));
  DIM := sum(SPS.DenFactors.multiplicities);
  HV := ToHVec(SPS.num);
  Return len(HV)-DIM;
EndDefine; -- PSerRegularityIndex


Define EvalHilbertFn(X, N)
  TopLevel ERR;
  TopLevel CurrentRing;
  QQt := RingQQt(1);
  t := RingElem(QQt, "t");
  If type(X)=type(tagged([],"Hilbert")) Then
    If N<len(untagged(X)[1]) Then Return untagged(X)[1,N+1];
    Else  Return LC(subst(untagged(X)[2],[[t, N]]));
    EndIf;
  EndIf;
  If type(X)=type(tagged([],"PSeries")) Then
    If IsZero(untagged(X).num) Then Return 0; EndIf;
    If $.IsStandard(X) Then
      SPS := untagged(HSSimplified(X));
      HV  := ToHVec(SPS.num);
      DIM := sum(SPS.DenFactors.multiplicities);
      If DIM = 0 Then
	If 0<=N And N<len(HV) Then Return HV[N+1]; Else Return 0; EndIf;
      EndIf;
      Return sum([HV[I]*$.TruncBin(N+DIM-I,DIM-1) | I In 1..len(HV)]);
    EndIf;
    // univariate and not standard
    HPDen := den(X);
//    If len(HPDen[1]) = 1 Then
    If NumIndets(RingOf(HPDen.factors[1])) = 1 Then
//      Using QQt Do
	D := 1;
-- 	Foreach W In HPDen Do
-- 	  D := D * sum([t^(I*W[1]) | I In 0..div(N, W[1])]);
-- 	  D := NR(D, [t^(N+1)]);
-- 	EndForeach;
	For i := 1 To len(HPDen.factors) Do
	  W := log(LT(HPDen.factors[i]));
	  For j:=1 To HPDen.multiplicities[i] Do
	    D := D * sum([t^(I*W[1]) | I In 0..div(N, W[1])]);
	    D := NR(D, [t^(N+1)]);
	  EndFor;
	EndFor;
	F := NR(D * HSNumToPoly(X), [t^(N+1)]);
	RatResult := CoeffOfTerm(t^N, F);
	Return RatResult.Num;
//      EndUsing;
    EndIf;
  EndIf;
  error(ERR.BAD_PARAMS, "(EvalHilbertFn)");
EndDefine; -- EvalHilbertFn


Define HFRegularityIndex(HF)
  HFSmall := untagged(HF)[1];
  HPoly := untagged(HF)[2];
  For D := len(HFSmall)-1 To 0 Step -1 Do
    If HFSmall[D+1] <> eval(HPoly, [D]) Then Return D+1; EndIf;
  EndFor;
  If HPoly = 0 Then error("cannot compute RegularityIndex for 0 function"); EndIf;
  DD := -1;
  While True Do
    If eval(HPoly, [DD])<>0 Then Return DD+1; EndIf;
    DD := DD-1;
  EndWhile;
  Return 0;
EndDefine; -- HFRegularityIndex


------[   pretty printing   ]--------

Define Tagged(X,T)
  Return tagged(X, T);
EndDefine;


Define Print_HSDen(D, HSDen)
  Factors := untagged(HSDen).factors;
  Exp := untagged(HSDen).multiplicities;
  HSRing := RingOf(untagged(HSDen).RemainingFactor);
  L := len(Factors);
  If $.DenIsStandard(HSDen) Then
    Print "(1-",indet(HSRing, 1),")" On D;
    If sum(Exp)<>1 Then  Print "^", sum(Exp) On D; EndIf
  Else
    If L<>1 Then  Print "( " On D;  EndIf;
    Print "(1-", LT(Factors[1]), ")" On D;
    If Exp[1]<>1 Then Print "^", Exp[1] On D; EndIf;
    For i := 2 To len(Factors) Do
      Print "*(1-", LT(Factors[i]), ")" On D;
      If Exp[i]<>1 Then Print "^", Exp[i] On D; EndIf;
    EndFor;
    If L<>1 Then  Print " )" On D;  EndIf;
  EndIf;
EndDefine;


Define Print_PSeries(D, PSer)
  QQt := RingQQt(1);
  If $.IsStandard(PSer) Then
    PSer := $.HSSimplified(PSer);
  Else
    PrintLn "---  Non-simplified HilbertPoincare' Series  ---" On D;
  EndIf;
  SpPolyNum := $sppoly.PolyToSPPoly(untagged(PSer).num);
  SpPolyNum := $sppoly.Tagged(reversed(untagged(SpPolyNum)),"SpPoly");
  Print "(", SpPolyNum, ")" On D;
  HP_Den := untagged(PSer).DenFactors;
  If len(untagged(HP_Den).factors)>0 Then Print " / ", tagged(HP_Den, "HSDen") On D;  EndIf;
EndDefine; -- Print_PSeries

Define Print_Hilbert(D, HF)
  For I:=1 To len(HF[1]) Do PrintLn "H(",I-1,") = ",HF[1,I] On D; EndFor;
  Print "H(t) = ", HF[2], "   for t >= ", len(HF[1]) On D;
EndDefine;


Define TruncBin(A,B)
  If A<B Or B<0 Then Return 0; EndIf;
  Return Bin(A,B);
EndDefine; -- TruncBin

------[ Hilbert multideg ]-----------------------------------------------

Define HilbertSeriesMultiDeg(QR, WM)
  Return $hp.PoincareMultiDeg(QR, WM);
EndDefine;

Define PoincareMultiDeg(QR, WM)
  TopLevel CurrentRing;
  P := CurrentRing;
  If type(QR) <> RING And not(IsQuotientRing(QR)) Then
    error("PoincareMultiDeg: first argument must be a quotient ring", QR);
  EndIf;
  If type(WM) <> MAT Then
    error("PoincareMultiDeg: second argument must be a matrix", WM);
  EndIf;
  If NumCols(WM) <> NumIndets(P) Then
    error("PoincareMultiDeg: wrong number of weights:" + Sprint(NumCols(WM)) + " " + Sprint(NumIndets(P)), WM);
  EndIf;
  ZeroesInFirstRow := [ W In GetRow(WM,1) | W<=0 ] <> [];
  If ZeroesInFirstRow Then WM := PositiveGrading4(WM); EndIf;
  K := CoeffRing(P);
  AuxRing ::= K[x[1..len(WM[1])]]; //, Weights(WM);
  PS := HilbertSeries(AuxRing/($contrib/specvar.CanonicalImage(DefiningIdeal(QR))));
  If not(ZeroesInFirstRow) Then
    Return PS;
  EndIf;
  PSNum := [ [ T[1], tail(T[2]) ] | T In untagged(PS).num];
  HSDen := [ tail(T) | T In untagged(PS).DenFactors];
  Return $.Make(PSNum, HSDen);
EndDefine; -- PoincareOfSeparated



-- EXAMPLES
-- Use R ::= QQ[x,y,z];
-- //PoincareMultideg(ideal(Indets())^2, LexMat(NumIndets()));

-- WM := mat([[1,0,0],[1,-1,0]]);
-- PoincareMultiDeg(R/ideal(Indets())^2, WM);

-- WM := mat([[1,0,0],[-1,1,1]]);
-- PoincareMultiDeg(R/ideal(0), WM);

-- WM := mat([[1,7,0],[0,-5,1]]);
-- PositiveGrading(WM);

------[ end of Hilbert multideg ]-------------------------------------------

----------------------------------------------------------------------
PrintTagged := Record[
		      PSeries := $.Print_PSeries,
		      HSDen := $.Print_HSDen,
		      Hilbert := $.Print_Hilbert
		      ];

EndPackage; -- Package $hp

