attachment:explicit_formula.sage of MSRItemp/ExplicitFormula

   1 import sys
   2 import mpmath
   3 
   4 psi = mpmath.digamma
   5 delta = log(2)/(2*pi)
   6 
   7 def von_mangoldt(n, R = RR):
   8     '''
   9     return Lambda(n), evaluated in the ring R
  10     '''
  11     if n == 1:
  12         return R(0)
  13     if not is_prime_power(n):
  14         return R(0)
  15     factors = factor(n)
  16     return R(factors[0][0]).log()
  17 
  18 def beurling_function(terms = 10, precision = 0):
  19     z = SR.var('z')
  20     B = 1 + 2 * (sin(pi * z)/pi)^2 * (1/z - sum( 1/(z+n)^2 for n in [1 .. terms] ) )
  21     return B
  22     
  23 def selberg_minus_function(terms = 10, precision = 0, delta = 1, alpha = -1, beta = 1):
  24     z = SR.var('z')
  25     B = beurling_function(terms = terms, precision = precision)
  26     if precision == 0:
  27         precision = 53
  28     Sminus = -(1/2)*B.subs(z = delta*(alpha - z)) - (1/2)*B.subs(z = delta*(z - beta))
  29     if precision == 53:
  30         return fast_callable(Sminus, domain=float, vars='z')
  31     else:
  32         return fast_callable(Sminus, domain=RealField(precision), vars='z')
  33 
  34 
  35 def selberg_plus_function(terms = 10, precision = 0, delta = 1, alpha = -1, beta = 1):
  36     z = SR.var('z')
  37     B = beurling_function(terms = terms, precision = precision)
  38     if precision == 0:
  39         precision = 53
  40     Sminus = (1/2)*(B.subs(z = delta*(z - alpha)) + B.subs(z = delta*(beta - z)))
  41     if precision == 53:
  42         return fast_callable(Sminus, domain=float, vars='z')
  43     else:
  44         return fast_callable(Sminus, domain=RealField(precision), vars='z')
  45         
  46 def selberg_minus_function_integer_width(delta = 1, alpha = -1, beta = 1, domain = float, symbolic = False):
  47     z = SR.var('z')
  48     terms = RDF(delta*(beta - alpha)).round() - 1
  49     Sminus = (sin(pi*delta * (alpha - z)))^2/pi^2 * (sum ( [1/(delta*(alpha - z) + n)^2 for n in srange(1, terms + 1)]) - 1/(delta*(z - beta)) - 1/(delta*(alpha - z)))
  50     if symbolic:
  51         return Sminus
  52     else:
  53         return fast_callable(Sminus, domain=domain, vars='z')
  54 
  55 def sinc_squared(delta = 1):
  56     def f(x):
  57         if x == 0:
  58             return RR(1)
  59         else:
  60             return (RR(delta * pi * x).sin()/RR(delta * pi * x))^2
  61     return f
  62 
  63 def sinc(delta = 1):
  64     def f(x):
  65         if x == 0:
  66             return RR(1)
  67         else:
  68             return (RR(delta * pi * x).sin()/RR(delta * pi * x))
  69     return f
  70 
  71 
  72 def tri(delta = 1):
  73     def f(x):
  74         return max( RR(1 - abs(x/delta)), RR(0))/RR(delta)
  75     return f
  76 
  77 def elliptic_curve_rank_bound(delta):
  78     """compute a function that will give an upper bound for the analytic rank of an elliptic curve of conductor N, assuming GRH.
  79     Uses the explicit formula with the test function sinc(delta x)."""
  80 
  81     f = sinc_squared(delta)
  82     f_hat = tri(delta)
  83 
  84     gamma_term1, error1 = mpmath.quad(lambda t : psi(1/4 + t * i/2).real * f(t), [-oo, oo], error = true)
  85     gamma_term1 = RR(gamma_term1.real)
  86     error1 = RR(error1)
  87 
  88     print error1
  89 
  90     if error1 > .1:
  91         print "warning: error in computing the first numerical intergral was probably too large."
  92 
  93 
  94     gamma_term2, error2 = mpmath.quad(lambda t : psi(1/4 + t * i/2 + 1/2) * f(t), [-oo, oo], error = true)
  95     gamma_term2 = RR(gamma_term2.real)
  96     error2 = RR(error2)
  97 
  98     print error2
  99 
 100     if error2 > .1:
 101         print "warning: error in computing the second numerical intergral was probably too large."
 102 
 103     N = SR.var('N')
 104 
 105     # the prime sum...
 106     # f_hat has support in -delta, delta, so our sum will be over the prime powers
 107     # p^alpha such that log(p^alpha)/2pi < delta, i.e. p < exp(2 pi delta)
 108 
 109     S = RR(0)
 110 
 111     print ceil(exp(2 * pi * delta))
 112     sys.stdout.flush()
 113 
 114     for n in sxrange(ceil(exp(2 * pi * delta))):
 115         x = von_mangoldt(n) * 2
 116 
 117         if x != 0:
 118             S = S + x/RR(n).sqrt() * f_hat(RR(n).log()/RR(2 * pi)) 
 119 
 120     S = S / RR(pi)
 121 
 122     bound = (gamma_term1/RR(2 * pi) + gamma_term2/RR(2 * pi) - f_hat(0) * RR(pi).log()/RR(pi) + f_hat(0)/RR(2 * pi) * log(N))/f(0) + S
 123     return N, bound
 124 
 125 def gamma_integrand_1(delta):
 126     f = sinc_squared(delta)
 127     return lambda t : psi(1/4 + t * i/2).real * f(t)