fortran - 这些 double 值如何精确到 20 位小数？

当精度是一个问题时，我正在测试一些非常简单的等价错误，并希望以扩展 double 执行操作(以便我知道大约 19 位数字的答案)，然后以 double 执行相同的操作(其中第 16 位会有舍入误差)，但不知何故我的 double 算术保持了 19 位的精度。

当我在扩展 double 中执行操作，然后将数字硬编码到另一个 Fortran 例程中时，我得到了预期的错误，但是当我将扩展 double 变量分配给这里的 double 变量时会发生什么奇怪的事情吗？

program code_gen
    implicit none
    integer, parameter :: Edp = selected_real_kind(17)
    integer, parameter :: dp = selected_real_kind(8)
    real(kind=Edp) :: alpha10, x10, y10, z10
    real(kind=dp) :: alpha8, x8, y8, z8

    real(kind = dp) :: pi_dp = 3.1415926535897932384626433832795028841971693993751058209749445

    integer :: iter
    integer :: niters = 10

    print*, 'tiny(x10) = ', tiny(x10)
    print*, 'tiny(x8)  = ', tiny(x8)
    print*, 'epsilon(x10) = ', epsilon(x10)
    print*, 'epsilon(x8)  = ', epsilon(x8)

    do iter = 1,niters
        x10 = rand()
        y10 = rand()
        z10 = rand()
        alpha10 = x10*(y10+z10)

        x8 = x10
        x8 = x8 - pi_dp
        x8 = x8 + pi_dp
        y8 = y10
        y8 = y8 - pi_dp
        y8 = y8 + pi_dp
        z8 = z10
        z8 = z8 - pi_dp
        z8 = z8 + pi_dp
        alpha8 = alpha10

        write(*, '(a, es30.20)') 'alpha8 .... ', x8*(y8+z8)
        write(*, '(a, es30.20)') 'alpha10 ... ', alpha10

        if( alpha8 .gt. x8*(y8+z8) ) then
            write(*, '(a)') 'ERROR(.gt.)'
        elseif( alpha8 .lt. x8*(y8+z8) ) then
            write(*, '(a)') 'ERROR(.lt.)'
        endif
    enddo
end program code_gen

其中 rand() 是找到的 gfortran 函数 here 。

如果我们只讨论一种精度类型(例如，double)，那么我们可以将机器 epsilon 表示为 E16 ，它大约是 2.22E-16 。如果我们简单地将两个实数 x+y 相加，那么得到的机器表达数是 (x+y)*(1+d1) 其中 abs(d1) < E16 。同样，如果我们然后将该数字乘以 z ，结果值实际上是 (z*((x+y)*(1+d1))*(1+d2)) ，它几乎是 (z*(x+y)*(1+d1+d2)) ，其中 abs(d1+d2) < 2*E16 。如果我们现在转向扩展 double ，那么唯一改变的是 E16 变成 E20 并且具有大约 1.08E-19 的值。

我希望以扩展 double 执行分析，以便我可以比较两个应该相等的数字，但有时会显示舍入误差会导致比较失败。通过分配 x8=x10 ，我希望创建扩展 double 值 x10 的 double “版本”，其中只有 x8 的前~16 位数字符合 x10 的值，但是在打印出这些值时，它表明所有20 位数字是相同的，预期的 double 舍入误差并没有像我预期的那样发生。

还应该注意的是，在此尝试之前，我编写了一个程序，该程序实际上编写了另一个程序，其中 x 、 y 和 z 的值被“硬编码”到小数点后 20 位。在此版本的程序中，.gt. 和 .lt. 的比较按预期失败，但我无法通过将扩展 double 值转换为 double 变量来复制相同的失败。

为了进一步“扰乱” double 值并添加舍入误差，我已经添加，然后从我的 double 变量中减去了 pi，这应该使剩余的变量具有一些 double 舍入误差，但我仍然没有看到在最终结果中。

最佳答案

正如您链接的 gfortran 文档所述，rand 的函数结果是默认的实数值(单精度)。这样的值可以由您的每个其他实数类型精确表示。

也就是说， x10=rand() 将单精度值分配给扩展精度变量 x10 。它正是这样做的。现在存储在 x10 中的相同值被分配给 double 变量 x8 ，但这仍然可以完全表示为 double 。

single-as-double 有足够的精度，使用 double 和扩展类型的计算返回相同的值。 [请参阅本答案末尾的注释。]

如果您希望看到精度损失的实际影响，请从使用扩展或 double 值开始。例如，不要使用 rand(返回单精度值)，而是使用内在的 random_number

call random_number(x10)

(这具有成为标准 Fortran 的优势)。与在(几乎)所有情况下都返回值类型而不管值的最终用途如何的函数不同，此子例程将为您提供与参数相对应的精度。您将(希望)从“硬编码”实验中看到更多内容。

或者，正如 agentp 所评论的，从 double 值开始可能更直观

call random_number(x8); x10=x8   ! x8 and x10 have the precision of double precision
call random_number(y8); y10=y8
call random_number(z8); z10=z8

并从该起点执行计算:然后这些额外的位将开始显示。

总之，当您执行 x8=x10 时，您将获得与 x8 对应的 x10 的前几位，但其中许多位以及 x10 中的许多位都为零。

当涉及到 pi_dp 扰动时，您再次将单精度(这次是文字常量)值分配给 double 变量。仅仅拥有所有这些数字并不能使它成为默认的真实文字。您可以使用 _Edp 后缀指定不同类型的文字，如其他答案中所述。

最后，人们还必须担心编译器对 regards to optimization 做了什么。

我的论点是，从单精度值开始，执行的计算可以精确地表示为 double 和扩展精度(具有相同的值)。对于其他计算，或从设置更多位的起点或表示(例如，在某些系统或其他编译器上，类型为 selected_real_kind(17) 的数字类型可能具有完全不同的特征，例如不同的基数)，它们不必是案件。

虽然这主要是基于猜测并希望它能解释观察结果。幸运的是，有一些方法可以测试这个想法。当我们谈论 IEEE 算术时，我们可以考虑不精确标志。如果在计算过程中没有提出该标志，我们会很高兴。

使用 gfortran 有一个编译选项 -ffpe=inexact ，它会发出不精确的标志信号。使用 gfortran 5.0，支持内部模块 ieee_exceptions，可以以便携式/标准方式使用。

您可以考虑使用此标志进行进一步实验:如果它被提升，那么您可以期望看到两个精度之间的差异。

关于fortran - 这些 double 值如何精确到 20 位小数？，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/34639858/