E = EllipticCurve('11a1')
n=5; prec=5
Lserie = E.padic_lseries(5)
alp = Lserie.alpha(prec)
MSP = E.modular_symbol()
MSM = E.modular_symbol(sign=-1)
p = Lserie._p
bounds = Lserie._prec_bounds(n,prec)
padic_prec = max(bounds[1:]) + 5
res_series_prec=min(p**(n-1),prec)
K = QQ
gamma = K(1+p)
R = PowerSeriesRing(K,'T',res_series_prec)
T = R(R.gen(),res_series_prec)
L = R(0)
one_plus_T_factor = R(1)
gamma_power = K(1)
teich = Lserie.teichmuller(padic_prec)
pt = 2
p_power = p**(n-1)
for j in range(p_power):
s = K(0)
for a in range(1,p):
b = teich[a]*gamma_power
if pt & 2 == 0:
s += teich[a]^pt * ( (1/alp)^n * MSP((b)/(p^n) ) - (1/alp)^(n+1) * MSP((b)/(p^(n-1))) )
else:
s += teich[a]^pt * ( (1/alp)^n * MSM((b)/(p^n) ) - (1/alp)^(n+1) * MSM((b)/(p^(n-1))) )
L += s * one_plus_T_factor
one_plus_T_factor *= (1+T)
gamma_power *= gamma
print 'L before p-adic approx is %s' % L
K = Qp(p,padic_prec,print_mode='series')
R = PowerSeriesRing(K,'T',res_series_prec)
Lprime = R(L,res_series_prec)
print 'Lprime is %s' % Lprime
aj = Lprime.list()
if len(aj) > 0:
aj = [aj[0].add_bigoh(padic_prec-2)] + [aj[j].add_bigoh(bounds[j]) for j in range(1,len(aj))]
L = R(aj,res_series_prec )
if pt & 2 == 0:
L = L / Lserie._E.real_components()
else:
L = L
print 'The 5-adic L-function of 11A1 twisted by the (pt)-th power of the Teichmuller character is %s ' % L