1

Note: LaTeX isn't supported on this site. I'm not sure if there is a better way to write math equations other than to write them in code.

I'm writing a Fortran program to estimate pi through the summation of series:

A = Sum of a_i from i=1 to N

where 

pi/4 = 1 - 1/3 + 1/5 - 1/7 + 1/9 - 1/11 + 1/13 ...

To compute pi through the series summation, the suggested approach is to set

a_i = (-1)^(i+1)/(2i-1)

To do this, I wrote the following Fortran program -

program cpi
double precision val, pi
integer i
num = 1000
val = 0
do i = 1, num
  val = val + ((-1)**(i+1))/(2*i-1)
end do
pi = val
print *, 'Estimated Value of PI:', pi
end program cpi

When I run this program, the output is

   Estimated Value of PI:   1.0000000000000000

I must have made a mistake (likely in the /(2*i-1)). I new to Fortran and don't know what I did wrong.

Axion004
  • 943
  • 1
  • 11
  • 37

1 Answers1

1

I see my mistake! I need to write out 1.d0 and 2.d0 instead of 1 and 2 so that the calculations are evaluated in double format. I also forgot to multiply pi = val*4.d0. Changing the cpi program to

program cpi
double precision val, pi, MFLOPS
integer i, T1, T2
num = 1000000
val = 0.d0
call system_clock(T1)       ! get time stamp
do i = 1, num
  val = val + ((-1.d0)**(i+1.d0))/(2.d0*i-1.d0)
end do
call system_clock(T2)       ! get time stamp
MFLOPS = num*2.d0/((T2-T1)*1.d8)     ! compute MFlop/sec rate
pi = val*4.d0
print *, 'Estimated Value of PI:', pi
print *, 'The calculated number of MFLOPS is:', MFLOPS
end program cpi

returns

Estimated Value of PI:   3.1415916535897743     
The calculated number of MFLOPS is:   3.0303030303030304E-002

I also added a MFLOPS calculation to see the computational speed.

Axion004
  • 943
  • 1
  • 11
  • 37
  • You're not evaluating in double format. 1.0, 2.0, and 4.0 are single precision literal constants. The right hand side of the expression is computed in single precision. You need to have 1.d0, 2.d0, and 4.d0. You also do not want to use a real-valued exponent in (-1.0)**(i+1.0) as this may not very efficient. – evets Jan 19 '19 at 19:45
  • Is there an alternative way to write (-1)^(i+1)? I updated my program to evaluate this as (-1.d0)**(i+1.d0). – Axion004 Jan 19 '19 at 19:53
  • I think this part (-1)**(i+1) is okay; just changing (2*i-1) to (2.0d0 * i - 1.0d0) makes the ratio to be evaluated in double precision. (And, comparing the computational speed may be interesting...) – roygvib Jan 19 '19 at 19:59
  • Yes, don't compute (-1)**(i+1) at all. Add `sgn` to the double precision declaration. Set it to 1 outside the loop. You now sum `val = val + sgn / (2*i-1)`, and then do `sgn = -sgn` for the next iteration. Note the denominator uses integers, because `sgn` is double precision you get automatic conversion of the denominator to the correct precision. – evets Jan 19 '19 at 20:02
  • Changing `val = val + sgn / (2*i-1)` and adding `sgn = -sgn` makes the program significantly slower in the MFLOPS calculation. – Axion004 Jan 19 '19 at 20:14
  • And I guess num = 10^8 is probably better, because 10^6 is too small (and so affected by "noise"...) On my computer (+ gfortran), the sgn approach was generally the best, and (-1)**(i+1)/(2.0d0 * i - 1.0d0) was fast when -Ofast was used. – roygvib Jan 19 '19 at 20:51
  • Don't know how your doing your timing or compiler and compiler options. With gfortran -O, I get 18.63 clock ticks per iteration for `val = val + (-1.d0)**(i+1) / (2.d0 * i - 1.d0)`. I get 16.5 clock ticks per iteration for `val = val + (-1.d0)**(i+1) / (2 * i - 1)`. I get 15.8 clock ticks per iteration for `val = val + s / (2 * i - 1); s = - s`. Flipping a sign bit will likely always be faster than computing '(-1.d0)**(i+1). – evets Jan 19 '19 at 20:57
  • That is strange, the timing calculation is done through `MFLOPS = num*2.d0/((T2-T1)*1.d8)`. The variable (number of MFLOPS) is higher when applying `val = val + s/(2.d0*i-1.d0)` and `s=-s`. – Axion004 Jan 19 '19 at 21:10
  • You may want to use the kind `real64` from Fortran's intrinsic module `iso_fortran_env` instead of `double precision`, which is a deprecated feature of Fortran, and possibly non-portable. See here for example: https://stackoverflow.com/questions/50595344/using-iso-fortran-env-to-set-a-functions-kind-value – Scientist Jan 20 '19 at 17:50
  • `real64` is not exactly a replacement for `double precision`. The standard simply says that `double precision`'s precision must be greater than that of default real. If you have some specific *precision* requirement, then use `selected_real_kind`; if you have some specific *storage size* requirement, the use one of the constants in `iso_fortran_env`, like `real64`. See [this](https://stackoverflow.com/a/52800267/2938526). – Rodrigo Rodrigues Jan 21 '19 at 19:49