I'm trying to compare C and Fortran code for performance. For calculating pi using Leibniz's series, I've got the following Fortran code
program pi_leibniz
implicit none
integer, parameter :: dp=selected_real_kind(15,307)
integer :: k=0, precision=9
real(dp), parameter :: correct = 0.7853981633974483d0, eps = epsilon(real(1,dp))
real(dp) :: sum = 0.0, delta
character(8) :: fmt
logical, parameter :: explicit = .false.
real :: start, finish
delta = 10.**(-precision-1)*0.25
if (delta<eps) then
delta=eps
precision=14
print *, "Precision specified too high, reverting to double precision (14 digits)"
endif
write(fmt,'(A,I0,A,I0,A)') '(f',precision+2,'.',precision,')'
call cpu_time(start)
do
sum = sum + real((-1)**k,dp)/real(2*k+1,dp)
k = k+1
if (abs(sum-correct)<delta) exit
if (explicit) print fmt, 4.*sum
enddo
call cpu_time(finish)
print fmt, 4.*sum
print '(A,I0,A,I0,A)', "converged in ", k, " iterations with ", precision, " digits precision"
print '(g0,a)', finish-start," s"
end program pi_leibniz
and a nearly identical C code:
#include <stdio.h>
#include <time.h>
#include <float.h>
#include <math.h>
int main(void){
int precision=9;
size_t k=0;
const double correct=0.7853981633974483;
double sum=0.0, delta = 0.25*pow(10.0,-(precision+1));
clock_t start,finish;
double sgn = 1.0;
if (delta < DBL_EPSILON){
delta = DBL_EPSILON;
precision = 14;
printf("Precision specified too high, reverting to double precision (14 digits)\n");
}
start = clock();
for(k=0; fabs(sum-correct) >= delta; k++, sgn=-sgn)
sum += sgn/(2*k+1);
finish = clock();
printf("%.*f\n",precision,4*sum);
printf("converged in %zu iterations with %d digits precision\n",k,precision);
printf("%f s\n",(finish-start)/(double)CLOCKS_PER_SEC);
return 0;
}
I compile both with GNU compilers and just a -O2 option. Edit: 64-bit.
The Fortran code runs happily up to the full double precision, calculating the first 15 digits of pi in a couple of seconds on my machine. The C code performs even slightly faster than Fortran up to 8 decimal places, converging to the same digits in the same number of iterations; however, with precision=9
the Fortran code converges to 3.141592653 in 2.27s/1581043254 iterations, while the C code takes 12.9s/9858058108 iterations (~6x) and the last digit is off by 1. With higher precision, the time for Fortran is of the same order while C takes ~2 minutes to compute the first 11 digits of pi.
What could be the reason for the discrepancy and how do I avoid whatever is slowing down the C code?
Edit: I did as @pmg suggested and changed the loop in the C code, making the convergence monotonic:
for(k=0; fabs(sum-correct) > delta; k+=2)
sum += 1.0/(2*k+1) - 1.0/(2*k+3);
While this speeds up convergence somewhat at lower precision, it actually makes the C program essentially hang even at precision=8
now (takes more than 3 minutes to compute).
Edit 2: Since computing at precision>8
results in integer overflow, it seems that the correct way was to declare k
was as integer(8) :: k
in Fortran and unsigned long
in C. With this modification, the Fortran code now performs almost exactly as the C code for 10/11 digits of pi and seems to 'hang' at higher precision.
How come, then, that using an essentially incorrect method still produced the correct result before, and took the same amount of time to calculate whether it was 10 or 15 digits of pi? Just for fun, it took 1611454902 iterations to 'converge' to 3.14159265358979, which happens to be exactly pi to 14 decimal places.
Your Fortran code is incorrect.
You are likely using a default integer to be 32-bit and using HUGE(k)
you will see that the maximum integer value k
can take is 2147483647. In this case you will have
integer overflow happening with the iteration count and (before then) with its evaluation in real(2*k+1,dp)
.
Much as you use selected_real_kind
to find a real kind fitting your requirements, you use should selected_int_kind
to find a suitable integer kind. If we trust the C version, then the iteration count could reach such a large number that k
should have kind selected_int_kind(11)
.
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With