用 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 代码以完全双精度运行，在我的机器上几秒钟内计算出 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 为 Fortran 中的 integer(8) :: k 和 unsigned long在 C 中。通过此修改，对于 pi 的 10/11 位，Fortran 代码现在几乎与 C 代码完全相同，并且似乎 'hang' 精度更高。

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

Answer 1

您的 Fortran 代码不正确。

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

中的评估发生整数溢出

就像您使用 selected_real_kind 来找到符合您要求的真实种类一样，您应该使用 selected_int_kind 来找到合适的整数种类。如果我们信任 C 版本，那么迭代次数可能会达到 k 应该有种类 selected_int_kind(11).

这样大的数字

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

Calculating Pi with Leibniz's series in C and Fortran

c

fortran

numeric