用C和Fortran中的Leibniz系列计算Pi

Question

我正在尝试比较C和Fortran代码的性能。为了使用Leibniz的系列计算pi ，我获得了以下Fortran代码

program pi_leibniz
implicit none

    integer, parameter :: dp=selected_real_kind(15,307)
    integer :: k=0, precision=9
    real(dp), parameter :: correct = 0.7853981633974483d0, eps = epsilon(real(1,dp)) 
    real(dp) :: sum = 0.0, delta
    character(8) :: fmt
    logical, parameter :: explicit = .false.
    real :: start, finish

    delta = 10.**(-precision-1)*0.25
    if (delta<eps) then
        delta=eps
        precision=14
        print *, "Precision specified too high, reverting to double precision (14 digits)"
    endif

    write(fmt,'(A,I0,A,I0,A)') '(f',precision+2,'.',precision,')'

    call cpu_time(start)

    do
        sum = sum + real((-1)**k,dp)/real(2*k+1,dp)
        k = k+1
        if (abs(sum-correct)<delta) exit        
        if (explicit) print fmt, 4.*sum 
    enddo

    call cpu_time(finish)

    print fmt, 4.*sum
    print '(A,I0,A,I0,A)', "converged in ", k, " iterations with ", precision, " digits precision"
    print '(g0,a)', finish-start," s"

end program pi_leibniz

和几乎相同的C代码：

program pi_leibniz
implicit none

    integer, parameter :: dp=selected_real_kind(15,307)
    integer :: k=0, precision=9
    real(dp), parameter :: correct = 0.7853981633974483d0, eps = epsilon(real(1,dp)) 
    real(dp) :: sum = 0.0, delta
    character(8) :: fmt
    logical, parameter :: explicit = .false.
    real :: start, finish

    delta = 10.**(-precision-1)*0.25
    if (delta<eps) then
        delta=eps
        precision=14
        print *, "Precision specified too high, reverting to double precision (14 digits)"
    endif

    write(fmt,'(A,I0,A,I0,A)') '(f',precision+2,'.',precision,')'

    call cpu_time(start)

    do
        sum = sum + real((-1)**k,dp)/real(2*k+1,dp)
        k = k+1
        if (abs(sum-correct)<delta) exit        
        if (explicit) print fmt, 4.*sum 
    enddo

    call cpu_time(finish)

    print fmt, 4.*sum
    print '(A,I0,A,I0,A)', "converged in ", k, " iterations with ", precision, " digits precision"
    print '(g0,a)', finish-start," s"

end program pi_leibniz

我使用GNU编译器和-O2选项进行编译。编辑：64位。

Fortran代码可以愉快地运行到全双精度，可以在我的机器上几秒钟内计算出pi的前15位。C代码的性能甚至比Fortran快一点，最多可以达到8个小数位，并且在相同的迭代次数中收敛到相同的数字。但是，随着precision=9 Fortran代码在2.27s / 1581043254迭代中收敛到3.141592653，而C代码需要12.9s / 9858058108迭代（〜6x）并且最后一位被关闭1。以更高的精度，Fortran的时间为同样的顺序，而C则需要2分钟才能计算出pi的前11位数字。

差异的原因可能是什么？如何避免使C代码变慢的原因？

编辑：我按照@pmg的建议做了，并更改了C代码中的循环，使收敛单调：

#include <stdio.h>
#include <time.h>
#include <float.h>
#include <math.h>


int main(void){
    int precision=9;
    size_t k=0;
    const double correct=0.7853981633974483;
    double sum=0.0, delta = 0.25*pow(10.0,-(precision+1));
    clock_t start,finish;

    double sgn = 1.0;

    if (delta < DBL_EPSILON){
        delta = DBL_EPSILON;
        precision = 14;
        printf("Precision specified too high, reverting to double precision (14 digits)\n");
    }

    start = clock();

    for(k=0; fabs(sum-correct) >= delta; k++, sgn=-sgn)
        sum += sgn/(2*k+1);

    finish = clock();

    printf("%.*f\n",precision,4*sum);
    printf("converged in %zu iterations with %d digits precision\n",k,precision);
    printf("%f s\n",(finish-start)/(double)CLOCKS_PER_SEC);

    return 0;
} 

尽管这会以较低的精度加快收敛速度​​，但实际上甚至使C程序实际上precision=8现在就挂起了（计算需要3分钟以上的时间）。

编辑2：由于计算precision>8结果会导致整数溢出，因此似乎正确的声明方式k与integer(8) :: k在Fortran和unsigned longC中相同。通过此修改，Fortran代码现在的性能几乎与pi的10/11位的C代码完全一样。并且似乎以更高的精度“挂起”。

那么，为什么使用本质上不正确的方法仍然可以产生正确的结果，并且花相同的时间来计算它是pi的10位还是15位？只是为了好玩，它花了1611454902迭代才能“收敛”到3.14159265358979，恰好是pi到小数点后14位。

Answer 1

您的 Fortran 代码不正确。

您可能使用 32 位的默认整数，并且使用时HUGE(k)您会看到最大整数值为k2147483647。在这种情况下，您将在迭代计数以及（在此之前）在real(2*k+1,dp).

就像您用来selected_real_kind查找符合您要求的实数类型一样，您也可以使用 shouldselected_int_kind来查找合适的整数类型。如果我们信任 C 版本，那么迭代计数可能会达到如此大的数字，k应该有 kind selected_int_kind(11)。