attachment:computing_modular_symbols_via_complex_integration.sage of days22/wuthrich

   1 def Gamma(E,r,an=None):
   2     alpha = r.numerator()
   3     beta = r.denominator()
   4     N = E.conductor()
   5     newcusp=0
   6     epsilon=var('epsilon')
   7     if gcd(beta,N) == 1:
   8         [c,d] = xgcd(N*alpha,beta)[1:3]
   9         c=-c
  10         d=-d
  11         Gamma = matrix(ZZ,[[-beta,alpha],[c*N,d]])
  12         assert Gamma.det() == 1
  13         #Gamma = Gamma0(N)([-beta,alpha,c*N,d])
  14         A = 0
  15         B = 1
  16         [C,D] = xgcd(-B,A)[1:3]
  17         M = matrix(ZZ,[[A,B],[C,D]])
  18         #M = SL2Z([A,C,B,D])
  19         M = M.transpose()
  20 
  21 
  22 
  23         Q = N/B
  24 
  25         [x,y] = [1,-1]
  26 
  27         #x = -x/A
  28         #y=-y
  29 
  30         wQ = matrix(ZZ,[[Q*x,y],[N,-A*Q]])
  31         wQ = (wQ*M)[0,0]*wQ
  32         assert wQ[1,1]==0
  33         h = (wQ*M)[1,1]
  34        
  35         ep_solutions_dict = solve(epsilon == 1/Q * imag(SL2Z(M.inverse()).acton(SL2Z(Gamma).acton(r+i*epsilon)) ), epsilon, solution_dict=True)
  36         ep_solutions = list()
  37         for j in range(0,len(ep_solutions_dict)):
  38             ep_solutions.append(ep_solutions_dict[j][epsilon])
  39         ep = max(ep_solutions)
  40         #print 'epsilon = ',ep, ' = ', round(ep,10)
  41         #assert ep <= 0, 'Error: epsilon negative!'
  42 
  43     else:
  44         [c,d] = xgcd(N*alpha,beta)[1:3]
  45         [a,b] = xgcd(d,-c*N)[1:3]
  46         newcusp = (a*alpha + b*beta)/(c*N*alpha + d*beta)
  47         Gamma = matrix(ZZ,[[a,b],[c*N,d]])
  48         #Gamma = Gamma0(N)([-beta,alpha,c*N,d])
  49         A = newcusp.numerator()
  50         B = newcusp.denominator()
  51         [C,D] = xgcd(-B,A)[1:3]
  52         M = matrix(ZZ,[[A,B],[C,D]])
  53         #M = SL2Z([A,C,B,D])
  54         M = M.transpose()
  55 
  56         Q = N/B
  57 
  58         [x,y] = xgcd(Q^2,N)[1:3]
  59 
  60         x = -x/A
  61         y=-y
  62 
  63         wQ = matrix(ZZ,[[Q*x,y],[N,-A*Q]])
  64         wQ = (wQ*M)[0,0]*wQ
  65         h = (wQ*M)[1,1]
  66    
  67         ep_solutions_dict = solve(epsilon == 1/Q * imag(SL2Z(M.inverse()).acton(SL2Z(Gamma).acton(r+i*epsilon)) ), epsilon, solution_dict=True)
  68         ep_solutions = list()
  69         for j in range(0,len(ep_solutions_dict)):
  70             ep_solutions.append(ep_solutions_dict[j][epsilon])
  71         ep = max(ep_solutions)
  72         #print 'epsilon = ',ep, ' = ', round(ep,10)
  73         if ep <= 0:
  74             print 'Error: epsilon negative!'
  75 
  76     #manually fix ep
  77     #ep = .05
  78     #
  79     T = SL2Z(M.inverse()).acton(SL2Z(Gamma).acton(r+I*ep))
  80     #delta = E.period_lattice().basis()[0] / (2*E.torsion_order() )
  81     #print 'tau1 = ', round(r,10) + I*round(ep,10)
  82     #print 'tau2 = ', round(real((T+h)/Q),10) + I*round(imag((T+h)/Q),10)
  83     integral1 = CompApprox(E, r + I*ep,an)
  84     integral2 = CompApprox(E, (T + h)/Q,an)
  85    
  86     eps = dict([d,prod(E.root_number(p) for p in d.prime_divisors() )] for
  87 d in E.conductor().divisors())
  88     AL_eigenvalue = eps[Q]
  89    
  90     integral = integral1 - AL_eigenvalue*(round(real(integral2),10) + I * round(imag(integral2),10))
  91     m = E.torsion_order() * real(integral) / E.period_lattice().basis()[0]
  92     integral_exact = m.round() / E.torsion_order()
  93     #print integral / E.period_lattice().basis()[0]
  94 
  95     return integral_exact
  96 
  97 def CompApprox(E, tau, an=None):
  98     delta = E.period_lattice().basis()[0] / (3 * E.torsion_order() )
  99     #delta = 10^-6
 100     #print 'delta = ',delta
 101     epsilon_local = imag_part(tau)
 102     M = max(ceil((-1/(2*pi*epsilon_local))*ln(delta*9999/10000*(1-exp(-2*pi*epsilon_local)))),100) #optional lower bound
 103     #print 'M: ',M
 104     if M>10^9:
 105         print 'Warning: M large'
 106     b=max(ceil(log(2*pi*M,2) - log(log(1+delta/10000 / M),2)), 25)
 107     CC=ComplexField(b)
 108     tau = CC(tau)
 109     #print 'bits: ',CC(b)
 110     if an==None or len(an) < M+1:
 111         an = E.anlist(M)
 112     #print 'length of an: ',len(an)
 113     return sum(   (an[n]/n)*exp(2*CC.pi()*CC(I)*n*tau)  for n in range(1, M+1) )