5 Seconds to Find Pie

Question

Pi times e (or Pie if you like ambiguous notation) to 100 decimal places is:

8.5397342226735670654635508695465744950348885357651149618796011301792286111573308075725638697104739439...

(OIES A019609) (argument for possible irrationality)

Your task is to write a program that takes in a positive integer N, and outputs Pi*e truncated to N decimal places. e.g. if N = 2, then the output should be 8.53.

This is an optimization problem, so the submission that can give the correct output for the highest value of N will win.

To ensure all submissions are judged using the same computing power, your code must be run on ideone, using any language they support. According to the ideone faq, there is a 5 second run time limit for not logged in users. This 5 second limit is the one you must use, not the 15 second limit for logged in users. (See faq for other limits like memory, code size, etc.)

Specifically, anyone not logged into ideone should be able to run your program on ideone for all values of N from 1 to some maximum Nmax, and see correct output nearly all of the time. without any Time limit exceeded or Memory limit exceeded, etc. errors. The submission with the largest Nmax wins.

(Whether the actual time taken is a smidge over 5 seconds doesn't matter as long as ideone doesn't give errors. "Nearly all of the time" is defined as more than 95% of the time for any particular N.)

Details

• You may use any mathematical method you like to calculate Pi*e, but you may not hardcode the output beyond the first dozen digits of Pi, e or Pi*e.
  • Your program should be able to work for any N, given unlimited resources.
  • You may use built in Pi or e constants if your language happens to have them.
• You may not access websites or resources external to your code (if ideone even allows this).
• Beyond hardcoding and accessing external resources, anything that ideone allows is almost certainly fine.
• Your input and output must (obviously) work with whatever ideone provides for i/o (stdin/stdout only it seems). It's fine if you need quotes around the input N or the output is something like ans = ..., etc.
• Please include a link to an ideone snippet of your code with your Nmax as input.
• If there happens to be a tie (only likely if multiple submissions hit the 64kB output character limit) the highest votes answer wins.

Answer 1

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 doesn't seem to have gmpy2 installed, which is unfortunate for at least two reasons. One, because it would make the calcuation a lot faster, and two, because it makes any formula requiring an arbitrary precision square root impractical.

The formula I use for π was listed by Ramanujan as Formula (39):

which converges at the rate of ~5.89 digits per term. To my knowledge, this is the fastest converging series of its kind that doesn't require the evaluation of an arbitrary precision square root. Formula (44) in the same paper (convergence rate ~7.98 digits per term) is most often referred to as the Ramanujan Formula.

The formula I use for e is the sum of inverse factorials. The number of terms required is calculated as Γ^-1(10ⁿ), using an approximation I found on mathoverflow. The Lambert W₀ component is found using Newton's Method.

The calculation of each of these summations is done via the Fast E-function Evalution (more generally known as binary-splitting), originally devised by Karatsuba. The method reduces a summation to n terms to a single rational value p/q. These two values are then multiplied to produce the final result.

Update:
Profiling revealed that over half of the time needed for the calculation was spent in the final division. Only the upper most log₂(10ⁿ) bits of q are needed to obtain full precision, so I trim a few off beforehand. The code now fills the Ideone output buffer in 3.33s.

Update 2:
Since this is an optimization challenge, I decided to write my own division routine to combat CPython's slowness. The implementation of divnr above uses Newton-Raphson Division. The general idea is to calculate d = 1/q · 2ⁿ using Newton's Method, where n is the number of bits the result requires, and compute the result as p · d >> n. Runtime is now 2.87s - and this is without even chopping off bits before the calculation; it's unnecessary for this method.

Answer 2

PARI/GP : 33000

This is basically the program given at OEIS, modified to take input and format output correctly. It should serve as a baseline to beat, if nothing else.

I assume this is accurate. I checked it at 100 and 20k against OEIS, and it matched for both. It's quite hard to find further digits online to check against.

For 33,000 it takes about 4.5s, so it could probably be bumped a small bit. I just got tired of fiddling with the input and ideone's slow-ass submit/compile/run loop.

{ 
    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

Answer 3

Python 3

Since the built in pi and e don't have enough digits, I decided to calculate my own.

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

Output for STDIN = 1000:

8.5397342226735669504281233688422467204743749305568824722710929852470173635361001388261308713809518841081669216573834376992066672804294594807026609638293539437286935503772101801433821053915371716284188665787967232464763808892618434263301810056154560438283877633957941572944822034479453916753507796910068912594560500836608215235605783723340714760960119319145912948480279651779184356994356172418603464628747082162475871780202868607325544781551065680583616058471475977367814338295574582450942453416002008665325253385672668994300796223139976640645190237481531851902147391807396201201799703915343423499008135819239684881566321559967077443367982975103648727755579256820566722752546407521965713336095320920822985129589997143740696972018563360331663471959214120971348584257396673542429063767170337770469161945592685537660073097456725716654388703941509676413429681372333615691533682226329180996924321063261666235129175134250645330301407536588271020457172050227357541822742441070313522061438812060477519238440079

Answer 4

Scala - 1790

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

We use Wetherfield's formula for π (and Machin-formula code crudely ported from here). We calculate e using the ordinary power series.

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