Leibniz formula for π/4

Question

I am asked to print the summation of the Leibniz formula up to the nth term of the series correct to 15 decimal places.In Calculus, the Leibniz formula for π is given by: 1 - 1/3 + 1/5 -1/7 + ... = π/4

This is my code

#include<stdio.h>
#include<math.h>
int main()
{
    int n,i;
    long double s=0;
    scanf("%d",&n);
    for(i=0;i<n;i++){
        s+=(long double)pow(-1,i)/(2*i+1);
    }
    printf("%Lf\n",s);
    return 0;
}

Can anyone tell me why i can't reach precision upto 15th decimal place? My aim is not to print the value of pi/4, i just have to print the summation for a given n

What `n` are you entering? Did you read http://en.wikipedia.org/wiki/Leibniz_formula_for_%CF%80#Inefficiency from your link? — Jim Balter, Oct 17 '14 at 20:00
Did you read the `Inefficiency` section of the Wikipedia article you linked to? — Mark Dickinson, Oct 17 '14 at 20:00
The formula is for π/4. And that particular series converges *very* slowly. — Keith Thompson, Oct 17 '14 at 20:01
Also, `long double` mostly the alias for `double`, so it won't give more than 17 significant digits. But that's an other story in this case. — m0nhawk, Oct 17 '14 at 20:03
@m0nhawk: That depends on the implementation. On some systems, `long double` is substantially more precise than `double`. — Keith Thompson, Oct 17 '14 at 20:04
You're missing `#include ` and `#include `; your compiler should have warned you about undeclared functions. `pow` works on `double` values, use `powl` for `long double`. (Not that the extra precision of `long double` is doing you any good in this case anyway.) — Keith Thompson, Oct 17 '14 at 20:07
@KeithThompson that's why I put *mostly*, thought, I don't have a statistics what's more popular. — m0nhawk, Oct 17 '14 at 20:10
Check the answer from the holder of [π record](http://stackoverflow.com/a/14283481/1030110). There are pretty fast-convergent equalities for π. — m0nhawk, Oct 17 '14 at 20:11

chux - Reinstate Monica · Accepted Answer · 2014-10-18T21:32:11.410

Q: why ... reach precision upto 15th decimal place?
A: To display 15 decimal places after the decimal point use format "%0.15f". To calculate 15 decimal places of convergence, at a minimum, n will need to be very large.

As mentioned by @user3386109, "error in the result is bounded by 1 / (2n+1)", so about 5e14 iterations will be needed. (Rough estimate: 10 days on my PC.) As typical double has a precision of about 1 part in pow(2,53) or 1 in 9e15, the limits of double calculation are being reached. The below code compares the order of calculation to reduce error, but at best, the error will still be at least 0.5 parts in 9e15.

As the terms of the series oscillate about the limit, when stopping after n iterations, a final iteration of 1/2 the next iteration could be added. this would gain about 1 bit of accuracy.

As others have mentioned, other methods exist to calculate π that converge faster.

Updated per the good observation by @user3386109.

In summing the terms, code can sum them in various orders. The 2 below methods illustrate that a modestly more stable result is achieved earlier when summing the small terms together first. I`d expect at best only a 1 or 2 bit better answer.

This was re-done using float as the small improvement does not show up until the last few bits are stable. With double and this slowly convergent series, that would take too long.

//Leibniz formula for pi/4
typedef float fp;

fp LeibnizForward(unsigned n) {
  volatile fp sum = 0.0;
  fp sign = 1.0;
  unsigned i = 1;
  while (n-- > 0) {
    sum += sign / i;
    sign = -sign;
    i = (i + 2);
  }
  return sum;
}

fp LeibnizReverse(unsigned n) {
  volatile fp sum = 0.0;
  fp sign = 1.0;
  unsigned i = 2 * n - 1;
  if (n % 2 == 0)
    sign = -sign;
  while (n-- > 0) {
    sum += sign / i;
    sign = -sign;
    i = (i - 2);
  }
  return sum;
}

void PiTest(unsigned n) {
  printf("%u\n", n);
  static const fp pic = 3.1415926535897932384626433832795;
  const char *format = "%s %0.9f\n";
  printf(format, "pi-", nextafterf(pic,0));
  printf(format, "pi ", pic);
  printf(format, "pi+", nextafterf(pic,4));
  fp pif = LeibnizForward(n) * 4;
  printf(format, "pif", pif);
  fflush(stdout);
  fp pir = LeibnizReverse(n) * 4;
  printf(format, "pir", pir);
  fflush(stdout);
}

int main(void) {
  PiTest(0);
  PiTest(1);
  PiTest(10);
  PiTest(100);
  PiTest(1000);
  PiTest(10000);
  PiTest(100000);
  PiTest(1000000);
  PiTest(10000000);
  PiTest(100000000);

  return 0;
}

0
pi- 3.141592503
pi  3.141592741
pi+ 3.141592979
pif 0.000000000
pir 0.000000000
1
pi- 3.141592503
pi  3.141592741
pi+ 3.141592979
pif 4.000000000
pir 4.000000000
10
pi- 3.141592503
pi  3.141592741
pi+ 3.141592979
pif 3.041839600
pir 3.041839600
25
pi- 3.141592503
pi  3.141592741
pi+ 3.141592979
pif 3.181576490
pir 3.181576729
100
pi- 3.141592503
pi  3.141592741
pi+ 3.141592979
pif 3.131592512
pir 3.131592751
1000
pi- 3.14 1592503
pi  3.14 1592741
pi+ 3.14 1592979
pif 3.14 0592575
pir 3.14 0592575
10000
pi- 3.141 592503
pi  3.141 592741
pi+ 3.141 592979
pif 3.141 498566
pir 3.141 492605
100000
pi- 3.1415 92503
pi  3.1415 92741
pi+ 3.1415 92979
pif 3.1415 85827
pir 3.1415 82489
1000000
pi- 3.14159 2503
pi  3.14159 2741
pi+ 3.14159 2979
pif 3.14159 5364
pir 3.14159 1549
10000000
pi- 3.14159 2503 previous float
pi  3.14159 2741 machine float pi
pi+ 3.14159 2979 next float
pif 3.14159 6794
pir 3.14159 2503
100000000 Additional iterations do not improve the result.
pi- 3.14159 2503
pi  3.14159 2741 
pi+ 3.14159 2979
pif 3.14159 6794
pir 3.14159 2503
1000000000
pi- 3.14159 2503
pi  3.14159 2741
pi+ 3.14159 2979
pif 3.14159 6794
pir 3.14159 2503

The arrows seem to indicate that you think the answer is correct. However, note for example with n=100000 you have pi=3.14159 but with reverse Leibniz, the answer is 3.14158 — user3386109, Oct 17 '14 at 21:19
Could you explain in two words the reason of such difference? I suppose it's all because of specific representation of `float-pointing` variables in memory, but don't sure in it.. — yulian, Oct 17 '14 at 21:20
@Yulian Khlevnoy No - not in 2 words. The effect is small yet measurable. This slowly converging series has other issues though that swamp the effect I was attempting to show. — chux - Reinstate Monica, Oct 17 '14 at 21:35
@chux, it seems to be very interesting case. Can you give several links or something else that will help to understand that effect? — yulian, Oct 17 '14 at 21:44
@Yulian Khlevnoy Here is 1: http://stackoverflow.com/questions/6699066/in-which-order-should-floats-be-added-to-get-the-most-precise-result and http://en.wikipedia.org/wiki/Kahan_summation_algorithm — chux - Reinstate Monica, Oct 17 '14 at 22:06

score 0 · Answer 2 · answered Oct 17 '14 at 20:04

0

Change your printf statement to print more of the values after the decimal point.

As an example, using

printf("%30.28Lf\n",s);

will print the number to 28 places after the decimal point. With n=25, I got the following output.

0.7953941713587578038252220991

answered Oct 17 '14 at 20:04

R Sahu

204,454
14
159
270

2

Which is still a very poor approximation to π/4. – Keith Thompson Oct 17 '14 at 20:05
2

The error in the result is bounded by `1 / (2n+1)`. So to get a result accurate to 15 decimal places, you would need `n=10^15`. But if you did that, the accumulated rounding errors would invalidate the answer. The Leibniz formula just flat out cannot be used to calculate a result accurate to 15 decimal places. – user3386109 Oct 17 '14 at 20:34

Leibniz formula for π/4

2 Answers2