我正在测试一些非常简单的等价错误,当精度是一个问题,并希望以扩展双精度执行操作(以便我知道答案将在〜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位数字是相同的,并且期望的双精度舍入错误没有发生。
还应当指出的是,这次尝试前,我写这实际上写另一节目里的价值观x
,y
,并且z
是“硬编码”到20位小数的程序。在此版本的程序中,.gt.
和.lt.
的比较失败,但我无法通过将扩展的双精度值作为双精度变量进行投射来复制相同的故障。
为了进一步“干扰”双精度值并添加舍入误差,我添加了,然后减去pi
从我的双精度变量应该留下剩余的变量与一些双精度舍入误差,但我在最终结果中仍然没有看到。
你的函数rand()是什么? – francescalus
这是一个'Fortran'内在函数。 https://gcc.gnu.org/onlinedocs/gfortran/RAND.html – drjrm3
它并不是一个Fortran内在的。正如你使用的是gcc,那么 - 正如你所说的那样有这个内在 - 请把它放在这个问题上。这是因为并不是所有编译器都有这样的事情(而不是所有的事情)都会做同样的事情,它对于答案是重要的。 – francescalus