本文的理论基础来自是电脑杂志1996年第10期，作者郭继展发表的一篇文章，作者提出一个公式：

{\color{Blue} \pi =16\arctan \left ( \frac{1}{5} \right )-4\arctan \left ( \frac{1}{239} \right )}

在展开成两个级数之和，然后整理得到：

PI=16x(1/5-1/(5^3/3)+1/(5^5/5)-1/(5^7/7)+...)-4x(1/239-1/(239^3/3)+1/(239^5/5)-1/(239^7/7)+...)

=4x(4x5/25-239/57121)/1-4x(4x5/25^2-239/57121^2)/3+4x(4x5/25^3-239/57121^3)/5-

Module PI_Calc_Mod

Implicit None

Integer , public , parameter :: PI_KIND = KIND(1)

Private :: GenPI , Rot

contains

Subroutine CalcPI( nBit , myPI )

Integer , Intent( IN ) :: nBit

integer(Kind=PI_KIND) :: myPI(nBit),bitA(nBit),bitB(nBit),BitAs(nBit),BitBs(nBit),BitPIA(nBit),BitPIB(nBit)

integer :: i

myPI(:)=0

bitA(:)=0

bitB(:)=0

BitAs(:)=0

BitBs(:)=0

BitPIA(:)=0

BitPIB(:)=0

BitAs(2)=3

BitAs(3)=2

BitBs(2)=4

call Rot( nBit , BitBs , bitB , BitBs , 239 )

call GenPI( nBit , myPI , BitBs , BitAs )

Do i = 1 , nBit

call Rot( nBit , BitAs , bitA , BitAs , 25 )

call Rot( nBit , BitAs , bitA , BitPIA , 2*i+1 )

call Rot( nBit , BitBs , bitB , BitBs , 57121 )

call Rot( nBit , BitBs , bitB , BitPIB , 2*i+1 )

if( mod(i,2) == 1 ) then

call GenPI( nBit , myPI , BitPIA , BitPIB )

else

call GenPI( nBit , myPI , BitPIB , BitPIA )

end if

End Do

End Subroutine CalcPI

Subroutine Rot( nBit , nA , nB , nC , nM )

Integer , Intent( IN ) :: nBit , nM

Integer(Kind=PI_KIND) , Intent( INOUT ) :: nA(nBit) , nB(nBit) , nC(nBit)

integer :: i , j

Do i = 2 , nBit

j = mod( nB(i-1) , nM ) * 10 + nA( i )

nC(i) = j / nM

nB(i) = mod( j , nM )

End Do

End Subroutine Rot

Subroutine GenPI( nBit , nPI , npA , npB )

Integer , Intent( IN ) :: nBit

Integer(Kind=PI_KIND) , Intent( INOUT ) :: nPI(nBit) , npA(nBit) , npB(nBit)

integer :: i

nPI(2:nBit)=nPI(2:nBit)+npB(2:nBit)

Do i = nBit , 2 , -1

if( nPI(i) >= npA(i) ) then

nPI(i) = nPI(i) - npA(i)

else

nPI(i) = nPI(i) + 10 - npA(i)

nPI(i-1) = nPI(i-1) - 1

end if

nPI(i-1) = nPI(i-1) + nPI(i) / 10

nPI(i) = mod( nPI(i) , 10 )

End Do

End Subroutine GenPI

End Module PI_Calc_Mod

Program www_fcode_cn

use PI_Calc_Mod

Implicit None

Integer , parameter :: N = 100

integer(Kind=PI_KIND) :: myPI(N)

call CalcPI( N , myPI )

write(*,'(50i1)') myPI

End Program www_fcode_cn