在 C 和 Fortran 中使用 Leibniz 级数计算 Pi答案

【问题标题】：Calculating Pi with Leibniz's series in C and Fortran在 C 和 Fortran 中使用 Leibniz 级数计算 Pi
【发布时间】：2019-12-11 04:19:18
【问题描述】：

我正在尝试比较 C 和 Fortran 代码的性能。对于使用Leibniz's series 计算 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 代码：

#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.27 秒/1581043254 次迭代中收敛到 3.141592653，而 C 代码需要 12.9 秒/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 和 C 中的 unsigned long。通过此修改，Fortran 代码现在执行几乎与 10/11 位 pi 的 C 代码完全相同，并且似乎以更高的精度“挂起”。

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

【问题讨论】：

对C代码的建议（我“不会说”Fortran）：展开for循环，摆脱sgn：for (k = 1; fabs(sum - correct) >= delta; k += 2) { sum += 1 / (2 * k - 1); sum -= 1 / (2 * k + 1); }
如果你说 k 在 Fortran 版本中达到 1581043254，那么 2*k+1 很可能会遇到整数溢出。你应该考虑selected_int_kind()。就像你做的那样selected_real_kind()。
在我之前的评论中，将1s 设为双倍：sum += 1.0 / (2 * k - 1); OOPS
除了 C：size_t 对于需要数十亿范围的计数器来说可能是错误的类型。 unsigned long long、uint_least64_t 或 uint_fast64_t 之一会更好。
您似乎没有比较 C 和 Fortran 的性能。相反，您正在测试您在 C 和 Fortran 中编程的能力，因为这 2 个代码并不等效。首先，delta 的表达式的 RHS 以单精度计算。其次，在 C 中，您正在翻转符号sgn = -sgn，而在 Fortran 中，您计算 (-1)**k。 Fortran 编译是否识别习语并优化 (-1)**k 不是由语言指定的。第三，在 Fortran 代码中，主计算循环中有 if (explicit)，而 C 版本没有此测试。你要害死我了。

标签： c fortran numeric

【解决方案1】：

您的 Fortran 代码不正确。

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

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

【讨论】：