0

I'm looking for a Leibniz approximation of PI in C.

I've tested this function but it does not seem to work:

float leibnizPi(float n){
    double pi=1.0;
    int i;
    int N;
    for (i=3, N=2*n+1; i<=N; i+=2)
        pi += ((i&2) ? -1.0 : 1.0) / i;
    return 4*pi;
}
pointerless
  • 753
  • 9
  • 20
Unicone
  • 228
  • 2
  • 8
  • 1
    I think the `&` should be `%` for modulo operation – bhow Oct 17 '17 at 19:33
  • 1
    What is your output? Show a [mcve]. – Jabberwocky Oct 17 '17 at 19:33
  • 1
    @bhow `i` is always odd, so your suggestion would produce a constant one. – Tom Karzes Oct 17 '17 at 19:39
  • @That's right, my mistake. It seems I did not read through it carefully enough. – bhow Oct 17 '17 at 19:43
  • 1
    The `i&2` part is correct. That will be true for i=3, 7, 11,... but not 5, 9, 13.... Why is `n` a `float` when it's clearly just a count? Otherwise it looks fine. Why do think there's an error? In other words, what inputs are you giving it, what outputs do you get, and what do you think they should be? – Lee Daniel Crocker Oct 17 '17 at 19:48
  • I actually think this code works as shown. It just takes many iterations for this formula to converge. Try setting `n` to 1000000 and you will see it's pretty close to pi. – Tom Karzes Oct 17 '17 at 19:48
  • @TomKarzes.: Sorry for the hasty answer which ofcourse lead to errors. I didn't correlate with the leibnit'z expansion. Sorry for the inconvenience. Thanks for the criticism. I at first didnt get my fault. – user2736738 Oct 17 '17 at 19:52

2 Answers2

2

I've tested you function and it seems to work fine. For n=1000, I get a result of 3.142592.

How I've tested the function:

#include <stdio.h>

float leibnizPi(float n) {
    double pi=1.0;
    int i;
    int N;
    for (i=3, N=2*n+1; i<=N; i+=2)
        pi += ((i&2) ? -1.0 : 1.0) / i;
    return 4*pi;
}

int main(void)
{
    printf("Pi: %f\n", leibnizPi(1000));
    return 0;
}
or523
  • 141
  • 6
  • 1
    I think this is right. I tried it with `n=1000000` and it seems to converge, it just takes a lot of terms. – Tom Karzes Oct 17 '17 at 19:50
0

OP's code interesting does its computation using double math and then after a potentially long loop returns a float. Rather than truncate to a float. Recommend to return a double. With float, the result is not expected to be more accurate than 1 part in maybe 223.

Belows shows that unnecessary impact which becomes evident after n > about 100,000 and final about 600,000.

float leibnizPif(float n){
    double pi=1.0;
    int i;
    int N;
    for (i=3, N=2*n+1; i<=N; i+=2)
        pi += ((i&2) ? -1.0 : 1.0) / i;
    return 4*pi;
}

double leibnizPid(float n){
    double pi=1.0;
    int i;
    int N;
    for (i=3, N=2*n+1; i<=N; i+=2)
        pi += ((i&2) ? -1.0 : 1.0) / i;
    return 4*pi;
}

int main() {

  double pi = acos(-1);
  float f1 = pi;
  float f2 = nextafterf(pi, 0);
  printf("%e\n", f1-f2);
  printf("                     pi  %.12e\n", pi);
  for (unsigned i=1; i<29; i++) {
    unsigned n = 1u << i;
    float f = leibnizPif(n);
    double d = leibnizPid(n);
    printf("%9u f % .5e %.12e   ", n, (f-pi)/pi, f);
    printf("d % .5e %.12e\n",  (d-pi)/pi, d);

  }
  printf("                     pi  %.12e\n", acos(-1));
  return 0;
}

Output

            error                pi
                     pi  3.141592653590e+00
        2 f  1.03474e-01 3.466666698456e+00   d  1.03474e-01 3.466666666667e+00
        4 f  6.30540e-02 3.339682579041e+00   d  6.30540e-02 3.339682539683e+00
        8 f  3.52602e-02 3.252365827560e+00   d  3.52602e-02 3.252365934719e+00
       16 f  1.87080e-02 3.200365543365e+00   d  1.87080e-02 3.200365515410e+00
       32 f  9.64357e-03 3.171888828278e+00   d  9.64354e-03 3.171888735237e+00
       64 f  4.89682e-03 3.156976461411e+00   d  4.89679e-03 3.156976358911e+00
      128 f  2.46747e-03 3.149344444275e+00   d  2.46748e-03 3.149344475125e+00
      256 f  1.23857e-03 3.145483732224e+00   d  1.23856e-03 3.145483689446e+00
      512 f  6.20513e-04 3.143542051315e+00   d  6.20487e-04 3.143541969477e+00
     1024 f  3.10574e-04 3.142568349838e+00   d  3.10546e-04 3.142568263114e+00
     2048 f  1.55377e-04 3.142080783844e+00   d  1.55349e-04 3.142080696509e+00
     4096 f  7.76643e-05 3.141836643219e+00   d  7.76934e-05 3.141836734621e+00
     8192 f  3.88840e-05 3.141714811325e+00   d  3.88514e-05 3.141714709002e+00
    16384 f  1.94559e-05 3.141653776169e+00   d  1.94269e-05 3.141653685021e+00
    32768 f  9.74187e-06 3.141623258591e+00   d  9.71375e-06 3.141623170237e+00
    65536 f  4.88485e-06 3.141607999802e+00   d  4.85695e-06 3.141607912146e+00
   131072 f  2.45634e-06 3.141600370407e+00   d  2.42849e-06 3.141600282926e+00
   262144 f  1.24208e-06 3.141596555710e+00   d  1.21425e-06 3.141596468272e+00
   524288 f  6.34955e-07 3.141594648361e+00   d  6.07127e-07 3.141594560935e+00
  1048576 f  3.31391e-07 3.141593694687e+00   d  3.03564e-07 3.141593607263e+00
  2097152 f  1.79610e-07 3.141593217850e+00   d  1.51782e-07 3.141593130427e+00
  4194304 f  1.03719e-07 3.141592979431e+00   d  7.58910e-08 3.141592892008e+00
  // No further improvement in float as precision is reached.
  8388608 f  2.78275e-08 3.141592741013e+00   d  3.79455e-08 3.141592772799e+00
 16777216 f  2.78275e-08 3.141592741013e+00   d  1.89728e-08 3.141592713194e+00
 33554432 f  2.78275e-08 3.141592741013e+00   d  9.48662e-09 3.141592683393e+00
 67108864 f  2.78275e-08 3.141592741013e+00   d  4.74328e-09 3.141592668491e+00
134217728 f  2.78275e-08 3.141592741013e+00   d  2.37136e-09 3.141592661040e+00
268435456 f  2.78275e-08 3.141592741013e+00   d  1.18567e-09 3.141592657315e+00
536870912 f  2.78275e-08 3.141592741013e+00   d  5.92798e-10 3.141592655452e+00
                     pi  3.141592653590e+00
                                              double keeps getting better

Modified code output

// with changes to use 64-bit integer math for the loops, 
// the result continue to get better returning double.
// It just takes longer and longer. 
1073741824 f  2.78275e-08 3.141592741013e+00   d  2.95864e-10 3.141592654519e+00
2147483648 f  2.78275e-08 3.141592741013e+00   d  1.47599e-10 3.141592654053e+00
4294967296 f  2.78275e-08 3.141592741013e+00   d  7.34476e-11 3.141592653821e+00
8589934592 f  2.78275e-08 3.141592741013e+00   d  3.63406e-11 3.141592653704e+00
chux - Reinstate Monica
  • 143,097
  • 13
  • 135
  • 256