5 secondes pour trouver la tarte

11

Pi fois e (ou Pie si vous aimez la notation ambiguë) à 100 décimales est:

8.5397342226735670654635508695465744950348885357651149618796011301792286111573308075725638697104739439...

( OIES A019609 ) ( argument pour une possible irrationalité )

Votre tâche consiste à écrire un programme qui accepte un entier positif N et génère Pi * e tronqué à N décimales. par exemple, si N = 2, alors la sortie devrait être 8.53.

Ceci est un problème d'optimisation, donc la soumission qui peut donner la sortie correcte pour la valeur la plus élevée de N gagnera.

Pour garantir que toutes les soumissions sont jugées en utilisant la même puissance de calcul, votre code doit être exécuté sur ideone , dans n'importe quelle langue prise en charge. Selon la FAQ ideone , il y a une limite de 5 secondes pour les utilisateurs non connectés. Cette limite de 5 secondes est celle que vous devez utiliser, pas la limite de 15 secondes pour les utilisateurs connectés. (Voir la FAQ pour d'autres limites comme la mémoire, la taille du code, etc.)

Plus précisément, toute personne non connectée à ideone devrait pouvoir exécuter votre programme sur ideone pour toutes les valeurs de N de 1 à un certain Nmax maximum, et voir la sortie correcte presque tout le temps . sans aucune erreur Time limit exceededou Memory limit exceeded, etc. La soumission avec le plus grand Nmax gagne.

(Que le temps réel pris soit un peu plus de 5 secondes n'a pas d'importance tant que l'idéone ne donne pas d'erreurs. " Presque tout le temps " est défini comme plus de 95% du temps pour un N. particulier)

Détails

  • Vous pouvez utiliser n'importe quelle méthode mathématique que vous aimez pour calculer Pi * e, mais vous ne pouvez pas coder en dur la sortie au-delà des douze premiers chiffres de Pi, e ou Pi * e .
    • Votre programme devrait pouvoir fonctionner pour n'importe quel N, avec des ressources illimitées.
    • Vous pouvez utiliser des constantes Pi ou e intégrées si votre langue les possède.
  • Vous ne pouvez pas accéder à des sites Web ou à des ressources externes à votre code (si ideone le permet même).
  • Au-delà du codage en dur et de l'accès aux ressources externes, tout ce que l'idéone permet est presque certainement bien.
  • Vos entrées et sorties doivent (évidemment) fonctionner avec tout ce que l'ideone fournit pour les E / S (stdin / stdout seulement il semble). C'est bien si vous avez besoin de guillemets autour de l'entrée N ou la sortie est quelque chose comme ans = ..., etc.
  • Veuillez inclure un lien vers un extrait idéone de votre code avec votre Nmax en entrée.
  • S'il se trouve qu'il y a égalité (uniquement si plusieurs soumissions atteignent la limite de 64 Ko de caractères de sortie), la réponse la plus élevée gagne.
Loisirs de Calvin
la source
3
Mmm ... tarte ambiguë.
Dennis
Cela peut très facilement être un code-golf et serait plutôt plus amusant s'il l'est.
Optimizer
2
@Optimizer Il pourrait s'agir de golf de code, mais ce serait à peu près comme tous les autres golf de code de génération de chiffres. Je voulais essayer un concours basé sur le temps qui pourrait être vérifié en ligne. (Bien qu'un problème plus complexe en termes de calcul aurait pu être mieux.)
Hobbies de Calvin
si c'est le code, le golf APL gagnerait probablement (moins la partie de précision arbitraire)
TwiNight
1
Je soupçonne que ces programmes seront entièrement liés à IO en essayant d'écrire les chiffres sur stdout. Cinq secondes, c'est long pour quelque chose comme y-cruncher .
Will

Réponses:

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 ne semble pas avoir gmpy2installé, ce qui est regrettable pour au moins deux raisons. Un, car cela rendrait la calcification beaucoup plus rapide, et deux, car il rend impossible toute formule nécessitant une racine carrée de précision arbitraire.

La formule que j'utilise pour π a été répertoriée par Ramanujan comme formule (39):

qui converge au taux de ~ 5,89 chiffres par trimestre. À ma connaissance, il s'agit de la série convergente la plus rapide du genre qui ne nécessite pas l'évaluation d'une racine carrée de précision arbitraire. Formule (44) dans le même papier (taux de convergence ~ 7.98 chiffres par trimestre) est le plus souvent appelée la formule Ramanujan.

La formule que j'utilise pour e est la somme des factorielles inverses. Le nombre de termes requis est calculé comme Γ -1 ( 10 n ), en utilisant une approximation que j'ai trouvée sur mathoverflow . Le composant Lambert W 0 est trouvé en utilisant la méthode de Newton.

Le calcul de chacune de ces sommations se fait via l' évaluation rapide de la fonction E (plus généralement connue sous le nom de division binaire), conçue à l'origine par Karatsuba. La méthode réduit une sommation à n termes à une seule valeur rationnelle p / q . Ces deux valeurs sont ensuite multipliées pour produire le résultat final.

Mise à jour: Le
profilage a révélé que plus de la moitié du temps nécessaire au calcul a été passé dans la division finale. Seuls les bits logarithmiques supérieurs de 2 (10 n ) de q sont nécessaires pour obtenir une précision totale, donc j'en coupe quelques-uns au préalable. Le code remplit désormais le tampon de sortie Ideone en 3,33 s .

Mise à jour 2:
Comme il s'agit d'un défi d' , j'ai décidé d'écrire ma propre routine de division pour lutter contre la lenteur de CPython. La mise en œuvre de ce qui divnrprécède utilise la division Newton-Raphson . L'idée générale est de calculer d = 1 / q · 2 n en utilisant la méthode de Newton, où n est le nombre de bits requis par le résultat, et de calculer le résultat comme p · d >> n . Le temps d'exécution est maintenant de 2,87 s - et cela sans même couper les bits avant le calcul; c'est inutile pour cette méthode.

primo
la source
4

PARI / GP: 33000

Il s'agit essentiellement du programme donné à OEIS , modifié pour prendre correctement les entrées et formater les sorties. Il devrait servir de référence à battre, si rien d'autre.

Je suppose que c'est exact. Je l'ai vérifié à 100 et 20 km contre OEIS, et il correspondait aux deux. Il est assez difficile de trouver d'autres chiffres en ligne pour vérifier.

Pour 33 000, cela prend environ 4,5 secondes, donc il pourrait probablement être légèrement heurté. Je me suis juste lassé de tripoter l'entrée et la boucle de soumission / compilation / exécution lente d'idéone.

{ 
    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);
}

Lien Ideone.com

Géobits
la source
Vos chiffres correspondent aux miens, donc je vais sortir sur un membre et dire qu'ils sont probablement corrects.
primo
Ce programme passe essentiellement tout son temps dans la boucle, générant des chiffres un par un. Si vous le prenez, Str(floor(frac(x)*10^m)cela va des centaines / milliers de fois plus vite.
Charles
2

Python 3

Puisque les pi et e intégrés n'ont pas assez de chiffres, j'ai décidé de calculer les miens.

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

Lien IDEOne

Sortie pour STDIN = 1000:

8.5397342226735669504281233688422467204743749305568824722710929852470173635361001388261308713809518841081669216573834376992066672804294594807026609638293539437286935503772101801433821053915371716284188665787967232464763808892618434263301810056154560438283877633957941572944822034479453916753507796910068912594560500836608215235605783723340714760960119319145912948480279651779184356994356172418603464628747082162475871780202868607325544781551065680583616058471475977367814338295574582450942453416002008665325253385672668994300796223139976640645190237481531851902147391807396201201799703915343423499008135819239684881566321559967077443367982975103648727755579256820566722752546407521965713336095320920822985129589997143740696972018563360331663471959214120971348584257396673542429063767170337770469161945592685537660073097456725716654388703941509676413429681372333615691533682226329180996924321063261666235129175134250645330301407536588271020457172050227357541822742441070313522061438812060477519238440079
Beta Decay
la source
Nmax est la plus grande valeur d'entrée que vous puissiez donner à votre programme avant que ideone ne l'exécute plus.
Hobbies de Calvin
1
@ Calvin'sHobbies Je pense que nmax est arbitrairement grand ...
Beta Decay
1
ideone ne vous donne pas une puissance de calcul infinie. Quelle est la plus grande valeur d'entrée que votre programme puisse exécuter sur ideone? (Bien qu'en fait, votre programme ne respecte pas la should be able to work for any N, given unlimited resourcesrègle. La plupart des sorties sont des zéros à environ N = 10000.)
Hobbies Calvin
Ce n'est pas python3: NameError: name 'xrange' not defined.
Bakuriu
2

Scala - 1790

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

Nous utilisons la formule de Wetherfield pour π (et le code de la formule de Machin grossièrement porté à partir d' ici ). Nous calculons e en utilisant la série de puissance ordinaire.

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
la source