浮点数学精度，c++ vs fortran

Question

我试图在 C++ 和 Fortran 中实现一个递归函数，它计算第 n 个勒让德多项式在 x 处的值。在 Fortran 中，我有

recursive function legendre(n, x) result(p)
  integer, intent(in) :: n
  real(8), intent(in) :: x
  real(8) :: p

  if(n == 0) then
    p = 1.0
  else if (n == 1) then
    p = x
  else
    p = (2.0*real(n,dp)-1.0)*x*legendre(n-1,x)-(real(n,dp)-1.0)*legendre(n-2,x)
    p = p / real(n,dp)
  end if 
end function legendre

然后在 C++ 中我有

double legendre(int n, double x) {
  double p;
  if(n == 0) return 1.0;
  else if(n == 1) return x;
  else {
    p = (2.0*(double)n - 1.0)*x*legendre(n-1,x)-((double)n - 1.0)*legendre(n-2,x);
    p /= (double)n;
    return p;
  }
}

这两个函数对我来说似乎完全一样，都使用双精度，但是 Fortran 函数的结果与 C++ 的结果有很大不同。例如，根据 WolframAlpha，legendre(7,-0.2345) = 0.28876207107499049178814404296875。上面的两个代码，在没有优化的情况下编译时会产生

Fortran：0.28876206852410113

C++：0.28876207107499052285

我知道由于浮点运算，答案不应该相同，但是这里的双精度值差异对我来说似乎有些大。Fortran 值与其他两个值相差这么远的原因是什么？

Answer 1

虽然变量在FORTRAN函数定义为双精度（8个字节），常量已指定为默认（单精度，4字节）值。

根据这个讨论，这意味着算术被执行到单精度精度：

即使将结果分配给的变量定义为 DP，Fortran 标准也要求使用 SP 对常量进行算术运算。那是因为您使用的是默认实数常量，因为常量末尾没有任何类型参数。按照规则，默认实数是 SP。

并且，在同一个讨论中进一步：

...从 1991 年 6 月发布的 Fortran 90 开始，禁止这种将 SP 常量“提升”到 DP 的做法。

因此，为了强制使用双精度数学，请将常量指定为 DP：而不是，例如1.0，指定1.0D0（对于其他人，依此类推）。