Borwein在Fortran中计算Pi的算法收敛太快 [英] Borwein’s algorithm for the calculation of Pi in Fortran is converging too fast

问题描述

以下在Fortran中采用四次收敛的 Borwein算法的实现公认计算出Pi，但是收敛太快了.从理论上讲， a 会四次收敛到1/π.因此，在每次迭代中，正确数字的数量将增加四倍.

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

但是在第二次迭代之后，Pi的值不再改变，正如人们看到的前100个数字一样:

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

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

(The last output is the correct value of Pi for comparison.)

实施中是否有错误?我不确定，是否始终保持 real128 的精度.

Is there an error in the implementation? I’m not sure, if real128 precision is always kept.

推荐答案

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

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

这篇关于Borwein在Fortran中计算Pi的算法收敛太快的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！