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

Question

我正在尝试比较 C 和 Fortran 代码的性能。使用 Leibniz's series 计算圆周率, 我有以下 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 代码:

#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;
} 

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

Fortran 代码以完全 double 运行，在我的机器上几秒钟内计算出 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 代码中的循环，使收敛变得单调:

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

虽然这在一定程度上加快了较低精度的收敛速度，但它实际上使 C 程序现在甚至在 precision=8 时基本上挂起(计算时间超过 3 分钟)。

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

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

Answer 1

您的 Fortran 代码不正确。

您可能会使用默认的 32 位整数并使用 HUGE(k) 您会看到 k 可以采用的最大整数值是 2147483647。在这种情况下你会有整数溢出发生在迭代计数和(在此之前)其在 real(2*k+1,dp) 中的评估中。

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