How to resolve the algorithm Arithmetic-geometric mean/Calculate Pi step by step in the Fortran programming language
How to resolve the algorithm Arithmetic-geometric mean/Calculate Pi step by step in the Fortran programming language
Table of Contents
Problem Statement
Almkvist Berndt 1988 begins with an investigation of why the agm is such an efficient algorithm, and proves that it converges quadratically. This is an efficient method to calculate
π
{\displaystyle \pi }
. With the same notations used in Arithmetic-geometric mean, we can summarize the paper by writing:
π
4
a g m
( 1 , 1
/
2
)
2
1 −
∑
n
1
∞
2
n + 1
(
a
n
2
−
g
n
2
)
{\displaystyle \pi ={\frac {4;\mathrm {agm} (1,1/{\sqrt {2}})^{2}}{1-\sum \limits {n=1}^{\infty }2^{n+1}(a{n}^{2}-g_{n}^{2})}}}
This allows you to make the approximation, for any large N:
π ≈
4
a
N
2
1 −
∑
k
1
N
2
k + 1
(
a
k
2
−
g
k
2
)
{\displaystyle \pi \approx {\frac {4;a_{N}^{2}}{1-\sum \limits {k=1}^{N}2^{k+1}(a{k}^{2}-g_{k}^{2})}}}
The purpose of this task is to demonstrate how to use this approximation in order to compute a large number of decimals of
π
{\displaystyle \pi }
.
Let's start with the solution:
Step by Step solution about How to resolve the algorithm Arithmetic-geometric mean/Calculate Pi step by step in the Fortran programming language
Source code in the fortran programming language
program CalcPi
! Use real128 numbers: (append '_rf')
use iso_fortran_env, only: rf => real128
implicit none
real(rf) :: a,g,s,old_pi,new_pi
real(rf) :: a1,g1,s1
integer :: k,k1,i
old_pi = 0.0_rf;
a = 1.0_rf; g = 1.0_rf/sqrt(2.0_rf); s = 0.0_rf; k = 0
do i=1,100
call approx_pi_step(a,g,s,k,a1,g1,s1,k1)
new_pi = 4.0_rf * (a1**2.0_rf) / (1.0_rf - s1)
if (abs(new_pi - old_pi).lt.(2.0_rf*epsilon(new_pi))) then
! If the difference between the newly and previously
! calculated pi is negligible, stop the calculations
exit
end if
write(*,*) 'Iteration:',k1,' Diff:',abs(new_pi - old_pi),' Pi:',new_pi
old_pi = new_pi
a = a1; g = g1; s = s1; k = k1
end do
contains
subroutine approx_pi_step(x,y,z,n,a,g,s,k)
real(rf), intent(in) :: x,y,z
integer, intent(in) :: n
real(rf), intent(out) :: a,g,s
integer, intent(out) :: k
a = 0.5_rf*(x+y)
g = sqrt(x*y)
k = n + 1
s = z + (2.0_rf)**(real(k)+1.0_rf) * (a**(2.0_rf) - g**(2.0_rf))
end subroutine
end program CalcPi
You may also check:How to resolve the algorithm Hofstadter Figure-Figure sequences step by step in the Racket programming language
You may also check:How to resolve the algorithm Sorting algorithms/Shell sort step by step in the PureBasic programming language
You may also check:How to resolve the algorithm Hello world/Graphical step by step in the TI-83 BASIC programming language
You may also check:How to resolve the algorithm Rosetta Code/Find bare lang tags step by step in the Rust programming language
You may also check:How to resolve the algorithm Compare length of two strings step by step in the EasyLang programming language