2012-12-24 49 views
0

我在程序中使用libresample。經過一段時間(大約50分鐘)後,它在一個工作站的lib函數lrsFilterUD()中崩潰。C乘法或加法浮點結果NaN

float lrsFilterUD(float Imp[], /* impulse response */ 
       float ImpD[], /* impulse response deltas */ 
       UWORD Nwing, /* len of one wing of filter */ 
       BOOL Interp, /* Interpolate coefs using deltas? */ 
       float *Xp, /* Current sample */ 
       double Ph, /* Phase */ 
       int Inc, /* increment (1 for right wing or -1 for left) */ 
       double dhb) 
{ 
    float a; 
    float *Hp, *Hdp, *End; 
    float v, t; 
    double Ho; 

    v = 0.0; /* The output value */ 
    Ho = Ph*dhb; 
    End = &Imp[Nwing]; 
    if (Inc == 1)  /* If doing right wing...    */ 
    {      /* ...drop extra coeff, so when Ph is */ 
     End--;   /* 0.5, we don't do too many mult's */ 
     if (Ph == 0)  /* If the phase is zero...   */ 
     Ho += dhb;  /* ...then we've already skipped the */ 
    }       /* first sample, so we must also */ 
         /* skip ahead in Imp[] and ImpD[] */ 

    if (Interp) 
     while ((Hp = &Imp[(int)Ho]) < End) { 
     t = *Hp;  /* Get IR sample */ 
     Hdp = &ImpD[(int)Ho]; /* get interp bits from diff table*/ 
     a = Ho - floor(Ho);  /* a is logically between 0 and 1 */ 
     t += (*Hdp)*a; /* t is now interp'd filter coeff */ 
     t *= *Xp;  /* Mult coeff by input sample */ 
     v += t;   /* The filter output */ 
     Ho += dhb;  /* IR step */ 
     Xp += Inc;  /* Input signal step. NO CHECK ON BOUNDS */ 
     } 
    else 
     while ((Hp = &Imp[(int)Ho]) < End) { 
     dprintf("while begin: Hp = %p, *Hp = %a, (int)Ho = %d, Imp[(int)Ho] = %a, &Imp[(int)Ho] = %p", Hp, *Hp, (int)Ho, Imp[(int)Ho], &Imp[(int)Ho]); 
     t = *Hp;  /* Get IR sample */ 
     dprintf("before t = %a, *Xp = %a, Xp = %p", t, *Xp, Xp); 
     t *= *Xp;  /* Mult coeff by input sample */ 
     dprintf("after2 t = %a, v = %a", t, v); 
     v += t;   /* The filter output */ 
     dprintf("v = %a", v); 
     Ho += dhb;  /* IR step */ 
     Xp += Inc;  /* Input signal step. NO CHECK ON BOUNDS */ 
     } 

    return v; 
} 

我登錄,* XP,XP之前和相乘後的T值:

while begin: Hp = 0xaf5daa8, *Hp = -0.009034, (int)Ho = 16384, Imp[(int)Ho] = -0.009034, &Imp[(int)Ho] = 0xaf5daa8 
before multiplication t = -0.009034, *Xp = 0.000000, Xp = 0xaebe9b8 
after multiplication t = nan 

此代碼運行多次,也有相同的T和XP值墜毀前:

before multiplication t = -0.009034, *Xp = 0.000000, Xp = 0xaebe9c8 
after multiplication t = -0.000000, v = 282.423676 

或添加另一種情況:

before addition t = -460.799988, v = 0.000000 
after addition v = nan 

什麼可能導致南?這是在Linux上用gcc 4.1.2編譯的。

更新:將變量打印爲%a。結果:

//t = 0x1.2806bap+2 
//Hp = 0xb3bb870 
t = *Hp; 
//t = nan 

更新2:沒有這樣的問題,如果代碼被編譯ICPC。那麼是否有編譯器的具體問題?

+0

這種代碼格式很煩人。 – 2012-12-24 13:08:36

+0

請發佈您從中獲取日誌條目的實際代碼。另外,確保'Ho'沒有超出範圍。 (什麼類型是「Ho」?) –

+0

通過記錄添加了整個函數的代碼。 – Alex

回答

5

顯然,-0.009034•0.000000不應該產生NaN。因此,問題中提供的代碼和數據不是實際計算的準確表示,或者計算實現有缺陷。

如果我們假設硬件和基本運算的實現是沒有缺陷的,那麼一些可能性,調查包括:

  • t*Xp日誌記錄失敗立即乘法之前登錄的t*Xp正確的價值觀或在乘法後立即得到正確的值t
  • t*Xp的值的顯示不正確。例如,用於顯示*Xp的格式顯示「0.000000」,儘管*Xp具有其他值,例如NaN。
  • Xp指向某處不恰當,導致*Xp不可靠(例如,由外部操作更改)。
  • 顯示的代碼不準確。例如,它來自舊的來源,已被更改,或者它是新來源,但是之前編譯的代碼正在執行。

注:使用浮點對象調試,你應該打印用的格式,如「%F」,尤其不能與數字編號的默認值。您應該使用「%a」打印,該打印使用十六進制表示法輸出浮點值的確切值。在許多情況下,您也可以使用「%.99g」,前提是您的C實現提供了將浮點值轉換爲小數的良好轉換。

+0

將變量值添加爲要發佈的「%a」。 – Alex

+0

我今天遇到的另一個原因是這個警告(不是海報代碼,而是我自己的):「隱式聲明函數'fabs'」,因爲我忘記了包含math.h –

1

Wiki,有三種操作,可以返回NaN的如下:

1. Operations with a NaN as at least one operand. 
2. Indeterminate forms 
     The divisions 0/0 and ±∞/±∞ 
     The multiplications 0×±∞ and ±∞×0 
     The additions ∞ + (−∞), (−∞) + ∞ and equivalent subtractions 
     The standard has alternative functions for powers: 
     The standard pow function and the integer exponent pown function define 0pow(0), 1pow(∞), 
     and ∞pow(0) as 1. 
     The powr function defines all three indeterminate forms as invalid operations and 
     so returns NaN. 
3. Real operations with complex results, for example: 
     The square root of a negative number. 
     The logarithm of a negative number 
     The inverse sine or cosine of a number that is less than −1 or greater than +1. 

現在,這可以幫助你解決你自己的問題。

+0

我看到了,但我的情況沒有人。 – Alex

+0

請確保因爲您的算術運算對值't','* Xp','Xp'在'0'和'FLT_MAX'之間執行得很好。 –

0

您必須打印每個計算結果的子結果 - 或使用isnan()函數在常規地點檢查並追蹤其來源。這可能是一些「壞」的數學,或者你首先在垃圾中餵食(未初始化的變量可能是NaN)

3

有第五可能性埃裏克Postpischil的,否則優秀的答案沒有提及:

  • 乘法在的x87寄存器執行,以及浮點棧溢出發生了起因到可能不相關的早期操作在你的程序執行中。當處理器處於這種故障狀態時,在x87寄存器上執行的全部計算產生NaN結果。

這兩個最常見的原因是調用返回浮點結果的函數,該結果在範圍中沒有原型(有很多調用約定,這會導致調用者無法將結果從FP中彈出堆棧)和不正確的手寫(可能是內聯)程序集。

失敗只發生在一段時間之後的事實爲這種可能性提供了一些證據;如果有一個很少使用的代碼路徑泄漏了浮點堆棧的一個元素,那麼在失敗清單出現之前需要使用一些次數,這可能會讓它直到現在才能通知它。

要診斷或排除這種可能性,您需要查看浮點狀態寄存器(FPSR)的位6(SF)。根據您使用的編譯器,檢查FPSR的確切方法可能會有所不同。