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

我正在测试一些非常简单的等价错误,当精度是一个问题,并希望以扩展的双精度执行操作(这样我知道答案将在~19位数内)然后以双精度执行相同的操作(其中在第16位数会出现舍入误差,但不知怎的,我的双精度算术保持了19位精度.

当我在扩展double中执行操作时,然后将数字硬编码到另一个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功能在这里.

如果我们只谈论一种精确类型(例如,取双倍),那么我们可以将机器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,其中只有x8符合值的前16位数x10,但在打印出值时,它显示所有20位数都是同样和预期的双精度舍入误差没有像我预期的那样发生.

还应当指出的是,这次尝试前,我写这实际上写了一个程序的另一个节目里的价值观x,y以及z被"硬编码"到20位小数.在此版本的程序中,按预期进行比较.gt.和.lt.失败,但我无法通过将扩展的双精度值作为双精度变量进行复制来复制相同的失败.

为了进一步"扰乱"双精度值并添加舍入误差,我已经pi从我的双精度变量中添加,然后减去,这应该留下剩余变量具有一些双精度舍入误差,但我仍然没有看到最后的结果.

作为链接状态的gfortran文档,函数结果rand是默认实际值(单精度).这样的值可以由您的每个其他真实类型精确表示.

也就是说,x10=rand()将单个精度值分配给扩展精度变量x10.它完全如此.现在存储的相同值x10被分配给双精度变量x8,但这仍然可以表示为双精度.

单对双精度具有足够的精度,使用double和extended类型的计算返回相同的值.[见本答复末尾的注释.]

如果您希望看到精度损失的真实效果,那么首先使用扩展或双精度值.例如,不使用rand(返回单个精度值),而是使用内在函数random_number

call random_number(x10)

(它具有标准Fortran的优点).与函数不同,它在(几乎)所有的情况下,无论返回值类型的最终用途的价值的,该子程序会给你相应的参数精度.你(希望)会从你的"硬编码"实验中看到很多.

或者,正如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支持内部模块,可以以便携/标准方式使用.

您可以将此标志视为进一步的实验:如果它被提升,那么您可以期望看到两个精度之间的差异.