//load "D:\\magmadoc\\OrbitSizePolFormUnitriangGroup.mag";

                                        ////////////
                                        //Symmetric
                                        ////////////

///////////////////////////////////////////////////////////
//Describes the action of the standard lower unitriangular group
//in SL(n,q) on the Symmetric Tensor of the natural module
///////////////////////////////////////////////////////////

OrbitSizeSymPol:=function(n,p,f)
    q:=p^f;
    m:=Binomial(n+1,2);
    F<w>:=GF(q);
    M:=RMatrixSpace(F,n,n);
    A:=MatrixAlgebra(F,n);
    B:=MatrixAlgebra(F,m);
    e:=func<i,j|MatrixUnit(A,i,j)>;
    bas:=[e(i,i): i in [1..n]] cat [e(i,j)+e(j,i):i in [j+1..n], j in [1..n]];
    W:=sub<M|bas>;
    bas:=Basis(W);
    T:=GL(m,q);
    gens:=[A!1+w^j*e(i+1,i):j in [0..f-1],i in [1..n-1]];
    gensym:=[B!&cat[Coordinates(W,W!(u*b*Transpose(u))):b in bas]:u in gens];
    K:=sub<T|gensym>;
    ok:=Orbits(K);
    os:={*#x:x in ok*};
    P<t>:=PolynomialRing(Integers());
    osp:=&+[Multiplicity(os,q^i)*t^i:i in [0..#os-1]];
    return osp;
end function;

///////////////////////////////////////////////
// We provide here a recursive answer to the same question BUT ONLY if p=char(F) is ODD
///////////////////////////////////////////////

OrbitSizeSymPolRec:=function(n)
    ZZ:=Integers();
    P<q,t>:=PolynomialRing(ZZ,2);
    if n eq 0 then return 1;
    elif n eq 1 then return q;
    else return (1+(q-1)*t^(n-1))*P!$$(n-1)+(q-1)*t^(n-1)*&+[t^i:i in [0..n-2]]*P!$$(n-2);
    end if;
end function;

/*
R<t>:=PolynomialRing(Integers());
Evaluate(OrbitSizeSymPolRec(4),[3,t]);
OrbitSizeSymPol(4,3,1);
*/



                                            //////////////
                                            //Alternating
                                            //////////////

///////////////////////////////////////////////////////////
//Describes the action of the standard lower unitriangular group
//in SL(n,q) on the Exterior Tensor Square of the natural module
///////////////////////////////////////////////////////////

OrbitSizeAltPol:=function(n,p,f)
    q:=p^f;
    m:=Binomial(n,2);
    F<w>:=GF(q);
    M:=RMatrixSpace(F,n,n);
    A:=MatrixAlgebra(F,n);
    B:=MatrixAlgebra(F,m);
    e:=func<i,j|MatrixUnit(A,i,j)>;
    bas:=[e(i,j)-e(j,i):i in [j+1..n], j in [1..n]];
    W:=sub<M|bas>;
    bas:=Basis(W);
    T:=GL(m,q);
    gens:=[A!1+w^j*e(i+1,i):j in [0..f-1],i in [1..n-1]];
    gensym:=[B!&cat[Coordinates(W,W!(u*b*Transpose(u))):b in bas]:u in gens];
    K:=sub<T|gensym>;
    ok:=Orbits(K);
    os:={*#x:x in ok*};
    P<t>:=PolynomialRing(Integers());
    osp:=&+[Multiplicity(os,q^i)*t^i:i in [0..#os-1]];
    return osp;
end function;

///////////////////////////////////////////////
// We provide here a recursive answer to the same question
///////////////////////////////////////////////

OrbitSizeAltPolRec:=function(n)
    ZZ:=Integers();
    P<q,t>:=PolynomialRing(ZZ,2);
    if n eq 0 then return 1;
    elif n eq 1 then return 1;
    else return P!$$(n-1)+(q-1)*t^(n-2)*&+[t^i:i in [0..n-2]]*P!$$(n-2);
    end if;
end function;

/*
R<t>:=PolynomialRing(Integers());
Evaluate(OrbitSizeAltPolRec(4),[3,t]);
OrbitSizeAltPol(4,3,1);
*/


                                ///////////////////////////////
                                //Hermitian any characteristic
                                ///////////////////////////////

///////////////////////////////////////////////////////////
//Describes the action of the standard lower unitriangular group
//in SL(n,q^2) on the Hermitian Tensor Square of the natural module
///////////////////////////////////////////////////////////

OrbitSizeHermitianPol:=function(n,p,f)
    q:=p^f;
    m:=n^2;
    F0:=GF(q);
    F<w>:=GF(q^2);
    M:=RMatrixSpace(F,n,n);
    A:=MatrixAlgebra(F,n);
    B:=MatrixAlgebra(F0,m);
    e:=func<i,j|MatrixUnit(A,i,j)>;
    bas1:=[e(i,i): i in [1..n]] cat [e(i,j)+e(j,i):i in [j+1..n], j in [1..n]] cat [w*e(i,j)+w^q*e(j,i):i in [j+1..n], j in [1..n]];
    Mr,res:=KMatrixSpace(M,GF(q));
    W:=sub<Mr|[res(b):b in bas1]>;
    bas:=Basis(W);
    T:=GL(m,q);
    gens:=[A!1+w^j*e(i+1,i):j in [0..2*f-1],i in [1..n-1]];
    gensconj:=[A!1+w^(q*j)*e(i,i+1):j in [0..2*f-1],i in [1..n-1]];
    gensym:=[B!&cat[Coordinates(W,W!res(gens[u]*b@@res*gensconj[u])):b in bas]:u in [1..#gens]];
    K:=sub<T|gensym>;
    ok:=Orbits(K);
    os:={*#x:x in ok*};
    P<t>:=PolynomialRing(Integers());
    osp:=&+[Multiplicity(os,q^i)*t^i:i in [0..#os-1]];
    return osp;
end function;

///////////////////////////////////////////////
// We provide here a recursive answer to the same question
///////////////////////////////////////////////

OrbitSizeHermitianPolRec:=function(n)
    ZZ:=Integers();
    P<q,t>:=PolynomialRing(ZZ,2);
    if n eq 0 then return 1;
    elif n eq 1 then return q;
    else return (1+(q-1)*t^(2*n-2))*P!$$(n-1)+(q^2-1)*t^(2*n-3)*&+[t^(2*i):i in [0..n-2]]*P!$$(n-2);
    end if;
end function;

/*
R<t>:=PolynomialRing(Integers());
Evaluate(OrbitSizeHermitianPolRec(4),[3,t]);
OrbitSizeHermitianPol(4,3,1);
*/

                                    /////////////////////////////////////////
                                    //Recursive Symmetric Even Characteristic
                                    /////////////////////////////////////////

OrbitSizeSymPolEvenPartRec:=function(n,j)
    ZZ:=Integers();
    P<q,t>:=PolynomialRing(ZZ,2);
    if j lt 0 or j gt n then return P!0;
    elif j eq 0 then return P!OrbitSizeAltPolRec(n);
    elif j eq 1 then return (q-1)*t^(n-1)*&+[P!$$(n-1,i):i in [0..n-1]];
    else tot:=P!$$(n-1,j-1)+(q-1)^2*t^(2*n-2-j)*P!OrbitSizeAltPolRec(n-2);
        if j gt 2 then tot +:=(q-1)*&+[t^(2*n-2-i):i in [2..j-1]]*P!$$(n-2,j-2);
        end if;
        if j le n-1 then tot +:=(q-1)^2*t^(2*n-2-j)*&+[P!$$(n-2,i):i in [j-1..n-2]]+
                            (q-1)*&+[t^(2*n-1-i):i in [j+1..n]]*P!$$(n-2,j-1);
        end if;
    return tot;
    end if;
end function;

OrbitSizeSymPolEvenRec:=function(n)
    ZZ:=Integers();
    P<v,t>:=PolynomialRing(ZZ,2);
    R<u>:=PolynomialRing(P);
    return &+[R!(P!Evaluate(OrbitSizeSymPolEvenPartRec(n,j),[v+1,t])*u^j):j in [0..n]];
end function;

/*
pol:=OrbitSizeSymPolEvenRec(3);
coe:=Coefficients(pol);
str:="";
for j in [1..#coelis] do
    tmp:="";
    for i in [1..#coelis[j]] do
        tmp cat:=Sprint(coelis[j,i]) cat ")t^" cat IntegerToString(i-1) cat "+";
    end for;
    str cat:="(" cat tmp cat ")u^" cat IntegerToString(j-1) cat "+";
end for;
*/

OrbitSizeQuadPol:=function(n,f)
    q:=2^f;
    m:=Binomial(n+1,2);
    F<w>:=GF(q);
    M:=RMatrixSpace(F,n,n);
    A:=MatrixAlgebra(F,n);
    B:=MatrixAlgebra(F,m);
    bi:=func<i,j|((2*n-i)*(i-1) div 2)+j>;
    e:=func<i,j|MatrixUnit(B,i,j)>;
    gen:=function(r,s)
        mat:=B!1+w^(2*r)*e(bi(s-1,s-1),bi(s,s))+w^r*&+[e(bi(s-1,j),bi(s,j)):j in [s..n]];
        if s gt 2 then mat +:=w^r*&+[e(bi(i,s-1),bi(i,s)):i in [1..s-2]];
        end if;
        return mat;
    end function;
    gensym:=[gen(r,s):r in [0..f-1],s in [2..n]];
    T:=GL(m,q);
    K:=sub<T|gensym>;
    ok:=Orbits(K);
    os:={*#x:x in ok*};
    P<t>:=PolynomialRing(Integers());
    osp:=&+[Multiplicity(os,2^i)*t^i:i in [0..#os*f-1]];
    return osp;
end function;

ConvertToQuadric:=function(vec,n,f)
    q:=2^f;
    F<w>:=GF(q);
    R:=PolynomialRing(F,n);
    AssignNames(~R,["x" cat IntegerToString(i):i in [1..n]]);
    bi:=func<i,j|((2*n-i)*(i-1) div 2)+j>;
    pol:=&+[vec[bi(i,j)]*R.i*R.j:j in [i..n],i in [1..n]];
    return pol;
end function;

ConvertToAlternatingForm:=function(vec,n,f)
    q:=2^f;
    F<w>:=GF(q);
    R:=PolynomialRing(F,n);
    A:=MatrixAlgebra(F,n);
    AssignNames(~R,["x" cat IntegerToString(i):i in [1..n]]);
    bi:=func<i,j|((2*n-i)*(i-1) div 2)+j>;
    e:=func<i,j|MatrixUnit(A,i,j)>;
    pol:=&+[vec[bi(i,j)]*(e(i,j)+e(j,i)):j in [i+1..n],i in [1..n]];
    return pol;
end function;


HammingWeight:=function(vec)
    dim:=Dimension(Parent(vec));
    return &+[1:i in [1..dim]|vec[i] ne 0];
end function;


OrbitSizeQuadPolRec:=function(n,f,i)
    ZZ:=Integers();
    P<q,t>:=PolynomialRing(ZZ,2);
    if i lt 0 or i gt n then return P!0;
    elif n eq 0 and i eq 0 then return P!1;
    elif n eq 1 and i eq 0 then return P!q;
    elif n ge 2 and i eq 2 then return 2*(q-1)^2*t^(f*(2*n-3)-1)*P!$$(n-2,f,0);
    elif n ge 2 and i ge 3 then
        tot:=t^f*P!$$(n-1,f,i-1);
        if n ge 3 then
            tot +:=(q-1)^3*t^(f*(3*n-4-i))*Evaluate(P!OrbitSizeAltPolRec(n-3),[q,t^f]);
        end if;
        if i ge 4 then
            tot +:=(q-1)*&+[t^(f*(2*n-1-j)):j in [3..i-1]]*P!$$(n-2,f,i-2);
        end if;
        if n ge i+1 then
            tot +:=(q-1)*&+[t^(f*(2*n-j)):j in [i+1..n]]*P!$$(n-2,f,i-1)+
                    (q-1)^2*t^(f*(2*n-1-i))*&+[P!$$(n-2,f,j-2):j in [i+1..n]];
        end if;
        return tot;
    elif n ge 2 and i eq 0 then
    return P!$$(n-1,f,0)+(q-1)*t^(f*(n-1))*Evaluate(P!OrbitSizeAltPolRec(n-1),[q,t^f])+
            (q-1)*&+[t^(f*(2*n-1-j)):j in [2..n]]*P!$$(n-2,f,0)+&+[P!$$(n,f,j):j in [2..n]];
    end if;
end function;

OrbitSizeQuadPolRecDense:=function(n,i)
    ZZ:=Integers();
    P<v,t,u>:=PolynomialRing(ZZ,3);
    R:=Parent(OrbitSizeAltPolRec(2));
    inj:=hom<R->P|[P.i:i in [1..2]]>;
    if i lt 0 or i gt n then return P!0;
    elif n eq 0 and i eq 0 then return P!1;
    elif n eq 1 and i eq 0 then return P!v+1;
    elif n ge 2 and i eq 2 then return 2*u*v^2*t^(2*n-3)*P!$$(n-2,0);
    elif n ge 2 and i ge 3 then
        tot:=t*P!$$(n-1,i-1);
        if n ge 3 then
            tot +:=v^3*t^(3*n-4-i)*Evaluate(inj(OrbitSizeAltPolRec(n-3)),[v+1,t,u]);
        end if;
        if i ge 4 then
            tot +:=v*&+[t^(2*n-1-j):j in [3..i-1]]*P!$$(n-2,i-2);
        end if;
        if n ge i+1 then
            tot +:=v*&+[t^(2*n-j):j in [i+1..n]]*P!$$(n-2,i-1)+
                    v^2*t^(2*n-1-i)*&+[P!$$(n-2,j-2):j in [i+1..n]];
        end if;
        return tot;
    elif n ge 2 and i eq 0 then
    return P!$$(n-1,0)+v*t^(n-1)*Evaluate(inj(OrbitSizeAltPolRec(n-1)),[v+1,t,u])+
            v*&+[t^(2*n-1-j):j in [2..n]]*P!$$(n-2,0)+&+[P!$$(n,j):j in [2..n]];
    end if;
end function;

/*
detailed check in case n=3, f=2 we analyze different contributions theoretically and directly
n:=3;
f:=2;
q:=2^f;
    m:=Binomial(n+1,2);
    F<w>:=GF(q);
    M:=RMatrixSpace(F,n,n);
    A:=MatrixAlgebra(F,n);
    B:=MatrixAlgebra(F,m);
    bi:=func<i,j|((2*n-i)*(i-1) div 2)+j>;
    e:=func<i,j|MatrixUnit(B,i,j)>;
    gen:=function(r,s)
        mat:=B!1+w^(2*r)*e(bi(s-1,s-1),bi(s,s))+w^r*&+[e(bi(s-1,j),bi(s,j)):j in [s..n]];
        if s gt 2 then mat +:=w^r*&+[e(bi(i,s-1),bi(i,s)):i in [1..s-2]];
        end if;
        return mat;
    end function;
    gensym:=[gen(r,s):r in [0..f-1],s in [2..n]];
    T:=GL(m,q);
    K:=sub<T|gensym>;
    ok:=Orbits(K);
    os:={*#x:x in ok*};
    P<t>:=PolynomialRing(Integers());
    osp:=&+[Multiplicity(os,2^i)*t^i:i in [0..#os*f-1]];
    ok1step:=[o:o in ok|&or[&and[v[i] eq 0:i in [1..3]]:v in o]];
    {*#o:o in ok1step*};
    ok2step:=[o:o in ok|&or[&and[v[i] eq 0:i in [2..3]] and v[1] ne 0:v in o]];
    {*#o:o in ok2step*};
    ok3step:=[o:o in ok|&or[&or[v[i] ne 0:i in [2..3]] and v[1] eq 0:v in o]];
    {*#o:o in ok3step*};
    ok4step:=[o:o in ok|&or[v[1] ne 0 and v[2] eq 0 and v[3] ne 0 and v[4] eq 0 and v[5] eq 0:v in o]];
    {*#o:o in ok4step*};
*/

/*
ltu:=function(n)
    pol:=OrbitSizeQuadPolRecDense(n,0);
    R:=Parent(pol);
    ldt:=LeadingTerm(pol,R.3);
    coeu:=Coefficients(ldt,R.3);
    pul:=coeu[#coeu];
    coet:=Coefficients(pul,R.2);
    coev:=Coefficients(pul,R.1);
    return pul;//,coet,coev;
end function;

ltt:=function(n)
    pol:=OrbitSizeQuadPolRecDense(n,0);
    R:=Parent(pol);
    ldt:=LeadingTerm(pol,R.2);
    coeu:=Coefficients(ldt,R.2);
    pul:=coeu[#coeu];
    coet:=Coefficients(pul,R.2);
    coev:=Coefficients(pul,R.1);
    return pul;//,coet,coev;
end function;

n:=9;
m:=Floor(n/2);
for i in [2..n] do
pol:=OrbitSizeQuadPolRecDense(n,i);
P:=Parent(pol);
Coefficient(pol,P.3,m);
end for;


*/