11
\$\begingroup\$

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.
\$\endgroup\$
  • 3
    \$\begingroup\$ Mmm... ambiguous pie. \$\endgroup\$ – Dennis Sep 23 '14 at 4:16
  • \$\begingroup\$ This can very easily be a code-golf and would rather be more fun if it is. \$\endgroup\$ – Optimizer Sep 23 '14 at 6:53
  • 2
    \$\begingroup\$ @Optimizer It could be code-golf but then it would be pretty much like every other digit generation code-golf. I wanted to try a time based contest that could be verified online. (Though a more computationally complex problem might have been better.) \$\endgroup\$ – Calvin's Hobbies Sep 23 '14 at 7:07
  • \$\begingroup\$ if this is code golf APL would probably win (minus the arbitrary precision part) \$\endgroup\$ – TwiNight Sep 23 '14 at 7:24
  • 1
    \$\begingroup\$ I suspect that these programs will be entirely IO bound trying to write out the digits to stdout. Five seconds is a long time to something like y-cruncher. \$\endgroup\$ – Will Sep 23 '14 at 13:43
12
\$\begingroup\$

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(10n), using an approximation I found on mathoverflow. The Lambert W0 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 log2(10n) 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 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 · 2n 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.

\$\endgroup\$
4
\$\begingroup\$

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

\$\endgroup\$
  • \$\begingroup\$ Your digits match mine, so I'm going to go out on a limb and say they're probably correct. \$\endgroup\$ – primo Sep 25 '14 at 15:03
  • \$\begingroup\$ This program spends essentially all of its time in the loop, generating digits one by one. If you just take Str(floor(frac(x)*10^m) it goes hundreds/thousands of times faster. \$\endgroup\$ – Charles Apr 22 '16 at 13:46
2
\$\begingroup\$

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
\$\endgroup\$
  • \$\begingroup\$ Nmax is the largest input value you can give your program before ideone will no longer run it. \$\endgroup\$ – Calvin's Hobbies Sep 23 '14 at 6:58
  • 1
    \$\begingroup\$ @Calvin'sHobbies I think that nmax is arbitrarily large though... \$\endgroup\$ – Beta Decay Sep 23 '14 at 7:35
  • 1
    \$\begingroup\$ ideone doesn't give you infinite computing power. What is the largest input value your program can run on ideone? (Though in fact your program does not follow the should be able to work for any N, given unlimited resources rule. Most of the output is zeros at around N = 10000.) \$\endgroup\$ – Calvin's Hobbies Sep 23 '14 at 7:47
  • \$\begingroup\$ That's not python3: NameError: name 'xrange' not defined. \$\endgroup\$ – Bakuriu Sep 23 '14 at 16:57
2
\$\begingroup\$

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) 
}
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.