Borwein 在 Fortran 中计算 Pi 的算法收敛太快

Question

在 Fortran 中使用四次收敛的 Borwein’s algorithm 的以下实现诚然可以计算 Pi，但收敛得太快了。理论上，a 四次收敛于 1/π。因此，在每次迭代中，正确数字的数量增加了四倍。

! 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

但是在第二次迭代之后，Pi 的值不再改变，正如前 100 位数字所见：

$ 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

(最后输出的是Pi的正确值，用于比较。)

执行中是否有错误？我不确定，是否始终保持 real128 精度。

Answer 1

在第二次迭代中，您已经收敛到 real128 可以支持的位数：

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$

因此，进一步的迭代显示没有变化，因为它已经尽可能准确

Borwein 在 Fortran 中计算 Pi 的算法收敛太快

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

fortran

pi