attachment:approximation_of_integral.sage of days22/wuthrich

   1 def CompApprox(E, tau, l = 0, an=None):
   2    
   3    # delta determines how close we would like to get to our answer
   4    delta = E.period_lattice().basis()[0] / (2 * E.torsion_order() )
   5    print 'delta: ',delta
   6    if l == 0:
   7        epsilon = imag_part(tau)
   8    else:
   9        epsilon = l * imag_part(tau)
  10       
  11    # How far we need to sum given delta (and epsilon)
  12    M = ceil((-1/(2*pi*epsilon))*ln(delta*9999/10000*(1-exp(-2*pi*epsilon))))
  13    print 'M: ',M
  14    if M>10^9:
  15        print 'Warning: M large'
  16       
  17    # Given delta and M, let's determine the bits of precision necessary
  18    # to carry out the computations to our wanted degree of accuracy.
  19   
  20    b=max(ceil(log(2*pi*M,2) - log(log(1+delta/10000 / M),2)), 20)
  21    CC=ComplexField(b)
  22    tau = CC(tau)
  23    print 'bits: ',CC(b)
  24    if an==None or len(an) < M+1:
  25        an = E.anlist(M)
  26    print 'length of an: ',len(an)
  27 
  28    #Speeding things up for large M
  29    partialSum = list()
  30    partialSum.append((an[0]/1)*exp(2*CC.pi()*CC(I)*1*tau))
  31 
  32    n=1
  33    while n < M+1:
  34        partialSum.append((an[n]/n)*exp(2*CC.pi()*CC(I)*n*tau) + partialSum[n-1])
  35        n+=1  
  36   
  37    #Give us our result   
  38    return partialSum[n-1]