Efficient Approximation of Pi in C: A Comparison with Fortran | Help Needed

[fluidistic]

So this isn't for homework. Basically I have a program in fortran 90 that computes an approximation of pi using a very non efficient method. So that it can do about 10^8 operations and still doesn't converge that much which is what I'm looking for.
I want to write this program in C too, but I've never written any program in C yet, not even the hello world. I want to compare the speed of both programs for the task of approximating pi via a product.
In fortran the code is

program pi
implicit none

real (8) ::p
integer ::n,m

write(*,*)'Enter n please'
read(*,*)n

p=1d0
do m=1,n
p=p*((2d0*m)**2d0)/((2d0*m)**2d0-1d0)
end do
write(*,*)"The approximation of pi is",2*p
end program

Executing it gives:

 Enter n please
1000000
 The approximation of pi is   3.1415918681920378

I'd like to have the same characteristics for the program in C. Namely, I want to be able to enter n whenever I execute the program.
I do not know how to assignate a value to a real number (float in C?) that can change. So #assignate wouldn't work I think. I'd need this in order to start the product, p must equal 1.
I've been through some pages of a tutorial on the internet, so far my program is:

#include <stdio.h>
using namespace std;
int n,m
float p
m=1
int main()
{
 prinf ("Enter n.\n");
 scanf( "%d", n );
  getchar();
  return 0;
  p=1;
  do{
    
  }   while(m<n);
  
}

I know the lines "m=1" and "p=1" are probably wrong. Any help is appreciated.

 

[camel-man]

First off, take out using namespace std, that is C++, also return 0; goes at the end of your main function returning success and control back to OS. Then your scanf needs an ampersand before your variable ex: scanf("%d",&n);

 

[Coin]

Also most of your lines lack semicolons at the end, also "m=1" is outside of any function, also your do {} while; loop, there is nothing in between { and }, this means the loop will never end because the result of the conditional will never change.

The C equivalent of

do m=1,n
<stuff>
end

is

for(int m = 1; m <= n; m++) {
<stuff>
}

And the C equivalent of a**b is pow(a,b); (you will need to include math.h)

I bet you can find a computer program that just plain translates fortran into C.

 

[fluidistic]

Ok guys thanks a lot!
My code is now

#include <stdio.h>
#include <math.h>

int main()
{
  int n,m;
float p;
 prinf ("Enter n.\n");
 scanf( "%d", &n );
  getchar();

  p=1;
  for( m = 1; m <= n; m++) {
p=p*(pow((2*m),2)/(pow((2*m),2)-1));
}
prinf("The approximation of pi is\n",p);

getchar();
}

To Coin: yes I know I had put nothing between the do while because I was unsure on how to assignate values to p, n, m, etc. I think it's good to learn the basics of C, I wouldn't like to use a program to do this. :)
When I try to compile, I get

/tmp/ccppllAb.o: In function `main':
pi.c:(.text+0x11): undefined reference to `prinf'
pi.c:(.text+0x91): undefined reference to `prinf'
collect2: ld returned 1 exit status

So it seems something's wrong with my prinf statement.

 

[TylerH]

It's "printf." :P

 

[camel-man]

Also on this line "prinf("The approximation of pi is\n",p);"
it is printf like Tyler said, also you want to display your variable so you have to use %f inside your quotations to display a float, you can truncate the decimal places if you want ex. printf("The approximation of pi is %.4f\n",p); and you don't need the getchar(); there you can put it at the end of your code to keep the console window open if you are using dev-c++ or some other IDE that doesn't keep it open, btw what compiler are you using? EDIT: nvm I see why you have it there to take in the '\n' character after first scanf

 

[fluidistic]

Ok guys thanks, my program now works.
It's currently:

#include <stdio.h>
#include <math.h>

int main()
{
  int n,m;
float p;
 printf ("Enter n.\n");
 scanf( "%d", &n );


  p=1;
  for( m = 1; m <= n; m++) {
p=p*(pow((2*m),2)/(pow((2*m),2)-1));
}
printf("The approximation of pi is %.8f\n",2*p);

getchar();
}

However when I execute it, I get:

Enter n.
1000000
The approximation of pi is 3.14132977

But with the same n under fortran, I reached 3.1415918681920378 which is much closer to the "real" value of pi.
Also in C, when I increase n even more than 10^6, I get the same result, 3.14132977. Hmm... what's going on?

 

[camel-man]

try making p a double instead of a float ex double p; then in your printf statements use %lf to display a double

 

[fluidistic]

Ok thanks.
Now I get

pi.c: In function ‘main’:
pi.c:7: error: two or more data types in declaration specifiers

Line 7 is p=1; line.

 

[camel-man]

post all of your code

 

[fluidistic]

#include <stdio.h>
#include <math.h>

int main()
{
  int n,m;
float double p;
 printf ("Enter n.\n");
 scanf( "%d", &n );  p=1;
  for( m = 1; m <= n; m++) {
p=p*(pow((2*m),2)/(pow((2*m),2)-1));
}
printf("The approximation of pi is %lf\n",2*p);

getchar();
}

 

[camel-man]

take off float on line 7 it should just be double

int main()
{
  int n,m;
double p;
 printf ("Enter n.\n");
 scanf( "%d", &n );


  p=1;
  for( m = 1; m <= n; m++) {
p=p*(pow((2*m),2)/(pow((2*m),2)-1));
}
printf("The approximation of pi is %.10lf\n",2*p);

getchar();
	return 0;
}

that gave me the same answer as your fortran program

 

[fluidistic]

Ok thank you very much for all.
For C, n=10^9 around 19 s. (3.1415926423)
For fortran, n=10^9 around 26 s.(3.1415926422888196)
I don't know if it's because fortran kept more digits, in all cases I'm very impressed.

 

[camel-man]

you can always move out to 16 decimal places also just change .10lf to .16lf

 

[jhae2.718]

You could try a long double for more precision.

 

[TylerH]

If you're thinking about doing some numerical stuff in C, you might want to check out GMP. It's a library that allows you to use arbitrary precision integers and reals.

 

[Mark44]

Ok guys thanks a lot!
My code is now

#include <stdio.h>
#include <math.h>

int main()
{
  int n,m;
float p;
 prinf ("Enter n.\n");
 scanf( "%d", &n );
  getchar();

  p=1;
  for( m = 1; m <= n; m++) {
p=p*(pow((2*m),2)/(pow((2*m),2)-1));
}
prinf("The approximation of pi is\n",p);

getchar();
}

To Coin: yes I know I had put nothing between the do while because I was unsure on how to assignate values to p, n, m, etc. I think it's good to learn the basics of C, I wouldn't like to use a program to do this. :)
When I try to compile, I get

/tmp/ccppllAb.o: In function `main':
pi.c:(.text+0x11): undefined reference to `prinf'
pi.c:(.text+0x91): undefined reference to `prinf'
collect2: ld returned 1 exit status

So it seems something's wrong with my prinf statement.

I would do the loop without using pow(), since all you're doing is squaring a number (which is equivalent to multiplying it by itself).

I have introduced another variable, but haven't declared it. Your code would need to declare this variable.

for( m = 1; m <= n; m++) 
{
    value = 2.0 * m;
    p = p * value * value / (value * value - 1.0) ;
}

This code has at least two advantages:
1) It's simpler, by computing an intermediate value first, which makes the 2nd line easier to read, and easier to spot syntax errors caused by unmatched parentheses.
2) It uses multiplication and division, which can be performed more quickly than a relatively more expensive call to pow(). pow() should never be called if all you're doing is squaring a number.

I would also change the type of p from float to double. A float only has 6 or 7 decimal places of precision, so that will affect your answer.

 

[fluidistic]

I now get (with 16 decimals) 3.1415926422818243, still 19 s.
I'll retry tomorrow and I'll check out your posts guys more in details, it's 2:26 am currently... must sleep.

 

[Mark44]

You could try a long double for more precision.

Some compilers support long double, and some don't, or at least didn't in the past. I haven't looked into this recently, so it could be that modern compilers support long doubles.

A compiler that supports long double stores the number in 10 bytes, which is conveniently the same size as the registers in the floating point unit.

 

[Mark44]

I now get (with 16 decimals) 3.1415926422818243[/color], still 19 s.

Is that 19 seconds? Still, if you're using float, you don't have more than 6 to 7 digits of precision after the decimal point. From the calculator, \pi is 3.1415926535897932384626433832795

I'll retry tomorrow and I'll check out your posts guys more in details, it's 2:26 am currently... must sleep.

 

[jhae2.718]

Some compilers support long double, and some don't, or at least didn't in the past. I haven't looked into this recently, so it could be that modern compilers support long doubles.

A compiler that supports long double stores the number in 10 bytes, which is conveniently the same size as the registers in the floating point unit.

I know gcc does*, and it looks like Visual Studio does on 64-bit systems, and maps long double to double on 32-bit systems.

*Another option if using gcc is the nonstandard __float128 type, which is an IEEE 754 binary128 implementation.

 

[camel-man]

Is that 19 seconds? Still, if you're using float, you don't have more than 6 to 7 digits of precision after the decimal point. From the calculator, \pi is 3.1415926535897932384626433832795

He is using double

 

[rcgldr]

For C, n=10^9 around 19 s. (3.1415926423)
For fortran, n=10^9 around 26 s.(3.1415926422888196)
I don't know if it's because fortran kept more digits.

The difference in speed could be due to the C compiler optimizing your code by looking for duplicated calculations (common expressions) so it only calculates the value (2 m)² one time using a temp variable, while the fortran compiler calculates it twice. To confirm this you could change the fortran program to:

program pi
implicit none

real (8) ::p,z
integer ::n,m

write(*,*)'Enter n please'
read(*,*)n

p=1d0
do m=1,n
z=(2d0*m)**2d0
p=p*(z)/(z-1d0)
end do
write(*,*)"The approximation of pi is",2*p
end program

Additionally, C compiler might use add instead of multiply time 2 and multiply instead of pow(...,2). The fortan program for this would be:

program pi
implicit none

real (8) ::p,y,z
integer ::n,m

write(*,*)'Enter n please'
read(*,*)n

p=1d0
do m=1,n
y=m+m
z=y*y
p=p*(z)/(z-1d0)
end do
write(*,*)"The approximation of pi is",2*p
end program

This still may be slower than the C version because the fortran code may store results in y and z, while the C compiler is using temporary variables that don't have to be stored separately. Even if you converted the second example to C, the C compliler could realize that y and z are only used to update p and ingore them, only using temporary values instead.

It's also possible that your fortran compiler has an optimization enabling switch that you're not using.

If you're thinking about doing some numerical stuff in C, you might want to check out GMP. It's a library that allows you to use arbitrary precision integers and reals.

apfloat might be faster than GMP. Apfloat includes a test program to do a very optimized calculation of pi; on my system, it can calculate 1,000,000 digits of pi in about 7 seconds.

Wiki link to GMP description and links to GMP libraries:

http://en.wikipedia.org/wiki/GNU_Multi-Precision_Library

Link to apfloat

http://www.apfloat.org

 

[fluidistic]

Ok guys thank you for all.
If I have an opportunity to program more in my degree, I might use C over fortran so that at least I learn the basics of C which seems important.
I'm using gcc as C compiler and gfortran for fortran. I'm guessing they are both the latest versions.
I'm using double precision for both fortran and c.
My code for c is

#include <stdio.h>
#include <math.h>

int main()
{
  int n,m;
  double p;
 printf ("Enter n.\n");
 scanf( "%d", &n );


  p=1;
  for( m = 1; m <= n; m++) {
p=p*(pow((2*m),2)/(pow((2*m),2)-1));
}
printf("The approximation of pi is %.16lf\n",2*p);

getchar();
}

I've redone the test, gave me around 19/20 s for c and 25 s for fortran.
I tried your suggestion for C Mark, with the code:

#include <stdio.h>
#include <math.h>

int main()
{
  int n,m,value;
  double p;
 printf ("Enter n.\n");
 scanf( "%d", &n );


  p=1;
for( m = 1; m <= n; m++) 
{
    value = 2.0 * m;
    p = p * value * value / (value * value - 1.0) ;
}
printf("The approximation of pi is %.16lf\n",2*p);

getchar();
}

For n=10^9, I waited around 1 minute and 20 s, it had still not returned the approximation of pi so I aborted it.
And for any code, the value returned for the approximation isn't very good not because I'm using float/double, but because the method does not converge fastly at all. (Wallis product, I think it's called).
rcgldr, I tried both of your suggestions for fortran code and they gave me around 25/26 s.
My cpu is a core 2 duo E6300 (so 1.86 Ghz), although the programs use only 1 core.
I'd be interested to know if anyone else get similar results on different computers.

 

[rcgldr]

#include <stdio.h>
  int n,m,value;
  double p;

Value should be a double:

#include <stdio.h>
  int n,m;
  double p, value;

I modified the code a bit, using a fixed value for n to make sure what it is using. I used visual studio to compile and ran it from a dos console window to compare results. My result is different than yours, both of my programs calculate pi ~= 3.1415926430662622. The first program took 40.44 seconds, the second version took 6.12 seconds.

This one took 40.44 seconds:

#include <stdio.h>
#include <math.h>

int main()
{
int n,m;
double p;

    n = 1000000000;
    p = 1.0;
    for(m = 1; m <= n; m++) {
        p=p*(pow((2*m),2)/(pow((2*m),2)-1));
    }
    printf("The approximation of pi is %.16lf\n",2*p);
}

This one took 6.12 seconds:

#include <stdio.h>
#include <math.h>

int main()
{
int n,m;
double p, y, z;

    n = 1000000000;
    p = 1.0;
    for(m = 1; m <= n; m++) {
        y = (double)(m+m);
        z = y*y;
        p=p*(z /(z-1));
    }
    printf("The approximation of pi is %.16lf\n",2*p);
}

 

[fluidistic]

Value should be a double:

#include <stdio.h>
  int n,m;
  double p, value;

Ok, right.

I modified the code a bit, using a fixed value for n to make sure what it is using. I used visual studio to compile and ran it from a dos console window to compare results. My result is different than yours, both of my programs calculate pi ~= 3.1415926430662622. The first program took 40.44 seconds, the second version took 6.12 seconds.

I copied and pasted your programs. Both programs give exactly the same to me:

The approximation of pi is 3.1415926422818243

and were both completed in about 20 s.

 

[rcgldr]

vs2005 and vc/c++ 4.0 (or perhaps Windows XP) are defaulting the FPU (floating point unit) precision control to double precision. I modded the code to force FPU to extended precision. vc 4.0 got the same answer as gcc. vs2005 partially unfolds the loop and maybe other optimizations to get the faster speed.

The (z - 1) term is where most of the truncation error is occurring.

Windows caculator, 4 tan^-1(1):
The approximation of pi is 3.1415926535897932384626433832795

vs2005, double precision: 6.1 seconds
The approximation of pi is 3.1415926430662622

vs2005, extended precision: 6.6 seconds (closest result)
The approximation of pi is 3.1415926528075948

vc/c++ 4.0, double precision 15.1 seconds
The approximation of pi is 3.1415926430662622

vc/c++ 4.0, extended precision 16.2 seconds
The approximation of pi is 3.1415926422818243

vc/c++ 2.2 extended precision 21.2 seconds (16 bit real mode code)
The approximation of pi is 3.1415926451151450

code to set precision control to extended precision. I changed "int" to "long" so I could use vc 2.2 which is 16 bit compiler.

#include <stdio.h>

static short fcw; /* floating point control word */

int main()
{
long n,m;
double p, y, z;

/* set precision control to extended precision */

__asm{
    fnstcw fcw
    or fcw,0300h
    fldcw fcw
}

    n = 1000000000;
    p = 1.0;
    for(m = 1; m <= n; m++) {
        y = (double)(m+m);
        z = y*y;
        p *= (z / (z-1.));
    }
    printf("The approximation of pi is %.16lf\n",2.*p);
}

For visual studio. this can also be done using _controlfp() function:

#include <float.h>
/* ... */
int main()
{
/* ... */
    _controlfp(_PC_64, _MCW_PC);
/* ... */
}