5 sekund na znalezienie ciasta

11

Pi razy e (lub Pie, jeśli lubisz niejednoznaczny zapis) do 100 miejsc po przecinku, wynosi:

8.5397342226735670654635508695465744950348885357651149618796011301792286111573308075725638697104739439...

( OIES A019609 ) ( argument za możliwą nieracjonalnością )

Twoim zadaniem jest napisanie programu, który przyjmuje dodatnią liczbę całkowitą N i wyprowadza Pi * e obcięte do N miejsc dziesiętnych. np. jeśli N = 2, to wynik powinien wynosić 8.53.

Jest to problem związany z optymalizacją, więc wygrana zostanie przesłana, która może dać poprawny wynik dla najwyższej wartości N.

Aby upewnić się, że wszystkie zgłoszenia są oceniane przy użyciu tej samej mocy obliczeniowej, kod musi być uruchamiany na ideone , w dowolnym obsługiwanym języku. Zgodnie z ideone faq , dla niezalogowanych użytkowników istnieje 5 sekundowy limit czasu działania. Tego 5-sekundowego limitu należy użyć, a nie 15-sekundowego limitu dla zalogowanych użytkowników. (Zobacz często zadawane pytania na temat innych ograniczeń, takich jak pamięć, rozmiar kodu itp.)

W szczególności każdy, kto nie jest zalogowany do ideone, powinien mieć możliwość uruchomienia programu na ideone dla wszystkich wartości N od 1 do pewnego maksymalnego momentu obrotowego Nmax i widzieć prawidłowe wyjście prawie przez cały czas . bez Time limit exceededlub Memory limit exceededitp błędy. Zgłoszenie z największym Nmax wygrywa.

(To, czy faktyczny czas zajmuje smidge w ciągu 5 sekund, nie ma znaczenia, dopóki ideone nie powoduje błędów. „ Prawie cały czas ” jest definiowany jako ponad 95% czasu dla dowolnego konkretnego N.)

Detale

  • Możesz użyć dowolnej metody matematycznej, którą chcesz obliczyć Pi * e, ale nie możesz zakodować wyjścia poza pierwszymi tuzinami cyfr Pi, e lub Pi * e .
    • Twój program powinien być w stanie pracować dla dowolnego N, biorąc pod uwagę nieograniczone zasoby.
    • Możesz użyć wbudowanych stałych Pi lub e, jeśli Twój język je posiada.
  • Nie możesz uzyskiwać dostępu do stron internetowych ani zasobów zewnętrznych w stosunku do twojego kodu (jeśli ideone na to zezwala).
  • Oprócz stałego kodowania i dostępu do zasobów zewnętrznych wszystko, co pozwala ideone, jest prawie na pewno w porządku.
  • Twoje dane wejściowe i wyjściowe muszą (oczywiście) współpracować z tym, co ideone zapewnia we / wy (wydaje się, że tylko stdin / stdout). W porządku, jeśli potrzebujesz cudzysłowów wokół wejścia N lub jeśli dane wyjściowe są podobne ans = ..., itp.
  • Podaj link do fragmentu ideonu kodu wraz z danymi wejściowymi Nmax.
  • Jeśli zdarzy się remis (możliwe tylko, jeśli wiele zgłoszeń osiągnie limit znaków wyjściowych 64kB), odpowiedź najwyższych głosów wygrywa.
Hobby Calvina
źródło
3
Mmm ... niejednoznaczne ciasto.
Dennis
To może być bardzo łatwy do gry w golfa i wolałby więcej zabawy.
Optymalizator
2
@Optimizer To może być golf golfowy, ale wtedy byłby prawie jak każda inna generacja cyfrowego golfa. Chciałem wypróbować konkurs oparty na czasie, który można zweryfikować online. (Chociaż bardziej skomplikowany obliczeniowo problem mógłby być lepszy.)
Calvin's Hobbies
jeśli jest to kod golfowy APL prawdopodobnie wygrałby (bez części o dowolnej precyzji)
TwiNight
1
Podejrzewam, że te programy będą całkowicie związane z IO, próbując zapisać cyfry na standardowe wyjście. Pięć sekund to długi czas na coś w rodzaju y-crunchera .
Czy

Odpowiedzi:

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,

Wygląda na to, że Ideone nie został gmpy2zainstalowany, co jest niefortunne z co najmniej dwóch powodów. Po pierwsze, ponieważ spowodowałoby to znacznie szybsze obliczanie, a po drugie, ponieważ sprawia, że ​​każda formuła wymagająca arbitralnej precyzji pierwiastka kwadratowego jest niepraktyczna.

Wzór, którego używam dla π, został wymieniony przez Ramanujana jako Wzór (39):

co zbiega się w tempie ~ 5,89 cyfr na termin. Według mojej wiedzy jest to najszybsza zbieżna seria tego rodzaju, która nie wymaga oceny pierwiastka kwadratowego o dowolnej precyzji. Wzór (44) w tym samym papierze (szybkość zbieżności ~ 7.98 cyfr na słowa) jest często nazywany w Ramanujana wzorze.

Wzór, którego używam dla e, jest sumą odwrotnych silni. Liczba wymaganych terminów jest obliczana jako Γ -1 ( 10 n ), przy użyciu przybliżenia, które znalazłem na matematycznym przepływie . Składnik Lambert W 0 można znaleźć przy użyciu metody Newtona.

Obliczenia każdego z tych podsumowań dokonuje się za pomocą Szybkiej oceny funkcji E (bardziej ogólnie nazywanej dzieleniem binarnym), pierwotnie opracowanej przez Karatsubę. Metoda redukuje sumowanie do n warunków do pojedynczej wartości wymiernej p / q . Te dwie wartości są następnie mnożone, aby uzyskać końcowy wynik.

Aktualizacja:
Profilowanie ujawniło, że ponad połowa czasu potrzebnego na obliczenia została poświęcona na ostateczny podział. Aby uzyskać pełną precyzję, potrzebne są tylko najwyższe z logów 2 (10 n ) bitów q , więc najpierw odcinam kilka. Kod wypełnia teraz bufor wyjściowy Ideone w 3,33 s .

Aktualizacja 2:
Ponieważ jest to wyzwanie związane z , postanowiłem napisać własną procedurę podziału, aby zwalczyć powolność CPython. Realizacja divnrpowyższa wykorzystuje Newton-Raphson Division . Ogólną ideą jest obliczenie d = 1 / q · 2 n przy użyciu metody Newtona, gdzie n jest liczbą bitów wymaganą przez wynik, i obliczenie wyniku jako p · d >> n . Czas działania wynosi teraz 2,87 s - i to nawet bez odcinania bitów przed obliczeniami; nie jest to konieczne w przypadku tej metody.

primo
źródło
4

PARI / GP: 33000

Jest to zasadniczo program podany w OEIS , zmodyfikowany w celu poprawnego pobierania i formatowania danych wyjściowych. Powinien służyć jako podstawa do pokonania, jeśli nic więcej.

Zakładam , że to jest poprawne. Sprawdziłem to przy 100 i 20k względem OEIS i pasowało do obu. Trudno jest znaleźć dalsze cyfry online, aby sprawdzić.

Dla 33 000 zajmuje to około 4,5 s, więc prawdopodobnie może być trochę zderzony. Właśnie zmęczyło mnie majstrowanie przy wejściu i powolnej pętli submission / compile / run 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 Ideone.com

Geobity
źródło
Twoje cyfry pasują do moich, więc wyjdę na kończynę i powiem, że prawdopodobnie są poprawne.
primo
Ten program spędza zasadniczo cały swój czas w pętli, generując cyfry jeden po drugim. Jeśli po prostu weźmiesz Str(floor(frac(x)*10^m), idzie setki / tysiące razy szybciej.
Charles
2

Python 3

Ponieważ wbudowane pi i e nie mają wystarczającej liczby cyfr, postanowiłem obliczyć własne.

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

Dane wyjściowe dla STDIN = 1000:

8.5397342226735669504281233688422467204743749305568824722710929852470173635361001388261308713809518841081669216573834376992066672804294594807026609638293539437286935503772101801433821053915371716284188665787967232464763808892618434263301810056154560438283877633957941572944822034479453916753507796910068912594560500836608215235605783723340714760960119319145912948480279651779184356994356172418603464628747082162475871780202868607325544781551065680583616058471475977367814338295574582450942453416002008665325253385672668994300796223139976640645190237481531851902147391807396201201799703915343423499008135819239684881566321559967077443367982975103648727755579256820566722752546407521965713336095320920822985129589997143740696972018563360331663471959214120971348584257396673542429063767170337770469161945592685537660073097456725716654388703941509676413429681372333615691533682226329180996924321063261666235129175134250645330301407536588271020457172050227357541822742441070313522061438812060477519238440079
Rozpad beta
źródło
Nmax to największa wartość wejściowa, jaką możesz podać swojemu programowi, zanim ideone go nie uruchomi.
Calvin's Hobbies
1
@ Calvin'sHobbies Myślę, że nmax jest jednak arbitralnie duży ...
Beta Decay
1
ideone nie daje nieskończonej mocy obliczeniowej. Jaka jest największa wartość wejściowa, którą Twój program może uruchomić na ideone? (Chociaż w rzeczywistości twój program nie przestrzega should be able to work for any N, given unlimited resourcesreguły. Większość danych wyjściowych to zera w okolicach N = 10000.)
Calvin's Hobbies
To nie python3: NameError: name 'xrange' not defined.
Bakuriu
2

Scala - 1790

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

Używamy formuły Wetherfielda dla π (i kod formuły Machina z grubsza przeniesiony stąd ). Obliczamy e za pomocą zwykłej serii mocy.

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
źródło