这些双精度值如何精确到 20 位小数？

Question

当精度成为问题时，我正在测试一些非常简单的等价错误，并希望以扩展双精度执行操作（以便我知道答案是 ~19 位数字），然后执行相同的操作双精度（第 16 位会有舍入误差），但不知何故我的双精度算法保持了 19 位的精度。

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

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。如果我们现在转向扩展双精度，那么唯一改变的是 E16 变为 E20 并且值约为 1.08E-19.

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

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

为了进一步 'perturb' 双精度值并添加舍入误差，我添加了，然后从我的双精度变量中减去 pi，这应该使剩余的变量具有一些双精度精度舍入错误，但我仍然没有在最终结果中看到它。

Answer 1

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

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

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

如果您希望看到精度损失的实际效果，请从使用扩展精度值或双精度值开始。例如，与其使用 rand（return 单个精度值），不如使用固有的 random_number

call random_number(x10)

（具有标准 Fortran 语言的优势）。与在（几乎）所有情况下 return 都是值类型而不管值的最终用途如何的函数不同，此子例程将为您提供与参数相对应的精度。您会（希望）从 "hard-coded" 实验中看到很多东西。

或者，正如 agentp 评论的那样，从双精度值开始可能更直观

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 扰动时，您再次将单精度（这次是文字常量）值分配给双精度变量。仅仅拥有所有这些数字并不能使它成为默认的真实文字。您可以使用 _Edp 后缀指定不同类型的文字，如其他答案中所述。

最后，还要担心编译器对做了什么。

我的论点是，从单精度值开始，执行的计算可以精确地表示为双精度和扩展精度（具有相同的值）。对于其他计算，或从具有更多位集或表示的起点开始（例如，在某些系统或其他编译器中，种类 selected_real_kind(17) 的数字类型可能具有完全不同的特征，例如不同的基数）需要不会吧。

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

对于 gfortran，有一个编译选项 -ffpe=inexact，它会生成不精确的标志信号。 gfortran 5.0 支持内部模块 ieee_exceptions，可以以 portable/standard 方式使用。

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

这些双精度值如何精确到 20 位小数？

How are these double precision values accurate to 20 decimals?

precision

fortran