Code Golf: Leibniz formula for Pi

Question

I recently posted one of my favourite interview whiteboard coding questions in "What's your more controversial programming opinion", which is to write a function that computes Pi using the Leibniz formula.

It can be approached in a number of different ways, and the exit condition takes a bit of thought, so I thought it might make an interesting code golf question. Shortest code wins!

Given that Pi can be estimated using the function 4 * (1 - 1/3 + 1/5 - 1/7 + ...) with more terms giving greater accuracy, write a function that calculates Pi to within 0.00001.

Edit: 3 Jan 2008

As suggested in the comments I changed the exit condition to be within 0.00001 as that's what I really meant (an accuracy 5 decimal places is much harder due to rounding and so I wouldn't want to ask that in an interview, whereas within 0.00001 is an easier to understand and implement exit condition).

Also, to answer the comments, I guess my intention was that the solution should compute the number of iterations, or check when it had done enough, but there's nothing to prevent you from pre-computing the number of iterations and using that number. I really asked the question out of interest to see what people would come up with.

Answer 1

J, 14 chars

4*-/%>:+:i.1e6

---

Explanation

• 1e6 is number 1 followed by 6 zeroes (1000000).
• i.y generates the first y non negative numbers.
• +: is a function that doubles each element in the list argument.
• >: is a function that increments by one each element in the list argument.

So, the expression >:+:i.1e6 generates the first one million odd numbers:

1 3 5 7 ...

• % is the reciprocal operator (numerator "1" can be omitted).
• -/ does an alternate sum of each element in the list argument.

So, the expression -/%>:+:i.1e6 generates the alternate sum of the reciprocals of the first one million odd numbers:

1 - 1/3 + 1/5 - 1/7 + ...

• 4* is multiplication by four. If you multiply by four the previous sum, you have π.

That's it! J is a powerful language for mathematics.

---

Edit: since generating 9! (362880) terms for the alternate sum is sufficient to have 5 decimal digit accuracy, and since the Leibniz formula can be written also this way:

4 - 4/3 + 4/5 - 4/7 + ...

...you can write a shorter, 12 chars version of the program:

-/4%>:+:i.9!

Answer 2

Language: Brainfuck, Char count: 51/59

Does this count? =]

Because there are no floating-point numbers in Brainfuck, it was pretty difficult to get the divisions working properly. Grr.

Without newline (51):

+++++++[>+++++++<-]>++.-----.+++.+++.---.++++.++++.

With newline (59):

+++++++[>+++++++>+<<-]>++.-----.+++.+++.---.++++.++++.>+++.

Answer 3

Perl

26 chars

26 just the function, 27 to compute, 31 to print. From the comments to this answer.

sub _{$-++<1e6&&4/$-++-&_}       # just the sub
sub _{$-++<1e6&&4/$-++-&_}_      # compute
sub _{$-++<1e6&&4/$-++-&_}say _  # print

28 chars

28 just computing, 34 to print. From the comments. Note that this version cannot use 'say'.

$.=.5;$\=2/$.++-$\for 1..1e6        # no print
$.=.5;$\=2/$.++-$\for$...1e6;print  # do print, with bonus obfuscation

36 chars

36 just computing, 42 to print. Hudson's take at dreeves's rearrangement, from the comments.

$/++;$\+=8/$//($/+2),$/+=4for$/..1e6
$/++;$\+=8/$//($/+2),$/+=4for$/..1e6;print

About the iteration count: as far as my math memories go, 400000 is provably enough to be accurate to 0.00001. But a million (or as low as 8e5) actually makes the decimal expansion actually match 5 fractional places, and it's the same character count so I kept that.

Answer 4

Ruby, 33 characters

(0..1e6).inject{|a,b|2/(0.5-b)-a}

Answer 5

Another C# version:

(60 characters)

4*Enumerable.Range(0, 500000).Sum(x => Math.Pow(-1, x)/(2*x + 1));  // = 3,14159

Answer 6

52 chars in Python:

print 4*sum(((-1.)**i/(2*i+1)for i in xrange(5**8)))

(51 dropping the 'x' from xrange.)

36 chars in Octave (or Matlab):

l=0:5^8;disp((-1).^l*(4./(2.*l+1))')

(execute "format long;" to show all the significant digits.) Omitting 'disp' we reach 30 chars:

octave:5> l=0:5^8;(-1).^l*(4./(2.*l+1))'
ans = 3.14159009359631

Answer 7

Oracle SQL 73 chars

select -4*sum(power(-1,level)/(level*2-1)) from dual connect by level<1e6

Answer 8

Language: C, Char count: 71

float p;main(i){for(i=1;1E6/i>5;i+=2)p-=(i%4-2)*4./i;printf("%g\n",p);}

Language: C99, Char count: 97 (including required newline)

#include <stdio.h>
float p;int main(){for(int i=1;1E6/i>5;i+=2)p-=(i%4-2)*4./i;printf("%g\n",p);}

I should note that the above versions (which are the same) keep track of whether an extra iteration would affect the result at all. Thus, it performs a minimum number of operations. To add more digits, replace 1E6 with 1E(num_digits+1) or 4E5 with 4E(num_digits) (depending on the version). For the full programs, %g may need to be replaced. float may need to be changed to double as well.

Language: C, Char count: 67 (see notes)

double p,i=1;main(){for(;i<1E6;i+=4)p+=8/i/(i+2);printf("%g\n",p);}

This version uses a modified version of posted algorithm, as used by some other answers. Also, it is not as clean/efficient as the first two solutions, as it forces 100 000 iterations instead of detecting when iterations become meaningless.

Language: C, Char count: 24 (cheating)

main(){puts("3.14159");}

Doesn't work with digit counts > 6, though.

Answer 9

Haskell

I got it down to 34 characters:

foldl subtract 4$map(4/)[3,5..9^6]

This expression yields 3.141596416935556 when evaluated.

Edit: here's a somewhat shorter version (at 33 characters) that uses foldl1 instead of foldl:

foldl1 subtract$map(4/)[1,3..9^6]

Edit 2: 9^6 instead of 10^6. One has to be economical ;)

Edit 3: Replaced with foldl' and foldl1' with foldl and foldl1 respectively—as a result of Edit 2, it no longer overflows. Thanks to ShreevatsaR for noticing this.

Answer 10

23 chars in MATLAB:

a=1e6;sum(4./(1-a:4:a))

Answer 11

F#:

Attempt #1:

let pi = 3.14159

Cheating? No, its winning with style!

Attempt #2:

let pi =
    seq { 0 .. 100 }
    |> Seq.map (fun x -> float x)
    |> Seq.fold (fun x y -> x + (Math.Pow(-1.0, y)/(2.0 * y + 1.0))) 0.0
    |> (fun x -> x * 4.0)

Its not as compact as it could possibly get, but pretty idiomatic F#.

Answer 12

common lisp, 55 chars.

(loop for i from 1 upto 4e5 by 4 sum (/ 8d0 i (+ i 2)))

Answer 13

Mathematica, 27 chars (arguably as low as 26, or as high as 33)

NSum[8/i/(i+2),{i,1,9^9,4}]

If you remove the initial "N" then it returns the answer as a (huge) fraction.

If it's cheating that Mathematica doesn't need a print statement to output its result then prepend "Print@" for a total of 33 chars.

NB:

If it's cheating to hardcode the number of terms, then I don't think any answer has yet gotten this right. Checking when the current term is below some threshold is no better than hardcoding the number of terms. Just because the current term is only changing the 6th or 7th digit doesn't mean that the sum of enough subsequent terms won't change the 5th digit.

Answer 14

Using the formula for the error term in an alternating series (and thus the necessary number of iterations to achieve the desired accuracy is not hard coded into the program):

public static void Main(string[] args) {
    double tolerance = 0.000001;
    double piApproximation = LeibnizPi(tolerance);
    Console.WriteLine(piApproximation);
}

private static double LeibnizPi(double tolerance) {
    double quarterPiApproximation = 0;

    int index = 1;
    double term;
    int sign = 1;
    do {
        term = 1.0 / (2 * index - 1);
        quarterPiApproximation += ((double)sign) * term;
        index++;
        sign = -sign;
    } while (term > tolerance);

    return 4 * quarterPiApproximation;
}

Answer 15

C#:

public static double Pi()
{
    double pi = 0;
    double sign = 1;
    for (int i = 1; i < 500002; i += 2)
    {
        pi += sign / i;
        sign = -sign;
    }
    return 4 * pi;
}

Answer 16

Perl :

$i+=($_&1?4:-4)/($_*2-1)for 1..1e6;print$i

for a total of 42 chars.

Answer 17

Ruby, 41 chars (using irb):

s=0;(3..3e6).step(4){|i|s+=8.0/i/(i-2)};s

Or this slightly longer, non-irb version:

s=0;(3..3e6).step(4){|i|s+=8.0/i/(i-2)};p s

This is a modified Leibniz:

1. Combine pairs of terms. This gives you 2/3 + 2/35 + 2/99 + ...
2. Pi becomes 8 * (1/(1 * 3) + 1/(5 * 7) + 1/(9 * 11) + ...)

Answer 18

F# (Interactive Mode) (59 Chars)

{0.0..1E6}|>Seq.fold(fun a x->a+ -1.**x/(2.*x+1.))0.|>(*)4.

(Yields a warning but omits the casts)

Answer 19

Here's a solution in MUMPS.

pi(N)
 N X,I
 S X=1 F I=3:4:N-2 S X=X-(1/I)+(1/(I+2))
 Q 4*X

Parameter N indicates how many repeated fractions to use. That is, if you pass in 5 it will evaluate 4 * (1 - 1/3 + 1/5 - 1/7 + 1/9 - 1/11)

Some empirical testing showed that N=272241 is the lowest value that gives a correct value of 3.14159 when truncated to 5 decimal points. You have to go to N=852365 to get a value that rounds to 3.14159.

Answer 20

Javascript:

a=0,b=-1,d=-4,c=1e6;while(c--)a+=(d=-d)/(b+=2)

In javascript. 51 characters. Obviously not going to win but eh. :P

Edit -- updated to be 46 characters now, thanks to Strager. :)

---

UPDATE (March 30 2010)

A faster (precise only to 5 decimal places) 43 character version by David Murdoch

for(a=0,b=1,d=4,c=~4e5;++c;d=-d)a-=d/(b-=2)

Answer 21

C# using iterator block:

static IEnumerable<double> Pi()
{
    double i = 4, j = 1, k = 4;
    for (;;)
    {
        yield return k;
        k += (i *= -1) / (j += 2);
    }
}

Answer 22

For the record, this Scheme implementation has 95 characters ignoring unnecessary whitespace.

(define (f)
  (define (p a b)
    (if (> a b)
      0
      (+ (/ 1.0 (* a (+ a 2))) (p (+ a 4) b))))
  (* 8 (p 1 1e6)))

Answer 23

Here's a recursive answer using C#. It will only work using the x64 JIT in Release mode because that's the only JIT that applies tail-call optimisation, and as the series converges so slowly it will result in a StackOverflowException without it.

It would be nice to have the IteratePi function as an anonymous lambda, but as it's self-recursive we'd have to start doing all manner of horrible things with Y-combinators so I've left it as a separate function.

public static double CalculatePi()
{
    return IteratePi(0.0, 1.0, true);
}

private static double IteratePi(double result, double denom, bool add)
{
    var term = 4.0 / denom;
    if (term < 0.00001) return result;    
    var next = add ? result + term : result - term;
    return IteratePi(next, denom + 2.0, !add);
}

Answer 24

Language: dc, Char count: 35

dc -e '9k0 1[d4r/r2+sar-lad274899>b]dsbxrp'

Answer 25

Language: C99 (implicit return 0), Char count: 99 (95 + 4 required spaces)

exit condition depends on current value, not on a fixed count

#include <stdio.h>

float p, s=4, d=1;
int main(void) {
  for (; 4/d > 1E-5; d += 2)
        p -= (s = -s) / d;
  printf("%g\n", p);
}

compacted version

#include<stdio.h>
float
p,s=4,d=1;int
main(void){for(;4/d>1E-5;d+=2)p-=(s=-s)/d;printf("%g\n",p);}

Answer 26

Most of the current answers assume that they'll get 5 digits accuracy within some number of iterations and this number is hardcoded into the program. My understanding of the question was that the program itself is supposed to figure out when it's got an answer accurate to 5 digits and stop there. On that assumption here's my C# solution. I haven't bothered to minimise the number of characters since there's no way it can compete with some of the answers already out there, so I thought I'd make it readable instead. :)

    private static double GetPi()
    {
        double acc = 1, sign = -1, lastCheck = 0;

        for (double div = 3; ; div += 2, sign *= -1)
        {
            acc += sign / div;

            double currPi = acc * 4;
            double currCheck = Math.Round(currPi, 5);

            if (currCheck == lastCheck)
                return currPi;

            lastCheck = currCheck;
        }
    }

Answer 27

Java

    double pi=0,s=-1,i=1;
    for (;i<1E6;i+=2)pi+=((1d/i)*(s=-s)); 
    pi*=4;

Answer 28

Ruby:

irb(main):031:0> 4*(1..10000).inject {|s,x| s+(-1)**(x+1)*1.0/(2*x-1)}
=> 3.14149265359003

Answer 29

C# cheating - 50 chars:

static single Pi(){
  return Math.Round(Math.PI, 5));
}

It only says "taking into account the formula write a function..." it doesn't say reproduce the formula programmatically :) Think outside the box...

C# LINQ - 78 chars:

static double pi = 4 * Enumerable.Range(0, 1000000)
               .Sum(n => Math.Pow(-1, n) / (2 * n + 1));

C# Alternate LINQ - 94 chars:

static double pi = return 4 * (from n in Enumerable.Range(0, 1000000)
                               select Math.Pow(-1, n) / (2 * n + 1)).Sum();

And finally - this takes the previously mentioned algorithm and condenses it mathematically so you don't have to worry about keep changing signs.

C# longhand - 89 chars (not counting unrequired spaces):

static double pi()
{
  var t = 0D;
  for (int n = 0; n < 1e6; t += Math.Pow(-1, n) / (2 * n + 1), n++) ;
  return 4 * t;
}

Answer 30

64 chars in AWK:

~# awk 'BEGIN {p=1;for(i=3;i<10^6;i+=4){p=p-1/i+1/(i+2)}print p*4}'
3.14159

Answer 31

#!/usr/bin/env python
from math import *
denom = 1.0
imm = 0.0
sgn = 1
it = 0
for i in xrange(0, int(1e6)):
    imm += (sgn*1/denom)
    denom += 2
    sgn *= -1    
print str(4*imm)

Answer 32

C++

double LeibnizPi( double tolerance )
{
    double sum = 1.0;
    for( int plus_div = 5, minus_div = -3, limit = 10 / tolerance; plus_div < limit ; plus_div += 4, minus_div -= 4 )
        sum += 1./plus_div + 1./minus_div;
    return 4 * sum;
}

Answer 33

After noting that

(= (- (/ 4 n)
      (/ 4 (+ n 2)))
   (/ 8 n (+ n 2)))

or, in a more familiar notation:

4    4      8
- - --- = ------
n   n+2   n(n+2)

Common Lisp, with a do* loop (62 essential characters):

(do* ((n 1 (+ n 4))
      (p 8/3 (+ p (/ 8 n (+ n 2)))))
     ((< (- pi p) 1e-6)
      p)

with a tail recursive function (70 essential characters):

(defun l (n p)
  (if (< (- pi p) 1e-6)
      p
      (l (+ n 4)
          (+ p (/ 8 n (+ n 2))))))
(l 1 0)

and with the extended loop (86 essential characters):

(loop for n from 1 by 4
      sum (/ 8 n (+ n 2)) into p
      until (< (- pi p) 1e-6)
      finally (return p))

all under the presumption that preliminary checks how far we have to go to get the desired accuracy are cheating.

Answer 34

double d = 1;
double s = 1;
double pi = 0;

while(4.0 / d > 0.000001){
    pi += s*4.0/d;
    d+=2;
    s = -s;        
}
printf("%f\n", pi);

Answer 35

Uh....as a general rule in numeric processing one should sum series from the smallest term toward the largest to avoid trouble with loss of precision. So in

fortran77

stripped down (248 characters)

      function pi(n)
      pi=0.
      t=10**(-n-0.5)
      i=int(4/t)
      i=i/2
      s=-1.                     
      do 1 j=i*2+1,1,-2
         pi = pi + s/j
         s=-s
 1    continue
      pi=abs(pi)*4              
      return
      end

With a scaffold and comments (600 characters)

      program leibnitz

      n=5
      p=int(pi(n)*10.**n)/10.**n
      write(6,*)p 

      stop
      end

c     Returns pi computed to <n> digits by the leibniz formula
      function pi(n)
      pi=0.
c     find the smallest term we need by insuring that it is too small to
c     effect the interesting digits.
      t=10**(-n-0.5)
      i=int(4/t)
      i=i/2
      s=-1.                     ! sign of term, might be off by one, but
      do 1 j=i*2+1,1,-2
         pi = pi + s/j
         s=-s
 1    continue
      pi=abs(pi)*4              ! we fix the sign problem here
      return
      end

output:

   3.1415901

It seems to work for arbitrary number of digits up to 6ish where the precision of real runs out. It is not optimized for speed or for minimum number of operations.

Answer 36

Java

void pi(){
    double x=1,y=1,d=1;
    for(;x<1E6;) { y=-y;d+=y/((2*x++)+1); }
    System.out.println(d*4);
}

Answer 37

Python 3 (40 bytes)

sum(8/(n*(n+2))for n in range(1,5**8,4))

This version uses optimization from @Svante's answer.

print +7 bytes

print(sum(8/(n*(n+2))for n in range(1,5**8,4)))

Python 2.x +1 byte

sum(8./(n*(n+2))for n in range(1,5**8,4))

print +6 bytes

print sum(8./(n*(n+2))for n in range(1,5**8,4))

http://codepad.org/amtxUxKp

Answer 38

Lua, 46 characters

p=4 for i=3,9^6,4 do p=p-8/i/(i+2)end print(p)

Answer 39

VB 117 chars:

Function Pi()
  Dim t = 0D
  For n = 0 To 1000000
    t += Math.Pow(-1, n) / (2 * n + 1)
  Next
  Return 4 * t
End Function

VB LINQ 115 chars (omitting the unnecessary line continuation):

Function Pi()
  Return 4 * Enumerable.Range(0, 1000000) _
             .Sum(Function(n) Math.Pow(-1, n) / (2 * n + 1))
End Function

And then call:

Sub Main()
  Console.WriteLine("{0:n5}", Pi)
End Sub

Answer 40

Another VB solution, using the rather cool aggregation syntax:

Public ReadOnly Pi As Double = 4 * Aggregate i In Enumerable.Range(0, 100000) _
                                   Select (-1) ^ i / (i * 2 + 1) Into Sum()

Expression only: 74 characters without unnecessary whitespaces.

Answer 41

Erlang, ~126 chars:

-module (pi).
-export ([pi/0]).

pi() -> 4 * pi(0,1,1).
pi(T,M,D) ->
    A = 1 / D,
    if A > 0.00001 
              -> pi(T+(M*A), M*-1, D+2);
        true  -> T
    end.

Answer 42

Here's mine in C++, probably the longest way of doing it :P

double pi(){
   bool add = true;
   double rPi = 0;
   for(long i = 1; i < 99999999; i=i+2)
   {
            double y = (double) i;
            double x = (double) 1;
            if(add)
            {
                   rPi = rPi + (x/y);
                   add = false;
            }
            else
            {
                    rPi = rPi - (x/y);
                    add = true;
            }
   }
            return (rPi * (double) 4);
   }

Answer 43

PHP, 99 chars

$output =

${!${''}=function(){$pi=1;for($i=0;$i<pow(10,6);$i++){$pi+=($i%2?1:-1)/(3+$i*2);}return $pi*4;}}();

var_dump($output);

I guess there is some pretty tricks i don't know to reduce this answer ! Took the one-line output trick from this post.

Answer 44

I just sort of wrote this right after reading interview question in the topic on controversial opinion. It ain't pretty but it took me about 3-4 minutes and I am checking for accuracy in each loop. C++. I'll wake up tomorrow and post a solution that doesn't suck :)

double get_pi(int acc)
{

  double pi;
  double dynamicpart;
  int operationCoeff = 1;
  int denom = 3;
  while(1)
  { 
      dynamicpart =
         1/denom + operationCoeff*(denom+2);
      pi = 4*(1-dynamicpart);
      if(!(pi*acc*10-(int)pi*acc*10)) break;
)
      denom+=2;
      operationCoeff = -operationCoeff;
  }



}

Answer 45

R - 27 chars

sum(4/seq(1,1e6,2)*c(1,-1))

Answer 46

1 Character: . Written in "MySuperDuperDomainSpecificLanguageThatOnlyReturnsThisOneAnswerAndNothingElse".

Yes this is meant as a joke, but seriously unless you disallow DSLs then EVERY Code Golf contest could be won by some goober who writes his own language that uses one character to return just that one result...