Borwein’s algorithm for the calculation of Pi in Fortran is converging too fast

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