我正在測試一些非常簡單的等價錯誤,當精度是一個問題,並希望以擴展雙精度執行操作(以便我知道答案將在〜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