#The computation of $P(x)$ and $\tilde{P}(s)$
# define the groups of orders 108 and 216 
K = CyclotomicField(9);
a = K.gen();
w = a^6;
e = a^4;
l = w - w^2;
s = matrix(K, [[1,0,0],[0,w,0],[0,0,w^2]])
t = matrix(K, [[0,1,0],[0,0,1],[1,0,0]])
u = matrix(K,[[e,0,0],[0,e,0],[0,0,e*w]])
v = matrix(K,[[1/l,1/l,1/l],[1/l,w/l,w^2/l],[1/l,w^2/l,w/l]])
G = MatrixGroup(s,t,u,v);
G216 = G.subgroup([s,t,u * v * u.inverse(),v]);
G108 = MatrixGroup(s,t,v);
sum108=0;

#compute P(x) and tildeP(x)
for a in G108:
    A = (a).matrix();
    sum108 = sum108 - 1/A.charpoly();
print("P(x) =",sum108);
sum216=0;

for a in G216:
    A = (a).matrix();
    sum216 = sum216 - 1/A.charpoly();
print("tildeP(x) =",sum216)











#moments of Kl3---------------------------------------------

p=3; #choose a prime p

import math;
#define Kl3
def Kl(a):
    import sage.rings.all as rings
    sum = 0;
    for x in range(1,p):
        for y in range(1,p):
            sum = sum + exp (2 * pi * I/p * (x + y + (a * x * y).inverse_mod(p)));
    return sum

def s2(a):
    v = conjugate(CC(Kl(a))) * p;
    return v
                
# moment of Sym^3Kl_3
mom3=0;
for a in range(1,p):
    a = ZZ(a);
    mom3 = mom3 + CC(Kl(a))^3 - 2 * CC(Kl(a)) * CC(s2(a)) + CC(p)^3;
a3= -(mom3 + 1 + p^2)
print(a3)

# moment of Sym^4Kl_3
mom4 = 0;
for a in range(1,p):
    a = ZZ(a);
    mom4 = mom4 + CC(Kl(a))^4 - 3 * CC(Kl(a))^2 * CC(s2(a)) + CC(s2(a))^2 + 2 * CC(Kl(a)) * CC(p)^3;
a4 = -p^(-3) * (mom4 + 1 + p^2 + p^4)
print(a4)

# moment of Sym^6Kl_3
mom6=0;
for a in range(1,p):
    a = ZZ(a);
    mom6 = mom6 + CC(Kl(a))^6 + 4 * CC(Kl(a))^3 * p^3 + p^6 - 5 * CC(Kl(a))^4 * CC(s2(a)) - 6 * CC(Kl(a)) * p^3 * CC(s2(a)) + 6 * CC(Kl(a))^2 * CC(s2(a))^2 - CC(s2(a))^3;
print(mom6)

# moment of Kl_{SL_3}^{V_{2,1}}
mom21 = 0;
for a in range(1,p):
    a = ZZ(a);
    mom21 = mom21 +(CC(Kl(a))^2 * CC(s2(a)) - CC(s2(a))^2 - p^3 * CC(Kl(a)));
a21 = -p^(-4) * (mom21 + p + p^2 + p^3)
print(a21)

# moment of Kl_{SL_3}^{V_{2,2}}
mom22=0;
for a in range(1,p):
    a = ZZ(a);
    mom22 = mom22 + (CC(Kl(a))^2 - CC(s2(a))) * (CC(s2(a))^2 - p^3 * CC(Kl(a)))- p^3 * CC(Kl(a)) * CC(s2(a))   ;
a22 = - p^(-5) * (mom22 + p^2 + p^3 + 2 * p^4 + 2 * p^6)
print(a22) 

# when N*k^2 is too big 
M = ModularSymbols(Gamma0(N),k,sign=1,base_ring=GF(13))
Hp = M.new_subspace().cuspidal_subspace().hecke_polynomial(p)
print (Hp(a22))








#moments of Kl4---------------------------------------------

#choose a prime
p=7;

import math

#compute Kl_4(a;p), as a list Klp
k = GF(p,'r')
Klp = [-1];
for a in range(1,p):
    klp = 0;
    for x in k:
        for y in k:
            for z in k:
                if x*y*z!=0:
                    u = x + y + z + a/(x * y * z);
                    u = u.lift();
                    klp = klp + exp(2 * pi * I/p * u);
                    klp = n(klp);
    Klp.append(klp);

#compute Kl_4(a;p^2), as a list Klp2
K = GF(p^2,'r')
Klp2 = [-1];
for a in range(1,p):
    klp2 = 0;
    for x in K:
        for y in K:
            for z in K:
                if x*y*z!=0:
                    u = x + y + z + a/(x * y * z);
                    u = (u.trace());
                    u = u.lift();
                    klp2 = klp2 + exp(2 * pi * I/p * u);
                    klp2 = n(klp2);
    Klp2.append(klp2)

#compute s_1(a) and s_2(a) as lists S1 and S2
S1 = [1];
S2 = [1];
for a in range(1,p):
    s1 = - Klp[a];
    s2 = - Klp2[a];
    s2 = (s1^2 - s2)/2;
    S1.append(s1);
    S2.append(s2);

#compute moments 

#compute moments m
mom3=0
for a in range(1,p):
    mom3 = mom3 + S1[a]^3 - 2 * S1[a] * S2[a] + p^3 * conjugate(S1[a]);

mom4=0;
for a in range(1,p):
    mom = mom + S1[a]^4 - 3 * S1[a]^2 * S2[a] + S2[a]^2 + 2 * p^3 * S1[a] * conjugate(S1[a]) - p^6;
print("mom_4^3(p) =",mom3,"and mom_4^3(p) =",mom4)

#compute trace a
a3 = -(mom+1 + p^2 + p^3)/p^4;
a4 = -(mom+1 + p^2 + p^3 + p^4 + 2 * p^6)/p^4;
print("a_4^3(p) =",a3, "and a_4^4(p) =",a4 )

# when N*k^2 is too big 
M = ModularSymbols(Gamma0(N),k,sign=1,base_ring=GF(13))
Hp = M.new_subspace().cuspidal_subspace().hecke_polynomial(p)
print (Hp(a44))





#moments of Kl5------------------------------------------------

#choose a prime
p=2;

#compute Kl_5(a;p) and Kl_5(a;p^2), then print s_1(a) and s_2(a)
import math


#compute Kl_5(a;p) as the list Klp
k = GF(p,'r')
Klp = [];
for a in range(1,p):
    klp = 0;
    for x in k:
        for y in k:
            for z in k:
                for w in k:
                    if x*y*z*w!=0:
                        u = x + y + z + w + a/(x * y * z * w);
                        u = (u.trace());
                        u = u.lift();
                        klp = klp + exp(2 * pi * I/p * u);
                        klp = n(klp);
    Klp.append(klp);
    
#compute Kl_5(a;p^2) as the list Klp2
k = GF(p^2,'r')
Klp2 = [];
for a in range(1,p):
    klp2 = 0;
    for x in k:
        for y in k:
            for z in k:
                for w in k:
                    if x*y*z*w!=0:
                        u = x + y + z + w + a/(x * y * z * w);
                        u = (u.trace());
                        u = u.lift();
                        klp2 = klp2 + exp(2 * pi * I/p * u);
                        klp2 = n(klp2);
    Klp2.append(klp2)


#compute s_1(a) and s_2(a) as lists S1 and S2
S1 = [];
S2 = [];
for a in range(1,p):
    b = a-1;
    s1 = Klp[b];
    s2 = Klp2[b];
    s2 = (s1^2 - s2)/2;
    S1.append(s1);
    S2.append(s2);
    
#compute moments m
mom=0;
for a in range(1,p):
    b = a - 1;
    mom = S1[b]^3 - 2 * S1[b]*S2[b] + p^2 * conjugate(S2[b]);
print("mom_5^3(p) =",mom)

#compute trace a
a = -(mom + 1 + p^2 + p^3 + p^4 + p^6)/p^5;
print("a_5^3(p) =",a)