Fortran 中 Sin 与 Cos 函数的精度答案

【问题标题】：Precision of Sin vs Cos functions in FortranFortran 中 Sin 与 Cos 函数的精度
【发布时间】：2018-01-26 08:46:24
【问题描述】：

考虑代码，

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 函数不完全相同为零？

【问题讨论】：

您的编译器的episilon(1d0) 的值是多少？
值是~2.2e-16？
听起来很有道理。假设sin 和cos 都精确到7E-017 之内。这将为您提供cos 的答案。 1-7E-017 等于 1. [你说你理解“通常的浮点问题”。您是否需要对这方面进行解释，或者您是否称其为轻微疏忽？]
您可以通过使用 PI = ACOS(-1.d0) 来避免 PI 定义中的乘法。像三角函数这样的内在函数不需要双精度特定版本。 ATAN 就足够了，就像您使用 sin 和 cos 而不是 dsin 和 dcos 一样。 Fortran 将根据参数链接正确的版本。同样要正式正确，您的文字应该是双精度的：print*, sin(PI/2.d0), cos(PI/2.d0)。

标签： floating-point fortran precision gfortran

【解决方案1】：

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

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

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

p/2 的余弦约为 6.123233995736765886•10^-17。 double 中可表示的最接近的值是 6.12323399573676603586882014729198302312846062338790031898128063403419218957424163818359375•10

-17。

所以你观察到的结果是最接近真实数学值的结果。

【讨论】：

为了支持这个答案，有兴趣的读者可以使用ieee_selected_real_kind() 为给定的所需特性选择一种类型参数。

【解决方案2】：

让我们看看为什么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 也可以替换为适当的种类参数。这些显示在我上面的程序中。

【讨论】：