attachment:Teich-twist.sage of days22/wuthrich

   1 def Teich_twist(E = EllipticCurve('11a1'),p=5,pt=2,n=5,prec=5):
   2     Lserie = E.padic_lseries(p)
   3     alp = Lserie.alpha(prec)
   4     MSP = E.modular_symbol()
   5     MSM = E.modular_symbol(sign=-1)
   6 
   7     bounds = Lserie._prec_bounds(n,prec)
   8     padic_prec = max(bounds[1:]) + 5
   9     res_series_prec=min(p**(n-1),prec)
  10 
  11     K = QQ
  12     gamma = K(1+p)
  13     R = PowerSeriesRing(K,'T',res_series_prec)
  14     T = R(R.gen(),res_series_prec)
  15     L = R(0) #going to be the L-series eventually
  16     one_plus_T_factor = R(1)
  17     gamma_power = K(1)
  18     teich = Lserie.teichmuller(padic_prec)
  19     p_power = p**(n-1)
  20 
  21     for j in range(p_power):
  22         s = K(0)
  23         for a in range(1,p):
  24             b = teich[a]*gamma_power
  25             if pt & 2 == 0:
  26                 s += teich[a]^pt * ( (1/alp)^n * MSP((b)/(p^n) ) - (1/alp)^(n+1) * MSP((b)/(p^(n-1))) )
  27             else:
  28                 s += teich[a]^pt * ( (1/alp)^n * MSM((b)/(p^n) ) - (1/alp)^(n+1) * MSM((b)/(p^(n-1))) )
  29         L += s * one_plus_T_factor
  30         one_plus_T_factor *= (1+T)
  31         gamma_power *= gamma
  32     print 'Ellptic Curve : %s' % E.cremona_label()
  33     print 'prime : %s ' % p
  34     print 'the power of Teichmuller : %s ' % pt
  35     print 'n : %s' % n
  36     print 'precision : %s' % prec
  37 
  38     Kp = Qp(p,padic_prec,print_mode='series')
  39     Rp = PowerSeriesRing(Kp,'T',res_series_prec)
  40     Lprime = Rp(L,res_series_prec)
  41     aj = Lprime.list()
  42     if len(aj) > 0:
  43            aj = [aj[0].add_bigoh(padic_prec-2)] + [aj[j].add_bigoh(bounds[j]) for j in range(1,len(aj))]
  44     L = Rp(aj,res_series_prec ) 
  45     if pt & 2 == 0:
  46         L = L / Lserie._E.real_components()
  47     else:
  48         L = L
  49 
  50     print 'The p-adic L-function of a given elliptic curve twisted by the given power of the Teichmuller character is'
  51     print L