5 segundos para encontrar a torta

11

Pi vezes e (ou Torta, se você gosta de notação ambígua) com 100 casas decimais é:

8.5397342226735670654635508695465744950348885357651149618796011301792286111573308075725638697104739439...

( OIES A019609 ) ( argumento para possível irracionalidade )

Sua tarefa é escrever um programa que receba um número inteiro positivo N e produza Pi * e truncado para N casas decimais. por exemplo, se N = 2, a saída deve ser 8.53.

Esse é um problema de otimização, portanto a submissão que pode fornecer a saída correta para o valor mais alto de N será vencedora.

Para garantir que todos os envios sejam julgados usando o mesmo poder de computação, seu código deve ser executado em ideone , usando qualquer idioma compatível. De acordo com o ideone faq , existe um limite de tempo de execução de 5 segundos para usuários não conectados. Esse limite de 5 segundos é o que você deve usar, não o limite de 15 segundos para usuários conectados. (Consulte o FAQ para outros limites, como memória, tamanho do código etc.)

Especificamente, qualquer pessoa que não esteja logada no ideone deve poder executar seu programa no ideone para todos os valores de N de 1 a um máximo de Nmax e ver a saída correta quase o tempo todo . sem erros Time limit exceededou Memory limit exceededetc. A finalização com o maior Nmax vence.

(Se o tempo real gasto é menor do que 5 segundos não importa, desde que a ideona não cometa erros. " Quase o tempo todo " é definido como mais de 95% do tempo para qualquer N. em particular)

Detalhes

  • Você pode usar qualquer método matemático que desejar para calcular Pi * e, mas não pode codificar a saída além das primeiras dezenas de dígitos de Pi, e ou Pi * e .
    • Seu programa deve poder trabalhar para qualquer N, com recursos ilimitados.
    • Você pode usar constantes Pi ou e incorporadas se o seu idioma as possuir.
  • Você não pode acessar sites ou recursos externos ao seu código (se alguém permitir isso).
  • Além da codificação e do acesso a recursos externos, tudo o que a ideona permite é quase certo.
  • Sua entrada e saída devem (obviamente) funcionar com o que a ideona fornecer para a E / S (stdin / stdout apenas parece). Tudo bem se você precisar de aspas na entrada N ou se a saída for algo como ans = ...etc.
  • Inclua um link para um trecho de código do seu código com o seu Nmax como entrada.
  • Se houver um empate (provavelmente apenas se vários envios atingirem o limite de caracteres de saída de 64 kB), a resposta mais alta vence.
Passatempos de Calvin
fonte
3
Mmm ... torta ambígua.
Dennis
Isso pode muito facilmente ser um código de golfe e, se for o caso, seria mais divertido.
Optimizer
2
@Optimizer Poderia ser código-golfe, mas seria muito parecido com qualquer outro código-geração de dígitos. Queria experimentar um concurso com base no tempo que pudesse ser verificado on-line. (Embora um problema computacionalmente mais complexo poderia ter sido melhor.)
de Calvino Hobbies
se este é APL golf código provavelmente iria ganhar (menos a parte arbitrária precisão)
TwiNight
1
Eu suspeito que esses programas serão totalmente vinculados à entrada / saída tentando gravar os dígitos no stdout. Cinco segundos é muito tempo para algo como triturador de y .
Will

Respostas:

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,

O Ideone não parece ter sido gmpy2instalado, o que é lamentável por pelo menos dois motivos. Um, porque tornaria a calcuação muito mais rápida, e dois, porque torna impraticável qualquer fórmula que exija uma raiz quadrada de precisão arbitrária.

A fórmula que eu uso para π foi listada por Ramanujan como Fórmula (39):

que converge à taxa de ~ 5,89 dígitos por termo. Que eu saiba, essa é a série convergente mais rápida do tipo que não requer a avaliação de uma raiz quadrada de precisão arbitrária. A fórmula (44) do mesmo artigo (taxa de convergência ~ 7,98 dígitos por termo) é mais frequentemente referida como a fórmula de Ramanujan.

A fórmula que eu uso para e é a soma dos fatoriais inversos. O número de termos requeridos é calculado como -1 ( 10 n ), usando uma aproximação que encontrei no excesso de matemática . O componente Lambert W 0 é encontrado usando o método de Newton.

O cálculo de cada um desses somatórios é feito através da Avaliação Rápida da Função E (mais conhecida como divisão binária), originalmente criada por Karatsuba. O método reduz um somatório para n termos a um único valor racional p / q . Esses dois valores são multiplicados para produzir o resultado final.

Atualização: a
criação de perfil revelou que mais da metade do tempo necessário para o cálculo foi gasto na divisão final. Somente os log ( 2 n ) bits mais altos de q são necessários para obter precisão total. O código agora preenche o buffer de saída Ideone em 3,33s .

Atualização 2:
Como esse é um desafio de , decidi escrever minha própria rotina de divisão para combater a lentidão do CPython. A implementação divnracima usa a Divisão Newton-Raphson . A idéia geral é calcular d = 1 / q · 2 n usando o Método de Newton, onde n é o número de bits que o resultado requer e calcula o resultado como p · d >> n . O tempo de execução agora é de 2,87s - e isso sem cortar os bits antes do cálculo; é desnecessário para este método.

primo
fonte
4

PARI / GP: 33000

Este é basicamente o programa fornecido no OEIS , modificado para receber e formatar a saída corretamente. Deve servir como uma linha de base a ser vencida, se nada mais.

I supor isso é preciso. Eu verifiquei em 100 e 20k contra o OEIS, e foi compatível com ambos. É muito difícil encontrar mais dígitos online para comparar.

Para 33.000, são necessários cerca de 4,5 segundos, portanto provavelmente poderia ser um pouco esbarrado. Eu apenas me cansei de mexer no input e no loop de submissão / compilação / execução do 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);
}

Link para Ideone.com

Geobits
fonte
Seus dígitos correspondem aos meus, então eu vou sair do ramo e dizer que provavelmente estão corretos.
primo
Este programa passa essencialmente todo o seu tempo no loop, gerando dígitos um a um. Se você pegar Str(floor(frac(x)*10^m), vai centenas / milhares de vezes mais rápido.
Charles
2

Python 3

Como os pi e e internos não têm dígitos suficientes, decidi calcular os meus.

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])

Link IDEOne

Saída para STDIN = 1000:

8.5397342226735669504281233688422467204743749305568824722710929852470173635361001388261308713809518841081669216573834376992066672804294594807026609638293539437286935503772101801433821053915371716284188665787967232464763808892618434263301810056154560438283877633957941572944822034479453916753507796910068912594560500836608215235605783723340714760960119319145912948480279651779184356994356172418603464628747082162475871780202868607325544781551065680583616058471475977367814338295574582450942453416002008665325253385672668994300796223139976640645190237481531851902147391807396201201799703915343423499008135819239684881566321559967077443367982975103648727755579256820566722752546407521965713336095320920822985129589997143740696972018563360331663471959214120971348584257396673542429063767170337770469161945592685537660073097456725716654388703941509676413429681372333615691533682226329180996924321063261666235129175134250645330301407536588271020457172050227357541822742441070313522061438812060477519238440079
Beta Decay
fonte
Nmax é o maior valor de entrada que você pode fornecer ao seu programa antes que o ideone não o execute mais.
Hobbies de Calvin
1
@ Calvin'sHobbies Acho que o nmax é arbitrariamente grande ...
Decay Beta 23/14
1
A ideona não lhe dá poder computacional infinito. Qual é o maior valor de entrada que seu programa pode executar em ideone? (Embora, de fato, seu programa não segue a should be able to work for any N, given unlimited resourcesregra maior parte da produção é zeros em torno N = 10000..)
de Calvino Hobbies
Isso não é python3: NameError: name 'xrange' not defined.
Bakuriu 23/09/14
2

Scala - 1790

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

Usamos a fórmula de Wetherfield para π (e o código da fórmula de Machin é transportado a partir daqui ). Calculamos e usando a série de potências comuns.

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
fonte