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