/* ---------------------------------------------------------------------- PROCEDURES ---------------------------------------------------------------------- */ LIB "zeroset.lib"; LIB "ring.lib"; LIB "poly.lib"; LIB "linalg.lib"; LIB "general.lib"; LIB "rinvar.lib"; LIB "multigrading.lib"; proc ComputeQP(W){ int K= size(W); int d = lcm(W); int M=0; for (int j=1; j<=K; j=j+1){ M= gcd(M, W[j]); } int d1 = d div M; int i; intvec W1; for (i=1; i<=K; i=i+1){ W1[i]= W[i] div M; } minpoly = rootofUnity(d1); /* We have used the version 3-1-6 of Singular. In this version, the function rootofUnity(n) doesn't return the correct primitive root if the argument n is 2*p*q, for every p, q primes distinct. In the later versions, the bug is been fixed. */ //minpoly = a8+a7-a5-a4-a3+a+1; /* 2*3*5 */ //minpoly = a12+a11-a9-a8+a6-a4-a3+a+1; /* 2*3*7 */ //minpoly = a10-a9+a8-a7+a6-a5+a4-a3+a2-a+1; /* 2*3*11 */ //minpoly = a24+a23-a19-a18-a17-a16+a14+a13+a12+a11+a10-a8-a7-a6-a5+a+1; /* 2*5*7 */ list Period; list Coeff; intvec W_aus; int t; int g; int period2=0; for (i=1; i<=K; i=i+1){ W_aus= W1; W_aus[i]=0; g=0; t=1; while (g!=1 and t<=K) { g = gcd(g,W_aus[t]); t = t+1; } period2 = lcm(period2, g); } Period[1]= 1; Coeff[1] = 1 / (product(W1) * factorial(K-1)); if (period2 == 1){ Period[2] = 1; Coeff[2]= sum(W1) / (2 * product(W1) * factorial(K-2)); } else{ Period[2]=0; Coeff[2]=0; } list L1; L1 = HilbertQuasiPolynomials(W1,Coeff,Period); list L; if (M == 1){ L= L1; } else{ for (i=1; i<=d; i=i+1){ L[i]=0; } poly q; for (i=1; i<=d1; i=i+1){ q=0; for (j=1; j<=K; j=j+1){ q= q + L1[i][j]/M^(K-j); } L[(M*i)] = q; } } return (L); } proc ComputeQPQuotient(W,I,L){ list HP = hilb(I,1,W); //It returns the numerator of not reduced Hilbert-PoincarĂ¨ series. list L_I; poly p; int i,j; int r; int d = lcm(W); for (i=1; i<=d; i=i+1){ p=0; for (j=1; j<= size(HP[1]); j=j+1){ r = i-j+1 mod d; while ( r <= 0 ){ r = r+d; } p= p+ HP[1][j]* (subst(L[r],x,x-j+1)); } L_I[i]=p; } return (L_I); } proc Vandermonde(M){ int K= size(M); matrix V[K][K]; int i,j; for (i=1; i<=K; i=i+1){ for (j=1; j<=K; j=j+1){ V[i,j]=(M[i])^(j-1); } } return (V); } proc FactorizePoly (W){ int K= size(W); int d = lcm(W); int i, j; list L; //In L[i] there is the number of times that the ith power of dth unity root appears in the factorization of the denominator of Hilbert-PoincarĂ© series associated to R. for (i=1; i<=d; i=i+1){ L[i]=0; } for (i=1; i<=K; i=i+1){ for (j=1; j<=W[i]; j=j+1){ L[(j*d)div W[i]]= L[(j*d) div W[i]] +1; } } return (L); } proc PartialSum (W,n,a){ list L = FactorizePoly (W); list PS; poly coeff1 =0; int i; for (int r=1; r<=n; r=r+1){ coeff1 =0; for (i=1; i<=size(L); i=i+1){ coeff1 = coeff1 + L[i]* (a^(i*r)); } PS[r]=coeff1; } return (PS); } proc FunctionHilbert(PS,n,a){ poly coeff=0; list F=1; int r; for (int j=1; j<=n; j=j+1){ coeff=0; for (r=1; r<=j; r=r+1){ coeff= coeff + PS[r] * F[j-r+1]; } F[j+1]= coeff/j; } return(F); } proc HilbertPoly (W,F,j,Coeff,Period){ poly x = var(1); int K= size(W); int d = lcm(W); list L = FactorizePoly (W); int i; list A; list B; //It contains the values of Hilbert function that need for the interpolation. int b=j; int s; if (K-(Period[1]+Period[2]) > 0){ for (i=0; i<= K-(2+Period[2]); i=i+1){ s= b + i*d; A = insert (A,s,size(A)); B = insert (B,F[s+1]- Coeff[1]*(s^(K-1)) - Coeff[2]*(s^(K-2)),size(B)); } matrix V = Vandermonde(A); matrix Y[K-(Period[1]+Period[2])][1]; for (i=1; i<=K-(Period[1]+Period[2]); i=i+1 ){ Y[i,1]=B[i]; } matrix D= inverse(V)*Y; } poly p= Coeff[1] * x^(K-1) + Coeff[2] * x^(K-2); for(i=1; i<=K-(Period[1]+Period[2]); i=i+1){ p= p+ D[i,1]* x^(i-1); } return (p); } proc HilbertQuasiPolynomials(W,Coeff,Period){ int K= size(W); int d = lcm(W); int i; int M=0; for (i=1; i<=K; i=i+1){ M= gcd(M, W[i]); } list PS = PartialSum(W, d*(K-(Period[1]+Period[2])),a); list F = FunctionHilbert(PS,d*(K-(Period[1]+Period[2])),a); list L; for (i=1; i<=d; i=i+1){ if( gcd(M,i)== M){ L= insert(L, HilbertPoly(W,F,i,Coeff,Period), size(L)); } else{ L= insert(L, 0, size(L)); } } return(L); }