Verifique la Hipótesis de Riemann para los primeros 1000 ceros.

Question

No sé si yo sería mejor publicar esto en la codificación de stackexchange o no, pero aquí va de todos modos:

Estoy esencialmente tratando de verificar la hipótesis de Riemann para los primeros 1000 ceros usando el de Euler-MacLaurin de expansión de $\zeta(s)$.

Pero aquí es donde mis problemas comienzan, en primer lugar mi código de python no parece estar funcionando, estoy esencialmente usando el método descrito en el Mike de Rubinstein - "Métodos de cálculo y los Experimentos en la Teoría Analítica de números" [https://arxiv.org/pdf/math/0412181v1.pdf] y aquí está el código que me han dado:

def CalculateRZ(s, tol):
    A = int(0)
    while (2*A - 1 < tol + 0.5 * math.log10(abs(s + 2*A - 2))):
        A = A + 1
    B = int(2 * A)
    N = int(math.ceil((10 / 2 * pi) * abs(s + B - 2)))

    sum1 = 0.0
    for n in range(1, N + 1):
        sum1 = sum1 + n ** (-s)

    ksum = 0.0
    for k in range(1, B + 1):
        ksum = ksum + Binomial(s + k - 2, k - 1) * (bernoulli(k) / k) * N ** (-s - k - 1)

    sum2 = (N ** (1 - s)) / (s - 1) + ksum

    EulerApx = sum1 + sum2

    return EulerApx

donde Binomial(n, k), y de bernoulli(k) son exactamente lo que parecen.

Esto parece perfectamente correcto para mí, sin embargo cuando me puse

x = 1/2 + 14.134725141734693790457251983562470270784257115699243175685567460149 * 1j

(es decir: el primer cero) y ejecutar

CalculateRZ(x, 100)

Puedo obtener como salida $\zeta(x) = -0.008282095160420333 + 0.009227302793957781j$ que puede ser algo cercano a cero, pero ciertamente no para de 100 dígitos de precisión y probablemente no lo suficiente exactitud para la Hipótesis de Riemann.

Además de esto, una vez que tengo este código de trabajo no estoy seguro de cómo ir sobre la comprobación de la Hipótesis de Riemann después.

Sería mejor crear una nueva función para $Z(t)$, entonces el uso de Newton Raphson para contar los ceros? Si es así ¿cómo puedo ir sobre la computación en la $Z(t)$? Debo hacer que sea como un retorno de función $\zeta(1/2+1t)e^{i\theta(t)}$ usando mi versión de zeta y una función que sólo hacia arriba calcula el $\theta(t)$? Donde he de poner mis puntos de inicio para la N-R método para asegurarse de que golpear a todos de cero una vez y sólo una vez?

Una vez que los ceros se cuentan entiendo que puedo integrar $\zeta$ sobre algunos rectángulo, y dividir ese número por $2\pi$ la cual me dará el número de singularidades que deben coincidir con el número de ceros encuentra el uso de Newton Raphson, por lo tanto muestra que todos los ceros de la mentira en la línea $Z(t)$, y por lo tanto la verificación de la RH. Pero, ¿cómo sería esta integral se calcula?

Sé que esto es una huuge pregunta, así que se agradece cualquier ayuda, incluso si no es una respuesta a la pregunta completa.

Muchas gracias por adelantado a todos los que se toman el tiempo de su día para que me ayude!

Answer 1

El mayor problema aquí es que el $1/2$ x = 1/2 + 14.1347 * 1j obtiene interpretarse como la división entera por lo que están tratando de calcular $\zeta(x)$ x = 14.1347 * 1j. Escribo como 0.5 lugar.

Entonces usted también debe arreglar el error mencionado en los comentarios $N^{-s-k-1} \to N^{-s-k+1}$. Ejecutando el código con estas revisiones da -1.17435267217e-14 - 3.36278987909e-16j que es bastante cercano a cero (en el doble de números de precisión, que es lo que estamos trabajando).

También tenga en cuenta que ya no eres el uso arbitrario de números de precisión en los cálculos, realmente no hay punto de correr con tol mayor que ~$20$.

La integridad de aquí es un completo código de trabajo

import math
from sympy import bernoulli

# Binomial coefficients for arbitrary n and integer k
def Binomial(n,k):
  res = 1.0
  for i in range(0,k):
    res *= (n-i)/(k-i)
  return res

def CalculateRZ(s, tol):
  A = int(0)
  while (2*A - 1 < tol + 0.5 * math.log10(abs(s + 2*A - 2))):
      A = A + 1
  B = int(2 * A)
  N = int(math.ceil((10 / (2 * math.pi)) * abs(s + B - 2)))

  sum1 = 0.0
  for n in range(1, N + 1):
      sum1 = sum1 + n ** (-s)

  ksum = 0.0
  for k in range(1, B + 1):
      b = Binomial(s + k - 2, k - 1)
      ksum = ksum + b * (bernoulli(k) / k) * N ** (-s - k + 1)

  sum2 = (N ** (1 - s)) / (s - 1) + ksum

  EulerApx = sum1 + sum2

  return complex(EulerApx)

s = 0.5 +  14.1347251417j
tol = 20 
print CalculateRZ(s, tol)