Fortran 中 Sin 和 Cos 函数的精度

Question

考虑代码，

PROGRAM TRIG_TEST
IMPLICIT NONE

DOUBLE PRECISION, PARAMETER :: PI=4.D0*DATAN(1.0D)

print *, sin(PI/2.0), cos(PI/2.0)

END PROGRAM TRIG_TEST

编译 gfortran 输出，

1.0000000000000000 6.1232339957367660E-017

我知道常见的浮点数问题，但是 sin 函数完全为 1，而 cos 函数不完全为零，这是有原因的吗？

Answer 1

让我们看看为什么 sin 的结果给出 1。在您的代码中

DOUBLE PRECISION, PARAMETER :: PI=4.D0*DATAN(1.0D0)

这是 pi 值的浮点近似值。它可能非常接近实际值，但事实并非如此。另外，sin 函数在计算 sin 时有错误。我们希望这些都是小错误。

我们预计 cos(pi/2) 的值为零。您的浮点计算 cos(PI/2) 与数学答案有大约 6.1232339957367660E-017 的误差。假设您的 sin 计算具有类似的误差幅度。

现在看看epsilon(0d)的值。这是 1d0+epsilon(0d0) 不等于 1d0 的最小数字。在模型中（您报告“~ 2.2e-16”）中的假定误差远小于此数字。

因此，1 是最接近浮点计算实际值的可表示数。

考虑该计划

  use, intrinsic :: iso_fortran_env, only : real128, real64
  implicit none

  real(real64), parameter :: PI=4.D0*ATAN(1._real64)
  real(real128), parameter :: PI_approx = PI   ! Not 4*ATAN(1._real128)

  print *, SIN(PI/2), COS(PI/2)
  print *, SIN(PI_approx/2), COS(PI_approx/2)

end

这会计算（可能）您的 PI/2 的罪，但精度更高（对 pi 使用相同的近似值）。在第二种情况下，我的编译器报告的值与 1 不同，但差异远小于 epsilon(0._real64)。

作为样式点，一般来说最好避免datan，而使用通用的atan。 double precision 也可以替换为适当的种类参数。这些都显示在我上面的程序中。

Answer 2

以下假设double是IEEE 754基本64位二进制格式。三角函数的常见实现不如格式支持的准确。但是，对于这个答案，我们假设它们 return 可能是最准确的结果。

π 不能精确表示在 double 中。最接近的可能值是 884279719003555 / 281474976710656 或 3.141592653589793115997963468544185161590576171875。我们称之为 p.

p/2 的正弦约为 1 − 1.8747•10⁻³³。 double 中可表示的两个值分别为 1 和 0.99999999999999988897769753748434595763683319091796875，大约为 1 − 1.11•10⁻¹⁶。其中较接近的是 1，因此最接近 p/2 正弦值的可表示值恰好是 1。

p/2的余弦约为6.123233995736765886•10⁻¹⁷。 double 中最接近的值是 6.12323399573676603586882014729198302312846062338790031898128063403419218957424163818359375•10^{−17[=135=].[=}

因此您观察到的结果是最接近真实数学值的可能结果。

Fortran 中 Sin 和 Cos 函数的精度

Precision of Sin vs Cos functions in Fortran

floating-point

precision

fortran

gfortran