1Package $hp -- Hilbert-Poincare Series
2
3export AffHilbertFn;
4export AffHilbertSeries;
5export EvalHilbertFn;
6export HVectorPkg;
7export HilbertFn;
8export HilbertPolyPkg;
9export HilbertSeries;
10export HilbertSeriesMultiDeg;
11export HilbertSeriesShifts;
12export RegularityIndex;
13export dim;
14export multiplicity;
15export TaggedHilbertFn;
16export TaggedHilbertSeries;
17//export HilbertSeriesMultiDeg;
18
19Define About()
20  PrintLn "    Author:  A.M.Bigatti";
21EndDefine; -- About
22
23--------------------------------------------------------------
24
25-- TAG
26--  TagForPSeries();
27--  TagForHilbertFn();
28--     "HSDen";   --  [H]ilbert [S]eries [Den]ominator
29-- TYPE PSeries := [SpPoly,HSDen];      --> TAGGED "PSeries"
30--      Hilbert := [INT,..,INT,SpPoly]  --> TAGGED "HilbertFn"
31--      HSDen   := LIST(SpPP);          --> TAGGED "HSDen"
32
33
34Define AffHilbertSeries(M)
35  If type(M) = RING Then
36    If   IsPolyRing(M) Then return AffHilbertSeries(M/ideal(M,[]));
37    Elif IsQuotientRing(M) Then P := BaseRing(M);
38    Else error("PolyRing or QuotientRing expected");
39    EndIf;
40  Else P := RingOf(M);
41  EndIf;
42--   If GradingDim(P)>1 Then
43--     error("Grading must be a positive ZZ-grading: "+Sprint(W));
44--   EndIf;
45  If not(IsStdGraded(P)) Then error("Ring must be standard graded"); EndIf;
46  If type(M) = RING Then --> IsQuotientRing(M)
47    PS := PoincareQuotient(LT(DefiningIdeal(M)));
48  Else
49    PS := HilbertSeries(LT(M));
50  EndIf;
51  NewDenFactors := untagged(HSDen(PS));
52  append(ref NewDenFactors.multiplicities, 1);
53  append(ref NewDenFactors.factors, 1-indet(RingQQt(1),1));
54  Return $.Make(untagged(PS).num, NewDenFactors);
55EndDefine; -- AffHilbertSeries
56
57
58Define PoincareQuotient(I)
59  return TaggedHilbertSeries(HilbertSeriesQuot(I));
60EndDefine; -- PoincareQuotient
61
62
63Define HilbertSeries(X)
64  If type(X) = RING Then
65    If IsPolyRing(X) Then return PoincareQuotient(ideal(X,[]));
66    Elif IsQuotientRing(X) Then
67      return PoincareQuotient(DefiningIdeal(X));
68    EndIf;
69  Elif type(X) = IDEAL Then
70    R := RingOf(X);
71    F0 := HilbertNumQuot(ideal(R,[]));
72    FX := HilbertNumQuot(X);
73    If RingOf(F0) = RingOf(FX) Then
74      F := F0-FX;
75    Else
76      F := F0 - apply(PolyAlgebraHom(RingOf(FX), RingOf(F0), indets(RingOf(F0))), FX);
77    EndIf;
78    return PSeries(F, DenMake(R));
79  Elif type(X) = MODULE Then
80    return HilbertSeriesModule(X);
81  EndIf;
82  error("RING, IDEAL, or MODULE expected");
83EndDefine; -- HilbertSeries
84
85
86Define HilbertSeriesShifts(M, ShiftsList)  // M: MODULE
87  If type(M) <> MODULE Then error("MODULE expected"); EndIf;
88  NC := NumCompts(M);
89  If len(ShiftsList)<>NC Then
90    error("ShiftsList has wrong length");
91  EndIf;
92  F := NewFreeModule(RingOf(M), ColMat(ShiftsList));
93  M := SubmoduleRows(F, GensAsRows(M));
94  return HilbertSeriesModule(M);
95EndDefine; -- HilbertSeriesShifts
96
97
98Define HilbertSeriesModule(M)
99  If not(IsHomog(M)) Then
100    PrintLn "WARNING! HilbertPoincare input not homogeneous: computing LT...";
101  EndIf;
102  P := RingOf(M);
103  F := ModuleOf(M); -- the FreeModule
104  ShiftsList := [ shifts(F)[i,1] | i in 1..NumCompts(F)];
105  LT_Ideals := [ ideal(RingOf(M), row) | row in GetRows(GensAsCols(LT(M))) ];
106  PSs  := [ HilbertSeries(I) | I In LT_Ideals ];
107  Den := untagged(HSDen(first(PSs)));
108  t := indet(RingQQt(1), 1);
109  PN := sum([t^(ShiftsList[i])*untagged(PSs[i]).num | i in 1..len(PSs)]);
110  return PSeries(PN, Den);
111EndDefine; -- HilbertSeriesModule
112
113
114Define HilbertFn(...)
115  If len(ARGV)=1 Then
116    return $.PSerToHilbert(HilbertSeries(ARGV[1]));
117  Else
118    return $.PSerToHilbert(HilbertSeries(ARGV[1]),ARGV[2]);
119  EndIf;
120EndDefine; -- Hilbert
121
122
123Define AffHilbertFn(...)
124  If len(ARGV)=1 Then
125    return $.PSerToHilbert(AffHilbertSeries(ARGV[1]));
126  Else
127    return $.PSerToHilbert(AffHilbertSeries(ARGV[1]),ARGV[2]);
128  EndIf;
129EndDefine; -- Hilbert
130
131
132Define HilbertPolyPkg(Q)
133  QQt := RingQQt(1);
134  return $.PSerToHilbertPoly(HilbertSeries(Q));
135EndDefine;
136
137
138Define HVectorPkg(Q)  return $.PSerHVector(HilbertSeries(Q)); EndDefine;
139
140Define multiplicity(M)
141  If type(M) = RING Then
142    If IsPolyRing(M) Then return 1;
143    Elif IsQuotientRing(M) Then
144      P := BaseRing(M);
145      If not(IsStdGraded(P)) Then
146	R := NewPolyRing(CoeffRing(P), SymbolRange("dim_indet",1,NumIndets(P)));
147	phi := PolyAlgebraHom(P, R, indets(R));
148	return multiplicity(R/ideal(apply(phi, gens(DefiningIdeal(M)))));
149      EndIf;
150      If IsHomog(DefiningIdeal(M)) Then return MultHom(M); Else return MultAff(M); EndIf;
151    Else MultErr := "Not defined for this type of ring"; error(MultErr);
152    EndIf;
153  EndIf;
154  If IsHomog(M) Then return MultHom(M); Else return MultAff(M); EndIf;
155EndDefine;
156
157
158-- Define MultHom(M) return $.PSerMultiplicity(HilbertSeries(M)); EndDefine;
159Define MultHom(M) return MultiplicityQuot(DefiningIdeal(M)); EndDefine;
160Define MultAff(M) return $.PSerMultiplicity(AffHilbertSeries(M)); EndDefine;
161
162Define dim(M)
163  If type(M) = RING Then
164    If IsPolyRing(M) Then return NumIndets(M); EndIf;
165    if not(IsQuotientRing(M)) then error("Not defined for this type of RING"); endif;
166    --> QuotientRing
167    P := BaseRing(M);
168    I := DefiningIdeal(M);
169    If not(IsStdGraded(P)) Then
170      If HasGBasis(I) Then I := LT(I); EndIf;
171      R := NewPolyRing(CoeffRing(P), SymbolRange("dim_indet",1,NumIndets(P)));
172      return dim(R/ideal(apply(PolyAlgebraHom(P,R,indets(R)), gens(I))));
173    EndIf;
174    If IsHomog(I) Then return DimHom(M); Else return DimAff(M); EndIf;
175  EndIf;
176  --> module
177  If IsHomog(M) Then return DimHom(M); Else return DimAff(M); EndIf;
178EndDefine;
179
180
181define DimHom(M) return DimQuot(DefiningIdeal(M)); enddefine;
182define DimAff(M) return DimQuot(LT(DefiningIdeal(M))); enddefine;
183
184
185Define RegularityIndex(X)
186  If tag(X) = TagForPSeries() Then return $.PSerRegularityIndex(X); EndIf;
187  If tag(X) = TagForHilbertFn() Then return $.HFRegularityIndex(X); EndIf;
188  RegularityIndexERR := sprint(type(TaggedHilbertSeries(0))) +
189	" or " +  sprint(type(TaggedHilbertFn(0))) + " expected";
190  error(RegularityIndexERR);
191EndDefine; -- RegularityIndex
192
193-------------------------------
194
195Define PSerHVector(PSer)
196  QQt := RingQQt(1);
197  return ToHVec(untagged(HSSimplified(PSer)).num);
198EndDefine;
199
200Define PSerDim(PSer)
201  return sum(untagged(HSDen(HSSimplified(PSer))).multiplicities);
202EndDefine;
203
204Define PSerMultiplicity(PSer)
205  return AsINT(sum($.PSerHVector(PSer)));
206EndDefine;
207
208-------------------------------
209
210Define HSDen(PSer)
211  return tagged(untagged(PSer).DenFactors,"HSDen");
212EndDefine; -- Den
213
214Define Make(N, HSDen)
215  P := RingOf(N);
216  If P<>RingQQt(NumIndets(P)) Then
217    QQt := RingQQt(NumIndets(P));
218    N := PolyAlgebraHom(P, QQt, indets(QQt))(N);
219  EndIf;
220  return TaggedHilbertSeries(record[num:=N, DenFactors:=HSDen]);
221EndDefine; -- Make
222
223Define IsStandard(PSer)
224  return $.DenIsStandard(untagged(PSer).DenFactors);
225EndDefine; -- IsStandard
226
227Define DenMakeStandard(N)
228//  return NewList(N,[1]);
229  QQt := RingQQt(1);
230  return Record[multiplicities := [N],
231		factors := [1-indet(QQt,1)],
232		RemainingFactor := one(QQt)];
233EndDefine; -- DenMakeStandard
234
235Define DenMake(P)
236  If IsStdGraded(P) Then return DenMakeStandard(NumIndets(P)); EndIf;
237//  return [ wdeg(x) | x In indets(P)];
238  QQt := RingQQt(GradingDim(P));
239  return Record[multiplicities := NewList(NumIndets(P),1),
240		factors := [1 - MakeTerm(QQt, wdeg(x)) | x In indets(P)],
241		RemainingFactor := one(QQt)];
242EndDefine; -- DenMake
243
244Define HSDenToPoly(Kx, PSer)
245//  return product([1-MakeTerm(Kx, PP) | PP In untagged(HSDen(PSer))]);
246  Factors := untagged(HSDen(PSer)).factors;
247  Exp := untagged(HSDen(PSer)).multiplicities;
248  QQt := RingOf(Factors[1]);
249  phi := PolyAlgebraHom(QQt, Kx, indets(Kx));
250  return product([(phi(Factors[i]))^Exp[i] | i In 1..len(Factors)]);
251EndDefine; -- HSDenToPoly
252
253Define HSNumToPoly(Kx, PSer)
254  QQt := RingOf(untagged(PSer).num);
255  return PolyAlgebraHom(QQt, Kx, indets(Kx))(untagged(PSer).num);
256EndDefine; -- HSNumToPoly
257
258Define ToRatFun(FrFld, PSer)
259  N := HSNumToPoly(BaseRing(FrFld), PSer);
260  D := HSDenToPoly(BaseRing(FrFld),PSer);
261  phi := EmbeddingHom(FrFld);
262  return phi(N)/phi(D);
263EndDefine; -- ToRatFun
264
265-------------------------------
266
267Define PSeries(P, DenExp_Or_Den)
268  If type(DenExp_Or_Den)=INT Then
269    DenExp := DenExp_Or_Den;
270    If DenExp<0 Then error("PSeries: Expected non-negative INT"); EndIf;
271    Den := $.DenMakeStandard(DenExp);
272  Else
273    Den := DenExp_Or_Den;
274  EndIf;
275  return $.Make(P, Den);
276EndDefine;
277
278-------------------------------
279
280Define HSSimplified(PSer)
281  QQt := RingQQt(1);
282  If not($.IsStandard(PSer)) Then
283    error("HSSimplified: Operator not available for non-standard HPSeries");
284  EndIf;
285  If IsZero(untagged(PSer).num) Then return
286    $.Make(zero(QQt), $.DenMakeStandard(0));
287  EndIf;
288  t := indet(QQt, 1);
289  HPNum    := untagged(PSer).num;
290//  HPDenExp := len(untagged(HSDen(PSer)));
291  HPDenExp := sum(untagged(PSer).DenFactors.multiplicities);
292  While IsDivisible(HPNum, 1-t) Do
293    HPNum    := HPNum/(1-t);
294    HPDenExp := HPDenExp - 1;
295  EndWhile;
296  return $.Make(HPNum, $.DenMakeStandard(HPDenExp));
297EndDefine;
298
299-------------------------------
300
301Define DenIsStandard(HSDen)
302  fact := 1-indet(RingOf(HSDen.RemainingFactor),1); // 1-t
303  Foreach f In HSDen.factors Do
304    If f<>fact Then return false; EndIf;
305  EndForeach;
306  return True;
307EndDefine; -- DenIsStandard
308
309--------------------------------
310--<<   Hilbert Functions    >>--
311--------------------------------
312
313Define PSerToHilbert(...)
314  TopLevel ERR;
315  If len(ARGV) = 2 Then
316    return $.EvalHilbertFn(ARGV[1],ARGV[2]);
317  Elif
318    len(ARGV) = 1 Then
319    return $.PSerToHilbertFn(ARGV[1]);
320  Else error(ERR.BAD_PARAMS_NUM, ": 1 or 2 expected (Hilbert)");
321  EndIf;
322EndDefine;
323
324-------------------------------
325
326Define AuxHilbertPoly(HSSimplifiedPS)
327  SPS := untagged(HSSimplifiedPS);
328  If IsZero(SPS.num) Then return 0; EndIf;
329  HV  := ToHVec(SPS.num);
330  DIM := sum(SPS.DenFactors.multiplicities);
331  return sum([HV[I]*binomial(indet(RingQQt(1),1)+DIM-I,DIM-1) | I In 1..len(HV)]);
332EndDefine; -- AuxHilbertPoly
333
334
335Define PSerToHilbertPoly(PS)
336  QQt := RingQQt(1);
337  SPS := $.HSSimplified(PS);
338  If untagged(HSDen(SPS)).factors = [] Then return 0; EndIf;
339  return $.AuxHilbertPoly(SPS);
340EndDefine; -- PSerToHilbertPoly
341
342
343Define ToHVec(SpP)
344  If IsZero(SpP) Then return []; EndIf;
345  t := indet(RingOf(SpP),1);
346  Skeleton := [t^i | i In 0..deg(SpP)];
347  HV := coefficients(SpP, Skeleton);
348  return [AsINT(n) | n In HV];
349EndDefine;
350
351Define PSerToHilbertFn(PS)
352  QQt := RingQQt(1);
353  If IsZero(untagged(PS).num) Then return TaggedHilbertFn([ [], 0]); EndIf;
354  SPS := untagged(HSSimplified(PS));
355//  DIM := len(untagged(HSDen(SPS)));
356  DIM := sum(SPS.DenFactors.multiplicities);
357  HV := ToHVec(SPS.num);
358  If DIM=0 Then return TaggedHilbertFn([HV,0]); EndIf;
359  REG := len(HV)-DIM;
360  If REG > 0 Then
361    TrBins := concat( NewList(REG-1,0),
362                   [$.TruncBin(I,DIM-1) | I In (DIM-1)..(DIM+REG)]);
363    HF1 :=[sum([HV[I]*TrBins[N-I+REG+1] | I In 1..(N+1)]) | N In 0..(REG-1)];
364  Else
365    HF1 := [];
366  EndIf;
367  return TaggedHilbertFn([HF1, PSerToHilbertPoly(SPS)]);
368EndDefine; -- PSerToHilbertFn
369
370
371Define PSerRegularityIndex(PS)
372  QQt := RingQQt(1);
373  If IsZero(untagged(PS).num) Then return 0; EndIf;
374  SPS := untagged(HSSimplified(PS));
375  DIM := sum(SPS.DenFactors.multiplicities);
376  HV := ToHVec(SPS.num);
377  return len(HV)-DIM;
378EndDefine; -- PSerRegularityIndex
379
380
381Define EvalHilbertFn(X, N)
382  TopLevel ERR;
383  QQt := RingQQt(1);
384  t := RingElem(QQt, "t");
385  If tag(X)=TagForHilbertFn() Then
386    If N<len(untagged(X)[1]) Then return untagged(X)[1,N+1];
387    Else
388      return AsINT(eval(untagged(X)[2], [RingElem(QQt,N)]));
389    EndIf;
390  EndIf;
391  If tag(X)=TagForPSeries() Then
392    If IsZero(untagged(X).num) Then return 0; EndIf;
393    If $.IsStandard(X) Then     //  standard grading
394      SPS := untagged(HSSimplified(X));
395      HV  := ToHVec(SPS.num);
396      DIM := sum(SPS.DenFactors.multiplicities);
397      If DIM = 0 Then
398	If 0<=N And N<len(HV) Then return HV[N+1]; Else return 0; EndIf;
399      EndIf;
400      return sum([HV[I]*$.TruncBin(N+DIM-I,DIM-1) | I In 1..len(HV)]);
401    EndIf;
402    // univariate and not standard
403    HPDen := untagged(X).DenFactors; //untagged(HSDen(X));
404//    If len(HPDen[1]) = 1 Then
405    If NumIndets(RingOf(HPDen.factors[1])) = 1 Then
406      TruncDen := 1;
407      For i := 1 To len(HPDen.factors) Do
408	w := deg(HPDen.factors[i]);
409	For j:=1 To HPDen.multiplicities[i] Do
410	  TruncDen := TruncDen * sum([t^(k*w) | k in 0..div(N, w)]);
411	  TruncDen := NR(TruncDen, [t^(N+1)]);
412	EndFor;
413      EndFor;
414      TruncNum := NR(untagged(X).num, [t^(N+1)]);
415      TruncHF := NR(TruncNum * TruncDen, [t^(N+1)]);
416      RatResult := CoeffOfTerm(TruncHF, t^N);
417      return AsINT(RatResult);
418    EndIf;
419  EndIf;
420  error(ERR.BAD_PARAMS, "(EvalHilbertFn)");
421EndDefine; -- EvalHilbertFn
422
423
424Define HFRegularityIndex(HF)
425  HFSmall := untagged(HF)[1];
426  HPoly := untagged(HF)[2];
427  For D := len(HFSmall)-1 To 0 Step -1 Do
428    If HFSmall[D+1] <> eval(HPoly, [D]) Then return D+1; EndIf;
429  EndFor;
430  If HPoly = 0 Then error("cannot compute RegularityIndex for 0 function"); EndIf;
431  DD := -1;
432  While True Do
433    If eval(HPoly, [DD])<>0 Then return DD+1; EndIf;
434    DD := DD-1;
435  EndWhile;
436  return 0;
437EndDefine; -- HFRegularityIndex
438
439
440------[   pretty printing   ]--------
441
442Define Tagged(X,T)
443  return tagged(X, T);
444EndDefine;
445
446
447Define TaggedHilbertFn(X)
448  return tagged(X, TagForHilbertFn());
449EndDefine;
450
451Define TaggedHilbertSeries(X)
452  return tagged(X, TagForPSeries());
453EndDefine;
454
455Define TagForHilbertFn() return "HilbertFn"; EndDefine;
456Define TagForPSeries() return "PSeries"; EndDefine;
457
458
459Define Print_HSDen(D, HSDen)
460  Factors := untagged(HSDen).factors;
461  Exp := untagged(HSDen).multiplicities;
462  HSRing := RingOf(untagged(HSDen).RemainingFactor);
463  L := len(Factors);
464  If $.DenIsStandard(HSDen) Then
465    Print "(1-",indet(HSRing, 1),")" On D;
466    If sum(Exp)<>1 Then  Print "^", sum(Exp) On D; EndIf
467  Else
468    If L<>1 Then  Print "( " On D;  EndIf;
469    Print "(1-", LT(Factors[1]), ")" On D;
470    If Exp[1]<>1 Then Print "^", Exp[1] On D; EndIf;
471    For i := 2 To len(Factors) Do
472      Print "*(1-", LT(Factors[i]), ")" On D;
473      If Exp[i]<>1 Then Print "^", Exp[i] On D; EndIf;
474    EndFor;
475    If L<>1 Then  Print " )" On D;  EndIf;
476  EndIf;
477EndDefine;
478
479
480Define Print_PSeries(D, PSer)
481  QQt := RingQQt(1);
482  If $.IsStandard(PSer) Then
483    PSer := $.HSSimplified(PSer);
484  Else
485    PrintLn "---  Non-simplified HilbertPoincare' Series  ---" On D;
486  EndIf;
487  SpPolyNum := $sppoly.PolyToSPPoly(untagged(PSer).num);
488  SpPolyNum := $sppoly.Tagged(reversed(untagged(SpPolyNum)),"SpPoly");
489  Print "(", SpPolyNum, ")" On D;
490  HP_Den := untagged(PSer).DenFactors;
491  If len(untagged(HP_Den).factors)>0 Then Print " / ", tagged(HP_Den, "HSDen") On D;  EndIf;
492EndDefine; -- Print_PSeries
493
494Define Print_HilbertFn(D, HF)
495  For I:=1 To len(HF[1]) Do PrintLn "H(",I-1,") = ",HF[1,I] On D; EndFor;
496  Print "H(t) = ", HF[2], "   for t >= ", len(HF[1]) On D;
497EndDefine;
498
499
500Define TruncBin(A,B)
501  If A<B Or B<0 Then return 0; EndIf;
502  return binomial(A,B);
503EndDefine; -- TruncBin
504
505------[ Hilbert multideg ]-----------------------------------------------
506
507Define HilbertSeriesMultiDeg(QR, WM)
508  return $hp.PoincareMultiDeg(QR, WM);
509EndDefine;
510
511Define PoincareMultiDeg(QR, WM)
512  If type(QR) <> RING And not(IsQuotientRing(QR)) Then
513    error("First argument must be a quotient ring", QR);
514  EndIf;
515  P := BaseRing(QR);
516  If type(WM) <> MAT Then
517    error("Second argument must be a matrix", WM);
518  EndIf;
519  If NumCols(WM) <> NumIndets(P) Then
520    error("Wrong number of weights:" + sprint(NumCols(WM)) + " " + sprint(NumIndets(P)), WM);
521  EndIf;
522  K := CoeffRing(P);
523  M := MakeTermOrd(WM);
524  AuxRing := NewPolyRing(K, SymbolRange("x",1,NumCols(M)), M, NumRows(WM));
525  phi := PolyAlgebraHom(P, AuxRing, indets(AuxRing));
526  PS := HilbertSeries(AuxRing/ideal(apply(phi,gens(DefiningIdeal(QR)))));
527  return PS;
528EndDefine; -- PoincareMultiDeg
529
530
531-- EXAMPLES
532-- Use R ::= QQ[x,y,z];
533-- //PoincareMultideg(ideal(indets())^2, LexMat(NumIndets()));
534
535-- WM := mat([[1,0,0],[1,-1,0]]);
536-- PoincareMultiDeg(R/ideal(indets())^2, WM);
537
538-- WM := mat([[1,0,0],[-1,1,1]]);
539-- PoincareMultiDeg(R/ideal(0), WM);
540
541-- WM := mat([[1,7,0],[0,-5,1]]);
542-- PositiveGrading(WM);
543
544------[ end of Hilbert multideg ]-------------------------------------------
545
546----------------------------------------------------------------------
547PrintTagged := Record[
548		      PSeries := $.Print_PSeries,
549		      HSDen := $.Print_HSDen,
550		      HilbertFn := $.Print_HilbertFn
551		      ];
552
553EndPackage; -- Package $hp
554
555