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 Pie, **but you may not hardcode the output beyond the first dozen digits of Pi, e or Pie**.
    • 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\$
11
  • 3
    \$\begingroup\$ Mmm... ambiguous pie. \$\endgroup\$
    – Dennis
    Sep 23, 2014 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, 2014 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\$ Sep 23, 2014 at 7:07
  • \$\begingroup\$ if this is code golf APL would probably win (minus the arbitrary precision part) \$\endgroup\$
    – TwiNight
    Sep 23, 2014 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, 2014 at 13:43

5 Answers 5

14
\$\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\$
2
  • \$\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, 2014 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, 2016 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\$
4
  • \$\begingroup\$ Nmax is the largest input value you can give your program before ideone will no longer run it. \$\endgroup\$ Sep 23, 2014 at 6:58
  • 1
    \$\begingroup\$ @Calvin'sHobbies I think that nmax is arbitrarily large though... \$\endgroup\$
    – Beta Decay
    Sep 23, 2014 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\$ Sep 23, 2014 at 7:47
  • \$\begingroup\$ That's not python3: NameError: name 'xrange' not defined. \$\endgroup\$
    – Bakuriu
    Sep 23, 2014 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\$
0
\$\begingroup\$

Java, 15,053

interface M{
  static void main(String[]a){
    int n=new Integer(new java.util.Scanner(System.in).nextInt()),i=1,x=n/100+9;
    var t=java.math.BigInteger.TEN.pow(n+=x).shiftLeft(1);
    var p=t;
    var s=t.TEN.pow(n-x);
    var e=s;
    for(;i<5e4;
      e=s.add(e.divide(s.valueOf(50001-i))))
      p=p.add(t=t.multiply(t.valueOf(i)).divide(t.valueOf(i+++i)));
    System.out.print((""+p.multiply(e)).substring(0,n-x));
  }
}

Outputs the digits without decimal dot after the leading 8.

Try it on ideone.

Increasing the 5e4/50001 will both increase the run-time and digit-accuracy. With 50,000 iterations, which runs in ~4-5 seconds usually, \$15\text{,}053\$ is the largest amount of digits we can reach with 100% accuracy. With an input of \$15\text{,}054\$, the last couple digits will be incorrect.

Uses a combination of my Java answer for the Calculate 500 digits of pi challenge and a port of this stackoverflow C# .NET answer with BigIntegers.

\$\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.