The following implementation of Borwein’s algorithm with quartic convergence in Fortran admittedly calculates Pi, but converges simply too fast. In theory, a converges quartically to 1/π. On each iteration, the number of correct digits is therefore quadrupled.
! pi.f90
program main
use, intrinsic :: iso_fortran_env, only: real128
implicit none
real(kind=real128), parameter :: CONST_PI = acos(-1._real128)
real(kind=real128) :: pi
integer :: i
do i = 1, 10
pi = borwein(i)
print '("Pi (n = ", i3, "): ", f0.100)', i, pi
end do
print '("Pi:", 11x, f0.100)', CONST_PI
contains
function borwein(n) result(pi)
integer, intent(in) :: n
real(kind=real128) :: pi
real(kind=real128) :: a, y
integer :: i
y = sqrt(2._real128) - 1
a = 2 * (sqrt(2._real128) - 1)**2
do i = 1, n
y = (1 - (1 - y**4)**.25_real128) / (1 + (1 - y**4)**.25_real128)
a = a * (1 + y)**4 - 2**(2 * (i - 1) + 3) * y * (1 + y + y**2)
end do
pi = 1 / a
end function borwein
end program main
But after the second iteration, the value of Pi does not change anymore, as one can see for the first 100 digits:
$ gfortran -o pi pi.f90
$ ./pi
Pi (n = 1): 3.1415926462135422821493444319826910539597974491424025615960346875056357837663334464650688460096716881
Pi (n = 2): 3.1415926535897932384626433832795024122930792206901249618089158148888330110426458929850923595950007439
Pi (n = 3): 3.1415926535897932384626433832795024122930792206901249618089158148888330110426458929850923595950007439
Pi (n = 4): 3.1415926535897932384626433832795024122930792206901249618089158148888330110426458929850923595950007439
Pi (n = 5): 3.1415926535897932384626433832795024122930792206901249618089158148888330110426458929850923595950007439
Pi (n = 6): 3.1415926535897932384626433832795024122930792206901249618089158148888330110426458929850923595950007439
Pi (n = 7): 3.1415926535897932384626433832795024122930792206901249618089158148888330110426458929850923595950007439
Pi (n = 8): 3.1415926535897932384626433832795024122930792206901249618089158148888330110426458929850923595950007439
Pi (n = 9): 3.1415926535897932384626433832795024122930792206901249618089158148888330110426458929850923595950007439
Pi (n = 10): 3.1415926535897932384626433832795024122930792206901249618089158148888330110426458929850923595950007439
Pi: 3.1415926535897932384626433832795027974790680981372955730045043318742967186629755360627314075827598572
(The last output is the correct value of Pi for comparison.)
Is there an error in the implementation? I’m not sure, if real128 precision is always kept.
At the second iteration you have converged to as many digits as real128 can support:
ian@eris:~/work/stack$ cat pi2.f90
Program pi2
Use, intrinsic :: iso_fortran_env, only: real128
Implicit None
Real( real128 ) :: pi
pi = 3.1415926535897932384626433832795024122930792206901249618089158148888330110426458929850923595950007439_real128
Write( *, '( f0.100 )' ) pi
Write( *, '( f0.100 )' ) Nearest( pi, +1.0_real128 )
Write( *, '( f0.100 )' ) Abs( acos(-1._real128) - Nearest( pi, +1.0_real128 ) )
End Program pi2
ian@eris:~/work/stack$ gfortran -std=f2008 -Wall -Wextra -fcheck=all pi2.f90
ian@eris:~/work/stack$ ./a.out
3.1415926535897932384626433832795024122930792206901249618089158148888330110426458929850923595950007439
3.1415926535897932384626433832795027974790680981372955730045043318742967186629755360627314075827598572
.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
ian@eris:~/work/stack$
Thus the further iterations show no change as it is already as accurate as it can be