Z:=IntegerRing();

/* factorize n into prime powers as u1^a1*u2^a2*u3^a3 */

A:=Factorization(n);
u:=[1,1,1];
a:=[0,0,0];
for i:=1 to #A do
    u[i]:=A[i][1];
    a[i]:=A[i][2];
end for;
u1:=u[1];
u2:=u[2];
u3:=u[3];
a1:=a[1];
a2:=a[2];
a3:=a[3];


/* given n and q, compute the number of q-linear subspaces of F_q^n */

N:=function(n,q);
    x:=0;
    for k:=1 to n do
        y:=1;
        for i:=0 to k-1 do
            y:=y*(q^(n-i)-1)/(q^(i+1)-1);
        end for;
        x:= x+y;
    end for;
    return x+1;
end function;


/* size of intermediate field F_q^d, where d=u1^j1*u2^j2*u3^j3 */

d:=function(j1,j2,j3);
x:=q^((u1^j1)*(u2^j2)*(u3^j3));
return x;
end function;


/* calculate the coefficient L for each intermediate field F_q^d such that x^L is in F_q^d, for any primitive element x in F_q^n */

L:=function(j1,j2,j3);
y:=(q^n-1)/(d(j1,j2,j3)-1);
return y;
end function;


/* Given the L coefficients and the trace representation of the cyclic code, compute the possible indices (see Thm. 2.6) */

IN:=function(k,j1,j2,j3);
p:=(LeastCommonMultiple(Z!k,Z!L(j1,j2,j3)))/k;
return p;
end function;


/* The number of nonzero subspaces of F_q^n that are maximally defined over F_q^d (see Thm. 3.3)*/

M:=function(j1,j2,j3);

if ((j1 ne a1) and (j2 ne a2) and (j3 ne a3)) then

w:=N((u1^(a1-j1))*(u2^(a2-j2))*(u3^(a3-j3)),q^((u1^j1)*(u2^j2)*(u3^j3)))-N(u1^(a1-j1-1),q^(u1^(j1+1)))-N(u2^(a2-j2-1),q^(u2^(j2+1)))-N(u3^(a3-j3-1),q^(u3^(j3+1)))+N((u1^(a1-j1-1))*(u2^(a2-j2-1)),q^((u1^(j1+1))*(u2^(j2+1))))+N((u2^(a2-j2-1))*(u3^(a3-j3-1)),q^((u2^(j2+1))*(u3^(j3+1))))+N((u1^(a1-j1-1))*(u3^(a3-j3-1)),q^((u1^(j1+1))*(u3^(j3+1))))-N((u1^(a1-j1-1))*(u2^(a2-j2-1))*(u3^(a3-j3-1)),q^((u1^(j1+1))*(u2^(j2+1))*(u3^(j3+1))));

else if ((j1 eq a1) and (j2 ne a2) and (j3 ne a3)) then

w:=N((u1^(a1-j1))*(u2^(a2-j2))*(u3^(a3-j3)),q^((u1^j1)*(u2^j2)*(u3^j3)))-N(u2^(a2-j2-1),q^(u2^(j2+1)))-N(u3^(a3-j3-1),q^(u3^(j3+1)))+N((u2^(a2-j2-1))*(u3^(a3-j3-1)),q^((u2^(j2+1))*(u3^(j3+1))));

else if ((j1 ne a1) and (j2 eq a2) and (j3 ne a3)) then

w:=N((u1^(a1-j1))*(u2^(a2-j2))*(u3^(a3-j3)),q^((u1^j1)*(u2^j2)*(u3^j3)))-N(u1^(a1-j1-1),q^(u1^(j1+1)))-N(u3^(a3-j3-1),q^(u3^(j3+1)))+N((u2^(a1-j1-1))*(u3^(a3-j3-1)),q^((u1^(j1+1))*(u3^(j3+1))));

else if ((j1 ne a1) and (j2 ne a2) and (j3 eq a3)) then

w:=N((u1^(a1-j1))*(u2^(a2-j2))*(u3^(a3-j3)),q^((u1^j1)*(u2^j2)*(u3^j3)))-N(u1^(a1-j1-1),q^(u1^(j1+1)))-N(u3^(a2-j2-1),q^(u2^(j2+1)))+N((u2^(a1-j1-1))*(u2^(a2-j2-1)),q^((u1^(j1+1))*(u2^(j2+1))));

else if ((j1 eq a1) and (j2 eq a2) and (j3 ne a3)) then

w:=N((u1^(a1-j1))*(u2^(a2-j2))*(u3^(a3-j3)),q^((u1^j1)*(u2^j2)*(u3^j3)))-N(u3^(a3-j3-1),q^(u3^(j3+1)));

else if ((j1 eq a1) and (j2 ne a2) and (j3 eq a3)) then

w:=N((u1^(a1-j1))*(u2^(a2-j2))*(u3^(a3-j3)),q^((u1^j1)*(u2^j2)*(u3^j3)))-N(u3^(a2-j2-1),q^(u2^(j2+1)));

else if ((j1 ne a1) and (j2 eq a2) and (j3 eq a3)) then

w:=N((u1^(a1-j1))*(u2^(a2-j2))*(u3^(a3-j3)),q^((u1^j1)*(u2^j2)*(u3^j3)))-N(u1^(a1-j1-1),q^(u1^(j1+1)));

else w:= N(1,q^n);

end if;
end if;
end if;
end if;
end if;
end if;
end if;

return w;
end function;


/* Create two matrices of size cols x rows to store the subcode indices and their multiplicities, where cols is the number of basic zeros and rows is the number of intermediate fields */

cols:=#i;
rows:=NumberOfDivisors(n);
Imat := Matrix(IntegerRing(), rows, cols, []);
Ml := Matrix(IntegerRing(), rows, 1, []);

/* Calculate the subcode indices and their multiplicities for a cyclic code with one basic zero (simplex case) */

I:={};
index:=0;

for j1:=0 to a1 do
for j2:=0 to a2 do
for j3:=0 to a3 do

    index:=index+1;
    
    if cols ge 2 then
        Ml[index, 1]:=N((u1^(a1-j1))*(u2^(a2-j2))*(u3^(a3-j3)),q^((u1^j1)*(u2^j2)*(u3^j3)));
    else
        Ml[index, 1]:= M(j1,j2,j3);
    end if;

    for z:=1 to cols do

        l:=IN(i[z],j1,j2,j3);
        I join:={l};
        Imat[index, z]:=l;

    end for;
end for;
end for;
end for;

I:=Exclude(I,q^n-1);

/* For the cases of 2 or 3 basic zeros, create two more matrices as above to keep the LCM of the indices and the multiples of the appearances */

doublesPerValue:=(cols-1)*rows;
LCMdoubles := Matrix(IntegerRing(), rows, cols*doublesPerValue, []);
Mdoubles := Matrix(IntegerRing(), rows, cols*doublesPerValue, []);

triplesPerValue:=rows^2;
LCMtriples := Matrix(IntegerRing(), rows, cols*triplesPerValue, []);
Mtriples := Matrix(IntegerRing(), rows, cols*triplesPerValue, []);


/* Compute the subcode indices and the subspace multiplicities for a cyclic code with two basic zeros*/

for z:=1 to cols do

    for k:=1 to rows do
        
        for offset:=1 to cols-1 do
        
            col:=((z+offset-1) mod cols) + 1;
        
            for row:=1 to rows do
                index:=(offset-1)*rows + row;
                LCMdoubles[k, (z-1)*doublesPerValue + index] := LeastCommonMultiple(Imat[k,z], Imat[row, col]);
                Mdoubles[k, (z-1)*doublesPerValue + index] := Ml[k, 1]*Ml[row, 1];
            end for;
            
        end for;

        /* Calculate the subcode indices and the subspace multiplicities for a cyclic code with three basic zeros*/
        
        col1:=((z+1-1) mod cols) + 1;
        col2:=((z+2-1) mod cols) + 1;
        
        for row1:=1 to rows do
        
            tmpLCM :=LeastCommonMultiple(Imat[k,z], Imat[row1, col1]);
            tmpM := Ml[k, 1]*Ml[row1, 1];
        
            for row2:=1 to rows do
            
                index:=(row1-1)*rows + row2;
                LCMtriples[k, (z-1)*triplesPerValue + index] := LeastCommonMultiple(tmpLCM, Imat[row2, col2]);
                Mtriples[k, (z-1)*triplesPerValue + index] := tmpM*Ml[row2, 1];
                
            end for;
            
        end for;

    end for;

end for;

/* Assign the suitable pair of matrices with respect to the number of basic zeros*/

if (cols gt 2) then
    LCMvalues :=LCMtriples;
    Mvalues:=Mtriples;
    cols:=cols*triplesPerValue;
elif (cols eq 2) then
    LCMvalues :=LCMdoubles;
    Mvalues:=Mdoubles;
    cols:=cols*doublesPerValue;
elif (cols eq 1) then
    LCMvalues :=Imat;
    Mvalues:=Ml;
end if;

/* For the indices stored in the first matrix, assign the maximum value among all the multiplicities in the second matrix*/

Ival :=Sort([val : val in I]);
Mmax := [];
nbrI := #Ival;

for index:=1 to nbrI do

    tmp := -1;

    for row:=1 to rows do

        for col:=1 to cols do
            
            if ((Ival[index] eq LCMvalues[row, col]) and (Mvalues[row, col] ge tmp)) then
                tmp:=Mvalues[row, col];
            end if;

        end for;

    end for;
    
    Mmax:=Append(Mmax, tmp);
    
end for;

/* Find index of value inside a sequence, only if it is between indices Start and Stop    */

indexOf:=function(value, sequence: Start:=1, Stop:=-1);
    if (Stop lt Start) then
        Stop := #sequence;
    end if;
    
    result := -1;

    for index:=Start to Stop do
        if (sequence[index] eq value) then
            result := index;
            break;
        end if;
    end for;
    
    return result;
end function;

/* Inclusion-Exclusion function for the subspace multiplicities, if there are more than one basic zeros*/

IncEx := function(sequence, index);
    
    Factors := [0 : i in [1..#sequence]];
    Factors[index] := 1;
    divs := Divisors(Z!sequence[index]);  /* Find divisors of the current subcode index */
    nbrDivs:=#divs - 1;  /*    Always skip last divisor since it is the index itself */
    
    count:=[1 : i in [1..nbrDivs]];
    factor:=[0 : i in [1..nbrDivs]];
    
    for divIndex:=nbrDivs to 1 by -1 do
    
        currentDiv:=divs[divIndex];
        factor[divIndex] := -count[divIndex];  /* Use as factor the divisor's count */
        
        /*    If the divisor is not prime, it may need to be included or excluded */
        if (factor[divIndex] ne 0) then
    
            divIndexInSequence := indexOf(currentDiv, sequence: Start:=1, Stop:=index-1);  /* Check if the divisor is in the sequence */
            
            if (divIndexInSequence ne -1) then  /*    If so, include or exclude it */
            
                Factors[divIndexInSequence] := factor[divIndex]; /*    Store the corresponding factor for future processing */
                subDivs := Divisors(divs[divIndex]); /*    Find the divisors of this divisor */
                nbrSubDivs:=#subDivs - 1;   /* Always skip last sub-divisor since it is the divisor itself */
                
                /*    Update the count of the other divisors corresponding to a sub-divisor */
                for subDivIndex:=nbrSubDivs to 1 by -1 do
                
                    currentSubDiv:=subDivs[subDivIndex];
                    subDivIndexInDivs := indexOf(currentSubDiv, divs: Start:=1, Stop:=divIndex-1);
                    
                    if (subDivIndexInDivs ne -1) then
                        count[subDivIndexInDivs]:=count[subDivIndexInDivs]+factor[divIndex];
                    end if;
                    
                end for;
                
            end if;
            
        end if;
    
    end for;
    
    return Factors;
end function;


Mcor := [0 : i in [1..nbrI]];  /* Initialise subcode multiplicities for inclusion-exclusion    */

Mcor[1]:=Mmax[1]-N(1,q^n-1); /* Number of cyclic subcodes */

/*    Number of subcodes for the quasi-cyclic cases */

if cols ne 1 then

    for index:=2 to nbrI do
    
        factors := IncEx(Ival, index);

        for factorIndex:=1 to index do
            Mcor[index]:=Mcor[index]+factors[factorIndex]*Mmax[factorIndex];
        end for;

    end for;
    
end if;

P:={};
for i:=1 to nbrI do
    P join:={[Ival[i],Mcor[i]]};
end for;


/* Check the sizes of q-cyclotomic cosets for each basic zero and return the result*/

C:={};
for j:=1 to #i do
     w := [i[j]*q^k mod (q^n-1) : k in [0..n-1] ];
    C join:={w[i]: i in [1..#w]};
end for;
if #C eq (#i)*n then
    print "The indices of the QC subcodes are:"; print I;
    print "The numbers of the subcodes corresponding to each index are:";print P;
else
    print "Some of the basic zeros have q-cyclotomic coset size less than n, please try with other exponents!";
end if;