5 segundos para encontrar pastel

11

Pi veces e (o Pie si te gusta la notación ambigua) a 100 decimales es:

8.5397342226735670654635508695465744950348885357651149618796011301792286111573308075725638697104739439...

( OIES A019609 ) ( argumento para una posible irracionalidad )

Su tarea es escribir un programa que tome un número entero positivo N y produzca Pi * e truncado a N decimales. por ejemplo, si N = 2, entonces la salida debería ser 8.53.

Este es un problema de optimización, por lo que ganará el envío que puede proporcionar la salida correcta para el valor más alto de N.

Para garantizar que todas las presentaciones se juzguen utilizando la misma potencia informática, su código debe ejecutarse en ideone , utilizando cualquier lenguaje que admitan. De acuerdo con las preguntas frecuentes de ideone , hay un límite de tiempo de ejecución de 5 segundos para los usuarios que no han iniciado sesión. Este límite de 5 segundos es el que debe usar, no el límite de 15 segundos para los usuarios registrados. (Consulte las preguntas frecuentes para conocer otros límites, como la memoria, el tamaño del código, etc.)

Específicamente, cualquier persona que no haya iniciado sesión en ideone debería poder ejecutar su programa en ideone para todos los valores de N desde 1 hasta un máximo de Nmax, y ver la salida correcta casi todo el tiempo . sin ninguna Time limit exceededo Memory limit exceeded, etc. errores. La presentación con el mayor Nmax gana.

(No importa si el tiempo real que se tarda es una pizca de más de 5 segundos, siempre y cuando ideone no dé errores. " Casi todo el tiempo " se define como más del 95% del tiempo para cualquier N. en particular)

Detalles

  • Puede usar cualquier método matemático que desee para calcular Pi * e, pero no puede codificar la salida más allá de la primera docena de dígitos de Pi, e o Pi * e .
    • Su programa debería poder funcionar para cualquier N, dados los recursos ilimitados.
    • Puede usar constantes Pi o e integradas si su idioma las tiene.
  • No puede acceder a sitios web o recursos externos a su código (si ideone lo permite).
  • Más allá de la codificación y el acceso a recursos externos, cualquier cosa que ideone permita es casi seguro que está bien.
  • Su entrada y salida debe (obviamente) funcionar con lo que ideone proporciona para E / S (stdin / stdout solo parece). Está bien si necesita comillas alrededor de la entrada N o la salida es algo así ans = ..., etc.
  • Incluya un enlace a un fragmento de ideone de su código con su Nmax como entrada.
  • Si sucede que hay un empate (solo es probable si las presentaciones múltiples alcanzan el límite de caracteres de salida de 64kB), la respuesta de votos más alta gana.
Pasatiempos de Calvin
fuente
3
Mmm ... pastel ambiguo.
Dennis
Esto puede ser fácilmente un código de golf y preferiría ser más divertido si lo es.
Optimizador
2
@Optimizer Podría ser code-golf, pero sería muy similar a cualquier otro golf de código de generación de dígitos. Quería probar un concurso basado en el tiempo que pudiera verificarse en línea. (Aunque un problema computacionalmente más complejo podría haber sido mejor.)
Calvin's Hobbies el
si este es el código de golf APL probablemente ganaría (menos la parte de precisión arbitraria)
TwiNight
1
Sospecho que estos programas estarán completamente obligados a escribir los dígitos en stdout. Cinco segundos es mucho tiempo para algo así como y-cruncher .
Será el

Respuestas:

12

Python - 65535

http://ideone.com/knKRhn

from math import exp, log

def divnr(p, q):
  """
    Integer division p/q using Newton-Raphson Division.
    Assumes p > q > 0.
  """

  sp = p.bit_length()-1
  sq = q.bit_length()-1
  sr = sp - sq + 1

  s = []
  t = sr
  while t > 15:
    s = [t] + s
    t = (t>>1) + 1
  # Base-case division
  r = (1 << (t<<1)) / (q >> sq-t)

  for u in s:
    r = (r << u-t+1) - (r*r * (q >> sq-u) >> (t<<1))
    t = u
  return (r * (p >> sq)) >> sr

def pibs(a, b):
  if a == b:
    if a == 0:
      return (1, 1, 1123)
    p = a*(a*(32*a-48)+22)-3
    q = a*a*a*24893568
    t = 21460*a+1123
    return (p, -q, p*t)
  m = (a+b) >> 1
  p1, q1, t1 = pibs(a, m)
  p2, q2, t2 = pibs(m+1, b)
  return (p1*p2, q1*q2, q2*t1 + p1*t2)

def ebs(a, b):
  if a == b:
    if a == 0:
      return (1, 1)
    return (1, a)
  m = (a+b) >> 1
  p1, q1 = ebs(a, m)
  p2, q2 = ebs(m+1, b)
  return (p1*q2+p2, q1*q2)

if __name__ == '__main__':
  n = input()

  pi_terms = int(n*0.16975227728583067)

  # 10^n == e^p
  p = n*2.3025850929940457

  # Lambert W_0(p/e) a la Newton
  k = log(p) - 1
  w = k - (k-1)/(k+1)
  while k > w:
    k = w
    w -= (k - p*exp(-k-1))/(k+1)

  # InverseGamma(e^p) approximation
  e_terms = int(p / w)

  pp, pq, pt = pibs(0, pi_terms)
  ep, eq = ebs(0, e_terms)

  z = 10**n
  p = 3528*z*ep*abs(pq)
  q = eq*abs(pt)

  pie = divnr(p, q)
  print pie,

Ideone no parece haberse gmpy2instalado, lo cual es lamentable por al menos dos razones. Uno, porque haría el cálculo mucho más rápido, y dos, porque hace que cualquier fórmula que requiera una raíz cuadrada de precisión arbitraria sea poco práctica.

La fórmula que uso para π fue listada por Ramanujan como Fórmula (39):

que converge a una tasa de ~ 5.89 dígitos por término. Que yo sepa, esta es la serie convergente más rápida de su tipo que no requiere la evaluación de una raíz cuadrada de precisión arbitraria. La fórmula (44) en el mismo documento (tasa de convergencia ~ 7.98 dígitos por término) se conoce con mayor frecuencia como la fórmula de Ramanujan.

La fórmula que uso para e es la suma de factoriales inversos. El número de términos requeridos se calcula como Γ -1 ( 10 n ), usando una aproximación que encontré en mathoverflow . El componente Lambert W 0 se encuentra utilizando el Método de Newton.

El cálculo de cada una de estas sumas se realiza a través de la evaluación rápida de la función E (más generalmente conocida como división binaria), ideada originalmente por Karatsuba. El método reduce una suma de n términos a un solo valor racional p / q . Estos dos valores se multiplican para producir el resultado final.

Actualización: el
perfil reveló que más de la mitad del tiempo necesario para el cálculo se gastó en la división final. Solo se necesitan los bits más altos de log 2 (10 n ) de q para obtener una precisión completa, por lo que recorto algunos de antemano. El código ahora llena el búfer de salida de Ideone en 3.33s .

Actualización 2:
Dado que este es un desafío de , decidí escribir mi propia rutina de división para combatir la lentitud de CPython. La implementación de lo divnranterior usa la División Newton-Raphson . La idea general es calcular d = 1 / q · 2 n usando el Método de Newton, donde n es el número de bits que requiere el resultado, y calcular el resultado como p · d >> n . El tiempo de ejecución ahora es 2.87s , y esto es sin cortar bits antes del cálculo; No es necesario para este método.

primo
fuente
4

PARI / GP: 33000

Este es básicamente el programa dado en OEIS , modificado para tomar la entrada y formatear la salida correctamente. Debería servir como línea de base para vencer, si nada más.

Yo supongo que esto es exacto. Lo comprobé a 100 y 20k contra OEIS, y coincidió con ambos. Es bastante difícil encontrar más dígitos en línea para verificar.

Para 33,000 toma alrededor de 4.5s, por lo que probablemente podría ser golpeado un poco. Me cansé de jugar con la entrada y el bucle de envío / compilación / ejecución lento de ideone.

{ 
    m=eval(input());
    default(realprecision, m+80); 
    x=Pi*exp(1);
    s="8.";
    d=floor(x);
    x=(x-d)*10;
    for (n=1, m, d=floor(x); 
         x=(x-d)*10; 
         s=Str(s,d));
    print(s);
}

Enlace de Ideone.com

Geobits
fuente
Tus dígitos coinciden con los míos, así que voy a arriesgarme y decir que probablemente sean correctos.
primo
Este programa pasa esencialmente todo su tiempo en el bucle, generando dígitos uno por uno. Si solo lo toma Str(floor(frac(x)*10^m), va cientos / miles de veces más rápido.
Charles
2

Python 3

Como el construido en pi y e no tiene suficientes dígitos, decidí calcular el mío.

import decimal
import math
decimal.getcontext().prec=1000000
decimal=decimal.Decimal;b=2500
print(str(sum([decimal(1/math.factorial(x)) for x in range(b)])*sum([decimal(1/16**i*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6))) for i in range(b)]))[0:int(input())+2])

IDEOne enlace

Salida para STDIN = 1000:

8.5397342226735669504281233688422467204743749305568824722710929852470173635361001388261308713809518841081669216573834376992066672804294594807026609638293539437286935503772101801433821053915371716284188665787967232464763808892618434263301810056154560438283877633957941572944822034479453916753507796910068912594560500836608215235605783723340714760960119319145912948480279651779184356994356172418603464628747082162475871780202868607325544781551065680583616058471475977367814338295574582450942453416002008665325253385672668994300796223139976640645190237481531851902147391807396201201799703915343423499008135819239684881566321559967077443367982975103648727755579256820566722752546407521965713336095320920822985129589997143740696972018563360331663471959214120971348584257396673542429063767170337770469161945592685537660073097456725716654388703941509676413429681372333615691533682226329180996924321063261666235129175134250645330301407536588271020457172050227357541822742441070313522061438812060477519238440079
Decaimiento Beta
fuente
Nmax es el valor de entrada más grande que puede darle a su programa antes de que ideone ya no lo ejecute.
Aficiones de Calvin
1
@ Calvin'sHobbies Creo que nmax es arbitrariamente grande, sin embargo ...
Beta Decay
1
ideone no te da una potencia informática infinita. ¿Cuál es el mayor valor de entrada que su programa puede ejecutar en ideone? (Aunque, de hecho, su programa no sigue la should be able to work for any N, given unlimited resourcesregla. La mayor parte de los resultados son ceros en torno a N = 10000.)
Calvin's Hobbies
Eso no es python3: NameError: name 'xrange' not defined.
Bakuriu
2

Scala - 1790

IDEOne en http://ideone.com/A2CIto .

Usamos la fórmula de Wetherfield para π (y el código de fórmula de Machin toscamente portado desde aquí ). Calculamos e usando la serie de potencia ordinaria.

object Main extends App {
  import java.math.{BigDecimal => JDecimal}
  import java.math.RoundingMode._
  import scala.concurrent.Future
  import scala.concurrent.Await
  import scala.concurrent.ExecutionContext.Implicits._
  import scala.concurrent.duration._
  val precision = 1800

  def acotPrecision(numDigits: Int)(x: BigDecimal) = {
    val x1 = x.underlying
    val two = JDecimal.valueOf(2)
    val xSquared = x1 pow 2
    val unity = JDecimal.ONE.setScale(numDigits, HALF_EVEN)
    var sum = unity.divide(x1, HALF_EVEN)
    var xpower = new JDecimal(sum.toString)
    var term = unity

    var add = false

    var n = JDecimal.valueOf(3).setScale(numDigits)
    while (term.setScale(numDigits, HALF_EVEN).compareTo(JDecimal.ZERO) != 0) {
      xpower = xpower.divide(xSquared, HALF_EVEN)
      term = xpower.divide(n, HALF_EVEN)
      sum = if (add) sum add term else sum subtract term
      add = !add
      n = n add two
    }
    sum
  }

  def ePrecision(numDigits: Int) = {
    val zero = JDecimal.ZERO
    var sum = zero
    var term = JDecimal.ONE.setScale(numDigits, HALF_EVEN)
    var n = JDecimal.ONE.setScale(numDigits, HALF_EVEN)
    while(term.setScale(numDigits, HALF_EVEN).compareTo(zero) != 0) {
      sum = sum add term
      term = term.divide(n, HALF_EVEN)
      n = n add JDecimal.ONE
    }
    sum
  }

  val acot = acotPrecision(precision) _

  def d(x: Int) = JDecimal.valueOf(x)

  def piFuture = Future(d(4) multiply (
    (d(83) multiply acot(107)) add (d(17) multiply acot(1710)) subtract (d(22) multiply acot(103697))
    subtract (d(24) multiply acot(2513489)) subtract (d(44) multiply acot(18280007883L))
   add (d(12) multiply acot(7939642926390344818L))
   add (d(22) multiply acot(BigDecimal("3054211727257704725384731479018")))
  ))

  def eFuture = Future(ePrecision(precision))

  Await.result(
    for (pi <- piFuture;
         e <- eFuture) yield println((pi multiply e).setScale(precision - 10, DOWN))
  , 5 seconds) 
}
James_pic
fuente