2016-01-06 136 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