LIB "zeroset.lib";
LIB "ring.lib";
LIB "poly.lib";
LIB "linalg.lib";
LIB "general.lib";
LIB "rinvar.lib";
LIB "multigrading.lib";



intvec W = 1,1; //weights vector

ring R = (0,a), (x,y), wp(W); //polynomial ring graded by W

int K= size(W);
int d = 1;
for (int i=1; i<=K; i=i+1){
    d= lcm(d, W[i]);
}


int M=0;
for (int j=1; j<=K; j=j+1){
    M= gcd(M, W[j]);
}

int d1 = d div M;


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;
    }
}

L;


ideal I = (y^2); //ideal

I = std(I);
list HP = hilb(I,1,W);

list L_I;
poly p;

int r;

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;
}


L_I;



/* ----------------------------------------------------------------------
 PROCEDURES
 ---------------------------------------------------------------------- */




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 = 1;
    int i, j;
    
    for (i=1; i<=K; i=i+1){
        d= lcm(d, W[i]);
    }
    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 (C,PS,n,a){
    
    poly coeff =0;
    
    for (int r=1; r<=n; r=r+1){
        coeff= coeff + PS[r] * C[n-r+1];
    }
    return(coeff/n);
    
}



proc HilbertPoly (W,PS,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.
    list C=1; //It is an auxiliar list in which there are the values of Hilbert function.
    int b=j;
    
    while(b < K-1){
        b = b+d;
    }
    for (int t=1; t<=b; t=t+1){
        
        C[t+1]= FunctionHilbert(C,PS,t,a);
    }
    int s;
    if (K-(Period[1]+Period[2]) > 0){
    for (i=0; i<= K-(2+Period[2]); i=i+1){
        s= b + i*d;
        for (t= s-i*d; t<=s; t=t+1){
            C[t+1]= FunctionHilbert(C,PS,t,a);
        }
        A = insert (A,s,size(A));
        B = insert (B,C[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+1)*(K-1),a);
    list L;
    for (i=1; i<=d; i=i+1){
        if( gcd(M,i)== M){
            L= insert(L, HilbertPoly(W,PS,i,Coeff,Period), size(L));
        }
        else{
            L= insert(L, 0, size(L));
        }
    }
    return(L);
}