Evaluar la función Zeta de Riemann en un número complejo

11

Introducción

Encontré esta pregunta que estaba cerrada porque no estaba clara, pero fue una buena idea. Haré todo lo posible para convertir esto en un desafío claro.

La función Riemann Zeta es una función especial que se define como la continuación analítica de

ingrese la descripción de la imagen aquí

al plano complejo Hay muchas fórmulas equivalentes que lo hacen interesante para el golf de código.

Desafío

Escriba un programa que tome 2 flotantes como entrada (la parte real e imaginaria de un número complejo) y evalúe la función Riemann Zeta en ese punto.

Reglas

  • Entrada y salida a través de la consola O función de entrada y valor de retorno
  • Los números complejos integrados no están permitidos, use flotadores (número, doble, ...)
  • Sin funciones matemáticas, excepto funciones + - * / pow logtrigonométricas reales (si desea integrar, use la función gamma, ... debe incluir esta definición de funciones en el código)
  • Entrada: 2 flotadores
  • Salida: 2 flotadores
  • Su código debe contener un valor que proporcione una precisión teóricamente arbitraria cuando se hace arbitrario grande / pequeño
  • El comportamiento en la entrada 1 no es importante (este es el único polo de esta función)

¡El código más corto en bytes gana!

Ejemplo de entrada y salida

Entrada:

2, 0

Salida:

1.6449340668482266, 0

Entrada:

1, 1

Salida:

0.5821580597520037, -0.9268485643308071

Entrada:

-1, 0

Salida:

-0.08333333333333559, 0

Jens Renders
fuente
1
¿Cuál es la precisión de salida requerida? No estoy seguro de entender que su código debe contener un valor que proporcione una precisión teóricamente arbitraria cuando se hace arbitrario grande / pequeño . ¿Te refieres a un valor máximo de bucle que cuando se incrementa sin límite proporciona una mayor precisión? ¿Puede ese valor estar codificado?
Luis Mendo
@DonMuesli Esto significa que la precisión depende de un parámetro, digamos N, al que le puede dar el valor que desee, pero para cualquier precisión dada, puede hacer que N sea lo suficientemente pequeño o grande para lograr esa precisión. La palabra teóricamente existe porque no debe preocuparse por la precisión limitada de la máquina o el lenguaje.
Jens Renders
Para aclarar aún más N: es suficiente que para cualquier enlace epse entrada xexista un Nque se calcule zeta(x)dentro eps; o debe existir una Nque dependa solo epsy garantice que para cualquier x(o tal vez para algo xmás que una función dada epsdesde el polo) logre el límite; o puede Ndepender de x, pero las respuestas deben explicar cómo calcular Ndado xy eps? (Mi teoría analítica de números no es demasiado, pero sospecho que las opciones 2 y 3 van a estar más allá de todos menos uno o dos pósters regulares).
Peter Taylor
@PeterTaylor N lo suficientemente grande: para cualquiera xy para cualquiera epsdebe existir un Ptal que para toda N>Pla salida esté más cerca que epsel valor exacto. ¿Está esto claro? ¿Necesito aclararlo para el caso con N lo suficientemente pequeño?
Jens Renders
No, eso está suficientemente claro.
Peter Taylor

Respuestas:

8

Python - 385

Esta es una implementación sencilla de la Ecuación 21 de http://mathworld.wolfram.com/RiemannZetaFunction.html. Utiliza la convención de Python para argumentos opcionales; si desea especificar una precisión, puede pasar un tercer argumento a la función; de lo contrario, utiliza 1e-24 de forma predeterminada.

import numpy as N
def z(r,i,E=1e-24):
 R=0;I=0;n=0;
 while(True):
  a=0;b=0;m=2**(-n-1)
  for k in range(0,n+1):
   M=(-1)**k*N.product([x/(x-(n-k))for x in range(n-k+1,n+1)]);A=(k+1)**-r;t=-i*N.log(k+1);a+=M*A*N.cos(t);b+=M*A*N.sin(t)
  a*=m;b*=m;R+=a;I+=b;n+=1
  if a*a+b*b<E:break
 A=2**(1-r);t=-i*N.log(2);a=1-A*N.cos(t);b=-A*N.sin(t);d=a*a+b*b;a=a/d;b=-b/d
 print(R*a-I*b,R*b+I*a)
RT
fuente
z(2,0)da un valor incorrecto, debe ser pi ^ 2/6.
GuillaumeDufay
4

Python 3 , 303 297 bytes

Esta respuesta se basa en la respuesta de Python de RT con varias modificaciones:

  • En primer lugar, Binomial(n, k)se define como p = p * (n-k) / (k+1)que cambia Binomial(n,k)a Binomial(n,k+1)con cada pasada del bucle.
  • En segundo lugar, se (-1)**k * Binomial(n,k)convirtió en lo p = p * (k-n) / (k+1)que voltea el signo en cada paso del ciclo for.
  • Tercero, el whilebucle ha sido cambiado para verificar inmediatamente si a*a + b*b < E.
  • En cuarto lugar, el bit a bit no operador ~se utiliza en varios lugares en los que ayudarían a jugar al golf, usando identidades tales como -n-1 == ~n, n+1 == -~n, y n-1 == ~-n.

Se hicieron varias otras modificaciones pequeñas para mejorar el golf, como colocar el forbucle en una línea y la llamada printen una línea con el código anterior.

Sugerencias de golf bienvenidas. Pruébalo en línea!

Editar: -6 bytes de una serie de pequeños cambios.

import math as N
def z(r,i,E=1e-40):
 R=I=n=0;a=b=1
 while a*a+b*b>E:
  a=b=0;p=1;m=2**~n
  for k in range(1,n+2):M=p/k**r;p*=(k-1-n)/k;t=-i*N.log(k);a+=M*N.cos(t);b+=M*N.sin(t)
  a*=m;b*=m;R+=a;I+=b;n+=1
 A=2**-~-r;t=-i*N.log(2);x=1-A*N.cos(t);y=A*N.sin(t);d=x*x+y*y;return(R*x-I*y)/d,(R*y+I*x)/d
Sherlock9
fuente
1

Axiom, 413 315 292 bytes

p(n,a,b)==(x:=log(n^b);y:=n^a;[y*cos(x),y*sin(x)]);z(a,b)==(r:=[0.,0.];e:=10^-digits();t:=p(2,1-a,-b);y:=(1-t.1)^2+t.2^2;y=0=>[];m:=(1-t.1)/y;q:=t.2/y;n:=0;repeat(w:=2^(-n-1);abs(w)<e=>break;r:=r+w*reduce(+,[(-1)^k*binomial(n,k)*p(k+1,-a,-b) for k in 0..n]);n:=n+1);[r.1*m-q*r.2,m*r.2+r.1*q])

Esto también implementaría la ecuación 21 de http://mathworld.wolfram.com/RiemannZetaFunction.html. Lo anterior debería ser el que interpretó la función Axioma z (a, b) aquí 16 veces más lenta que la siguiente función Zeta (a, b) [ ese debería ser el que se compiló] todos sin golf y comentados [1 segundo para Zeta () contra 16 segundos para z () para un valor de 20 dígitos después del punto flotante]. Para la pregunta de dígitos, uno elegiría la precisión llamando a dígitos (); función, por ejemplo dígitos (10); z (1,1) debe imprimir 10 dígitos después del punto, pero dígitos (50); z (1,1) debe imprimir 50 dígitos después del punto.

-- elevImm(n,a,b)=n^(a+i*b)=r+i*v=[r,v]
elevImm(n:INT,a:Float,b:Float):Vector Float==(x:=log(n^b);y:=n^a;[y*cos(x),y*sin(x)]::Vector Float);

--                      +oo               n
--                      ---              ---
--             1        \       1        \            n 
--zeta(s)= ---------- * /     ------  *  /    (-1)^k(   )(k+1)^(-s)
--          1-2^(1-s)   ---n  2^(n+1)    ---k         k  
--                       0                0


Zeta(a:Float,b:Float):List Float==
  r:Vector Float:=[0.,0.]; e:=10^-digits()

  -- 1/(1-2^(1-s))=1/(1-x-i*y)=(1-x+iy)/((1-x)^2+y^2)=(1-x)/((1-x)^2+y^2)+i*y/((1-x)^2+y^2)    

  t:=elevImm(2,1-a,-b);
  y:=(1-t.1)^2+t.2^2;
  y=0=>[] 
  m:=(1-t.1)/y; 
  q:=t.2/y
  n:=0
  repeat
     w:=2^(-n-1)
     abs(w)<e=>break  --- this always terminate because n increase
     r:=r+w*reduce(+,[(-1)^k*binomial(n,k)*elevImm(k+1,-a,-b) for k in 0..n])
     n:=n+1
  -- (m+iq)(r1+ir2)=(m*r1-q*r2)+i(m*r2+q*r1)
  [r.1*m-q*r.2,m*r.2+r.1*q]

this is one test for the z(a,b) function above:

(10) -> z(2,0)
   (10)  [1.6449340668 482264365,0.0]
                                              Type: List Expression Float
(11) -> z(1,1)
   (11)  [0.5821580597 520036482,- 0.9268485643 3080707654]
                                              Type: List Expression Float
(12) -> z(-1,0)
   (12)  [- 0.0833333333 3333333333 3,0.0]
                                              Type: List Expression Float
(13) -> z(1,0)
   (13)  []
RosLuP
fuente