-2

I wrote scripts to calculate pi in python, perl and c. They all use the same algorithm (trapezoidal reimann sum of a circle with n subintervals) and the python and perl programs always get the same result when n is the same. However, the c program doesn't get the same answer like it should, it actually overestimates pi, which is impossible. What is wrong with the c program?

Python:

#!/usr/bin/python
n = 1000000
def f(x):
    return (1-(float(x)**2))**float(0.5)
val = 0
for i in range(n):
i = i+1
val = val+f(float(i)/float(n))
val = val*2
pi = (float(2)/n)*(float(1)+val)
print pi

Perl:

#!/usr/bin/perl
$n = 1000000;
$h = 0;
for($a = 1; $a < $n; ++$a){
    $t = $a/$n;
    $val = (1-(($t)**2))**0.5;
    $h+=$val;
}
$h = $h*2;
$pi = (2/$n)*(1+$h);
printf "%.11f", $pi;

C:

#include <stdio.h>
#include <math.h>
int main()
{
    int count, n = 1000000;
    double pi;
    double val = 0.0, p = 0.50, x = 1.0, four = 4.0;
    for(count = 1; count < n; ++count)
    {
      val += (pow(x-(((double)count/(long double)n)*((double)count/(long 
double)n)), p));
    }
    pi = (val + x) * (four/(long double)n);
    printf("%.11f", pi);
}

Results: C = 3.14159465241 Perl and Python = 3.14159265241

ikegami
  • 367,544
  • 15
  • 269
  • 518
Paul
  • 88
  • 9
  • 1
    Can you add the results you get as well? – cs95 Jun 29 '17 at 22:21
  • Perl and Python get 3.14159265241, C gets 3.14159465241 – Paul Jun 29 '17 at 22:24
  • 1
    I'm trying to understand why you have an extra `i = i+1` in your python loop and why it isn't there anywhere else. – cs95 Jun 29 '17 at 22:26
  • Don't worry about the python one, it works. The C one is what matters, the python one is there just to show the structure. I actually wrote the python script a while back so it might be kinda weird. – Paul Jun 29 '17 at 22:28
  • `$a = ++$a` wat. – melpomene Jun 29 '17 at 22:29
  • `$ ./a.out ; 3.14159465241` It works! – cs95 Jun 29 '17 at 22:31
  • 5
    Why are your programs all structured differently (use of variables vs no variables, choice of names, arithmetic)? It makes it harder to compare them. – melpomene Jun 29 '17 at 22:32
  • 3.14159265359 is pi. 3.14159465241 is an overestimate, which should be impossible. It should get the same result as the perl and python codes but it doesn't which is what I am confused about – Paul Jun 29 '17 at 22:33
  • Your C code isn't the same. It's barely similar. I took your Perl code and transliterated it to C (matching statement for statement but with C syntax) and got the same (correct) results as your Perl code. Your C code has an error in its implementation of the algorithm. – lurker Jun 29 '17 at 22:33
  • I translated your Perl code to C literally (use of global variables everywhere and nonsensical statements included): http://ideone.com/KqXIYG. It produces 3.14159265241, same as the Perl version. – melpomene Jun 29 '17 at 22:38
  • 1
    Why is there a `four = 4.0` in your C code when the number `4` doesn't appear in your Python or Perl code? – melpomene Jun 29 '17 at 22:41
  • Because that was the first C program I ever wrote – Paul Jun 29 '17 at 22:42
  • 1
    Writing C for the first time doesn't mean you have to introduce a `4` into the algorithm. Your C code clearly isn't using the same algorithm if it needs a `4`. – melpomene Jun 29 '17 at 22:43
  • 1
    Oh now I remember, that was me trying to fix it and trying to change everything to make it work and then forgetting to change it back. My bad – Paul Jun 29 '17 at 22:45
  • The Python code will certainly not work as is, sonce there is no identation after the `for` statement. Indentation matters, in Python. Notice also that `float(0.5)` is the same as `0.5`, and there is a function `math.sqrt`, that will likely be faster and more accurate. You could safely replace `float(i)/float(n)` with `i/n` in Python 3 (integer division yields a float), and since you pass a float to `f`, `float(x)` is useless if `f`. Anyway, the whole of `f` could and should be inlined. `float(1)` and `float(2)` are as well useless. There are a few other similar oddities. –  Jun 30 '17 at 04:38

1 Answers1

3

Well, comparing your methods, it became obvious your final operation for calculating pi was incorrect.

Replace pi = (val + x) * (four/(long double)n); with these two lines:

val = val * (long double)2.0;
pi = (val + x) * ((long double)2.0/(long double)n);

Compiling and running gives:

3.14159265241

Which I believe is the output you're looking for.

cs95
  • 379,657
  • 97
  • 704
  • 746