Theta Series Automated.sagews

var_string = ""
deg = 4
for i in range(deg): 
    var_string += "a"
    var_string += str(i)
    var_string += ", "
var_string += "a"
var_string += str(deg)

a = var(var_string)
RRR.<x,y> = SR['x,y']

f = sum([a[i] * x^(deg-i) * y^i for i in range(deg+1)])

A, B, C = var('A, B, C')
QF = A*x^2 + B*x*y + C*y^2
Lap_f = C*derivative(f,x,2) - B * derivative(f,x,y) + A * derivative(f,y,2)

mat_var = matrix([[Lap_f.coefficients()[j].coefficient(a[i]) for i in range(deg+1)] for j in range(3)])


#Choose values for the quadratic form here. N is the discriminant.
AA = 1
BB = 1
CC = 3
N = (BB^2 - 4*AA*CC).abs()


#Redefining values. Before, we computed using symbolic expressions. Now, will work in rational polynomials in x and y.
mat = mat_var.subs(A==AA, B==BB, C==CC)
A = AA
B = BB
C = CC
P.<x, y> = QQ['x,y']
QF = A*x^2 + B*x*y + C*y^2


t1 = 2
t2 = 3
harm_bas = mat.right_kernel().basis()


a = t1*harm_bas[0] + t2*harm_bas[1]
f = sum([a[i] * x^(deg-i) * y^i for i in range(deg+1)])


R.<q> = QQ[['q']]
bound = ModularForms(Gamma1(N), 1+deg).sturm_bound() + 2

F = RealField(10)
lim = bound



cand_mod = sum([sum([QQ(f(x=a,y=b)) * q^(A*a^2 + B*a*b + C*b^2) for a in range(-lim,lim)]) for b in range(-lim,lim)]) + O(q^bound)
cand_mod_coeffs = cand_mod.padded_list(bound)

M = ModularForms(Gamma1(N), 1+deg, prec=bound)
bas = M.basis()

mod_bas_list = [bas[i].padded_list(bound) for i in range(len(bas))]
mod_bas_mat = matrix(mod_bas_list)



is_mod_Bcoeffs = mod_bas_mat.solve_left(vector(cand_mod_coeffs))
is_mod = sum([is_mod_Bcoeffs[i] * bas[i] for i in range(len(bas))])
is_mod

print "\n"

cand_mod

4*q + 28*q^3 + 64*q^4 - 196*q^5 - 128*q^9 + 484*q^11 + 448*q^12 - 1372*q^15 + 1024*q^16 - 3136*q^20 + 668*q^23 + 7104*q^25 - 3164*q^27 - 2212*q^31 + 3388*q^33 - 2048*q^36 - 8452*q^37 + 7744*q^44 + 6272*q^45 - 7672*q^47 + 7168*q^48 + 9604*q^49 + O(q^53) 4*q + 28*q^3 + 64*q^4 - 196*q^5 - 128*q^9 + 484*q^11 + 448*q^12 - 1372*q^15 + 1024*q^16 - 3136*q^20 + 668*q^23 + 7104*q^25 - 3164*q^27 - 2212*q^31 + 3388*q^33 - 2048*q^36 - 8452*q^37 + 7744*q^44 + 6272*q^45 - 7672*q^47 + 7168*q^48 + 9604*q^49 + O(q^53)