5 Sekunden, um Kuchen zu finden

11

Pi mal e (oder Pie, wenn Sie mehrdeutige Notation mögen) auf 100 Dezimalstellen ist:

8.5397342226735670654635508695465744950348885357651149618796011301792286111573308075725638697104739439...

( OIES A019609 ) ( Argument für mögliche Irrationalität )

Ihre Aufgabe ist es, ein Programm zu schreiben, das eine positive ganze Zahl N aufnimmt und Pi * e auf N Dezimalstellen abgeschnitten ausgibt. zB wenn N = 2, dann sollte der Ausgang sein 8.53.

Dies ist ein Optimierungsproblem, daher gewinnt die Übermittlung, die die richtige Ausgabe für den höchsten Wert von N liefern kann.

Um sicherzustellen, dass alle Einsendungen mit derselben Rechenleistung beurteilt werden, muss Ihr Code auf ideone ausgeführt werden und eine beliebige Sprache verwenden, die sie unterstützen. Gemäß der ideone-FAQ gibt es ein Laufzeitlimit von 5 Sekunden für nicht angemeldete Benutzer. Diese 5-Sekunden-Grenze müssen Sie verwenden, nicht die 15-Sekunden-Grenze für angemeldete Benutzer. ( Weitere Einschränkungen wie Speicher, Codegröße usw. finden Sie in der FAQ .)

Insbesondere sollte jeder, der nicht bei ideone angemeldet ist, in der Lage sein, Ihr Programm auf ideone für alle Werte von N von 1 bis zu einem maximalen Nmax auszuführen und fast immer die richtige Ausgabe zu sehen . ohne irgendwelche Time limit exceededoder Memory limit exceededusw. Fehler. Die Einreichung mit dem größten Nmax gewinnt.

(Ob die tatsächlich benötigte Zeit ein Smidge über 5 Sekunden ist, spielt keine Rolle, solange Ideone keine Fehler liefert. " Fast die ganze Zeit " wird als mehr als 95% der Zeit für ein bestimmtes N definiert.)

Einzelheiten

  • Sie können eine beliebige mathematische Methode verwenden, um Pi * e zu berechnen, aber Sie dürfen die Ausgabe nicht über die ersten Dutzend Stellen von Pi, e oder Pi * e hinaus fest codieren .
    • Ihr Programm sollte bei unbegrenzten Ressourcen für jedes N funktionieren können.
    • Sie können eingebaute Pi oder e Konstanten verwenden, wenn Ihre Sprache diese hat.
  • Sie dürfen nicht auf Websites oder Ressourcen außerhalb Ihres Codes zugreifen (sofern ideone dies überhaupt zulässt).
  • Abgesehen von der Hardcodierung und dem Zugriff auf externe Ressourcen ist alles, was ideone zulässt, mit ziemlicher Sicherheit in Ordnung.
  • Ihre Eingabe und Ausgabe muss (offensichtlich) mit dem Ideone funktionieren, das für I / O vorgesehen ist (nur stdin / stdout, wie es scheint). Es ist in Ordnung, wenn Sie Anführungszeichen um den Eingang N benötigen oder der Ausgang so etwas wie ans = ...usw. ist.
  • Bitte fügen Sie einen Link zu einem Ideone-Snippet Ihres Codes mit Ihrem Nmax als Eingabe hinzu.
  • Wenn es ein Unentschieden gibt (nur wahrscheinlich, wenn mehrere Einsendungen die 64-KB-Ausgabezeichengrenze erreichen), gewinnt die Antwort mit den höchsten Stimmen.
Calvins Hobbys
quelle
3
Mmm ... mehrdeutiger Kuchen.
Dennis
Dies kann sehr leicht ein Code-Golf sein und würde eher mehr Spaß machen, wenn es ist.
Optimierer
2
@Optimizer Es könnte Code-Golf sein, aber dann wäre es so ziemlich wie jedes andere Code-Golf der Zifferngenerierung. Ich wollte einen zeitbasierten Wettbewerb ausprobieren, der online überprüft werden kann. (Obwohl ein rechnerisch komplexeres Problem vielleicht besser gewesen wäre.)
Calvins Hobbys
Wenn dies Code Golf ist, würde APL wahrscheinlich gewinnen (abzüglich des willkürlichen Präzisionsteils)
TwiNight
1
Ich vermute, dass diese Programme vollständig an E / A gebunden sind und versuchen, die Ziffern in stdout zu schreiben. Fünf Sekunden sind eine lange Zeit für so etwas wie einen Y-Cruncher .
Wird

Antworten:

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 scheint nicht gmpy2installiert zu sein, was aus mindestens zwei Gründen bedauerlich ist. Erstens, weil es die Berechnung viel schneller machen würde, und zweitens, weil es jede Formel, die eine Quadratwurzel mit beliebiger Genauigkeit erfordert, unpraktisch macht.

Die Formel, die ich für π verwende, wurde von Ramanujan als Formel (39) aufgeführt:

die mit einer Rate von ~ 5,89 Stellen pro Term konvergiert . Meines Wissens ist dies die am schnellsten konvergierende Reihe ihrer Art, für die keine Quadratwurzel mit beliebiger Genauigkeit ausgewertet werden muss. Formel (44) in dem gleichen Papier (Konvergenzrate ~ 7,98 Digits pro Term) wird am häufigsten als bezeichnen die Ramanujans Formel.

Die Formel, die ich für e verwende, ist die Summe der inversen Fakultäten. Die Anzahl der erforderlichen Terme wird mit Γ -1 ( 10 n ) berechnet , wobei eine Näherung verwendet wird, die ich bei mathoverflow gefunden habe . Die Lambert W 0 -Komponente wird nach der Newtonschen Methode gefunden.

Die Berechnung jeder dieser Summierungen erfolgt über die schnelle E-Funktionsbewertung (allgemein bekannt als binäre Aufteilung), die ursprünglich von Karatsuba entwickelt wurde. Das Verfahren reduziert eine Summation auf n Terme auf einen einzelnen rationalen Wert p / q . Diese beiden Werte werden dann multipliziert, um das Endergebnis zu erhalten.

Update: Die
Profilerstellung ergab, dass mehr als die Hälfte der für die Berechnung benötigten Zeit in der Endabteilung verbracht wurde. Es werden nur die obersten log 2 (10 n ) Bits von q benötigt, um die volle Genauigkeit zu erhalten, daher schneide ich vorher einige ab. Der Code füllt jetzt den Ideone-Ausgabepuffer in 3,33 Sekunden .

Update 2:
Da dies eine , habe ich beschlossen, meine eigene Divisionsroutine zu schreiben, um die Langsamkeit von CPython zu bekämpfen. Die Implementierung von divnroben verwendet Newton-Raphson Division . Die allgemeine Idee besteht darin, d = 1 / q · 2 n unter Verwendung der Newtonschen Methode zu berechnen , wobei n die Anzahl der Bits ist, die das Ergebnis benötigt, und das Ergebnis als p · d >> n zu berechnen . Die Laufzeit beträgt jetzt 2,87 Sekunden - und dies ohne Abhacken von Bits vor der Berechnung; es ist für diese Methode nicht notwendig.

primo
quelle
4

PARI / GP: 33000

Dies ist im Grunde das bei OEIS bereitgestellte Programm , das so geändert wurde, dass die Eingabe korrekt übernommen und formatiert wird. Es sollte als Basis dienen, um zu schlagen, wenn nichts anderes.

Ich gehe davon aus, dass dies korrekt ist. Ich habe es bei 100 und 20.000 gegen OEIS überprüft und es passte für beide. Es ist ziemlich schwierig, online weitere Ziffern zu finden, gegen die man prüfen kann.

Für 33.000 dauert es ungefähr 4,5 Sekunden, so dass es wahrscheinlich ein wenig gestoßen werden könnte. Ich habe es einfach satt, mit der Eingabe und der langsamen Submit / Compile / Run-Schleife von ideone herumzuspielen.

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

Ideone.com Link

Geobits
quelle
Ihre Ziffern stimmen mit meinen überein, also werde ich mich auf die Beine stellen und sagen, dass sie wahrscheinlich richtig sind.
Primo
Dieses Programm verbringt im Wesentlichen die gesamte Zeit in der Schleife und generiert nacheinander Ziffern. Wenn Sie nur nehmen Str(floor(frac(x)*10^m), geht es hunderte / tausende Male schneller.
Charles
2

Python 3

Da die eingebauten pi und e nicht genügend Ziffern haben, habe ich beschlossen, meine eigenen zu berechnen.

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 Link

Ausgabe für STDIN = 1000:

8.5397342226735669504281233688422467204743749305568824722710929852470173635361001388261308713809518841081669216573834376992066672804294594807026609638293539437286935503772101801433821053915371716284188665787967232464763808892618434263301810056154560438283877633957941572944822034479453916753507796910068912594560500836608215235605783723340714760960119319145912948480279651779184356994356172418603464628747082162475871780202868607325544781551065680583616058471475977367814338295574582450942453416002008665325253385672668994300796223139976640645190237481531851902147391807396201201799703915343423499008135819239684881566321559967077443367982975103648727755579256820566722752546407521965713336095320920822985129589997143740696972018563360331663471959214120971348584257396673542429063767170337770469161945592685537660073097456725716654388703941509676413429681372333615691533682226329180996924321063261666235129175134250645330301407536588271020457172050227357541822742441070313522061438812060477519238440079
Beta-Zerfall
quelle
Nmax ist der größte Eingabewert, den Sie Ihrem Programm geben können, bevor ideone es nicht mehr ausführt.
Calvins Hobbys
1
@ Calvin'sHobbies Ich denke, dass nmax willkürlich groß ist ...
Beta Decay
1
ideone bietet Ihnen keine unendliche Rechenleistung. Was ist der größte Eingabewert, den Ihr Programm auf ideone ausführen kann? (Obwohl Ihr Programm tatsächlich nicht der should be able to work for any N, given unlimited resourcesRegel folgt . Der größte Teil der Ausgabe besteht aus Nullen bei etwa N = 10000.)
Calvins Hobbys
Das ist nicht python3 : NameError: name 'xrange' not defined.
Bakuriu
2

Scala - 1790

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

Wir verwenden die Wetherfield-Formel für π (und den von hier grob portierten Machin- Formelcode ). Wir berechnen e mit den gewöhnlichen Potenzreihen.

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
quelle