2016-01-06 88 views
3

我正在测试一些非常简单的等价错误,当精度是一个问题,并希望以扩展双精度执行操作(以便我知道答案将在〜19位数),然后以双精度执行相同的操作(第16位数字会出现舍入误差),但不知怎的,我的双精度算法保持了19位数的精度。这些双精度值如何精确到20位小数?

当我在extended 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功能发现here

如果我们只讲一种精度类型(例如双倍),那么我们可以将机器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位数字是相同的,并且期望的双精度舍入错误没有发生。

还应当指出的是,这次尝试前,我写这实际上写另一节目里的价值观xy,并且z是“硬编码”到20位小数的程序。在此版本的程序中,.gt..lt.的比较失败,但我无法通过将扩展的双精度值作为双精度变量进行投射来复制相同的故障。

为了进一步“干扰”双精度值并添加舍入误差,我添加了,然后减去pi从我的双精度变量应该留下剩余的变量与一些双精度舍入误差,但我在最终结果中仍然没有看到。

+0

你的函数rand()是什么? – francescalus

+0

这是一个'Fortran'内在函数。 https://gcc.gnu.org/onlinedocs/gfortran/RAND.html – drjrm3

+0

它并不是一个Fortran内在的。正如你使用的是gcc,那么 - 正如你所说的那样有这个内在 - 请把它放在这个问题上。这是因为并不是所有编译器都有这样的事情(而不是所有的事情)都会做同样的事情,它对于答案是重要的。 – francescalus

回答

6

作为链接状态的gfortran文档,rand的函数结果是默认实际值(单精度)。这样的价值可以完全由您的其他实际类型来表示。

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

使用double和extended类型的计算返回相同的值时,double-a-double中有足够的精度。 [请参阅本答案末尾的注释。]

如果您希望看到精确度损失的实际影响,请使用扩展或双精度值开始。例如,而不是使用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你得到相应于那些x10x8前几位,但许多那些位和那些遵循x10都是零。

当涉及到pi_dp扰动时,您再次将一个精度(这次是一个文字常量)赋值给一个双精度变量。只有拥有所有这些数字并不会使其成为默认真实文本以外的任何其他数字。正如其他答案中所述,您可以使用_Edp后缀指定不同类型的文字。

最后,人们还不得不担心编译器用regards to optimization做什么。


我的论点是,从单精度值开始,所执行的计算可精确地以双精度和扩展精度(具有相同的值)表示。对于其他计算,或者从具有更多位集的起点或表示(例如,在某些系统或其他编译器中,类型为selected_real_kind(17)的数字类型可能具有完全不同的特性,例如不同的基数),而这些特性不必是案件。

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

与gfortran有编译选项-ffpe=inexact这将使不准确的标志信号。使用gfortran 5.0,支持固有模块ieee_exceptions,可用于便携/标准方式。

你可以考虑这个标志进行进一步的实验:如果它被提出,那么你可以期望看到两个精度的差异。

+0

现象解释,你100%正确。 – drjrm3

+0

我认为重点是从相同的双表示值开始,并显示扩展精度影响计算结果。 (并抛出一些划分,我想它会..)我会做'x8 = rand(); x10 = x8“,那么你就会知道这是事实。 – agentp

+0

@agentp对第一部分是合理的解释(我在答案中增加了一些内容)。对于第二部分,您的'x8'和'x10'仍然会以相同的(单精度)值开始 - 或者我没有正确读取您的内容? – francescalus