2012-01-25 38 views
2

我想知道如何在沒有分割硬件和沒有浮點硬件的二進制中實現IEEE-754 32位單精度浮點除法?如何在沒有分割硬件和沒有浮點硬件的二進制中實現浮點除法

我有移動硬件,加,減,乘。

我已經實現了使用16位字的浮點乘法,加法和減法。

我正在一個專有的多核處理器上實現這些指令,並在彙編中編寫我的代碼。之前,我正在使用matlab來驗證我的算法。

我知道我需要減去指數,但是如何執行unsigned division對mantissas?

+1

你沒有包含了代表你的意思和你使用的語言的任何示例代碼。他們需要回答您的查詢。 – Lion

+1

如果你傾倒了一千條堆棧跟蹤或整個程序,那隻會有太多的毛病。 – BoltClock

+0

@BoltClock有時候更少。 – Veridian

回答

3

取決於你想如何複雜。保持它相當簡單,你可以嘗試用倒數近似來劃分。

而不是計算:(n/d)你會發現:n *(1/d)。

要做到這一點,你需要使用一些方法來計算互惠,例如,Newton-Raphson它使用牛頓的方法來連續計算除數的倒數的更準確的估計,直到它對你的目的「足夠」精確做最後的乘法步驟。

編輯

剛剛看到您的更新。這可能會,也可能不會對你有用!

2

如果您的硬件具有足夠快的整數倍數(比如說4-5個週期),則使用迭代方法計算recip = 1/divisor可能是最快的方法。然後,您將計算商數爲dividend * recip。如果硬件提供帶有雙倍寬度結果的整數乘法(即完整乘積)或者提供兩個32位整數乘積的高32位的「mulhi」指令,則對於必要的定點計算非常有用。

您將需要重新調整定點計算以保留接近32位的中間結果,以得出精確到​​舍入後最終結果所需的24位的結果。

最快的方法是有可能產生一個9位的開始近似「約」〜1 /除數後跟一個立方收斂迭代的倒數:

e = 1 - divisor * approx 
e = e * e + e 
recip = e * approx + approx 

它是最簡單的預先計算的起始近似和sstore它由一個256字節的數組構成,由除數的位23:16(即尾數的8個最重要的小數位)索引。由於所有近似值都是0x100到0x1FF(對應於0.5到0.998046875)範圍內的數字,因此不需要存儲每個值的最高有效位,因爲它將爲'1'。只需將0x100添加到檢索的表格元素中即可獲得初始近似值的最終值。

如果你不能承受256字節的表存儲,另一種產生精確到9位的開始近似值的方法將是一個3次多項式,近似爲1 /(1 + f),其中f是小數部分除數尾數,m。由於使用IEEE-754,已知m在[1.0,2.0)中,f在[0.0,1.0)中。

可以通過用除數反推初步商以建立餘數並選擇最終商以使餘數最小化來實現正確舍入到最接近或偶數(如果需要)。

下面的代碼演示了上面討論的近似原理,使用簡單的倒數情況,根據IEEE-754最近或偶數模式進行適當舍入,並適當處理特殊情況(零,非正規,無窮大,NaNs)。它假定一個32位IEEE-754單精度浮點數可以按位從一個32位無符號整型數據傳輸。所有操作都是在32位整數上執行的。

unsigned int frcp_rn_core (unsigned int z) 
{ 
    unsigned int x, y; 
    int sign; 
    int expo; 

    sign = z & 0x80000000; 
    expo = (z >> 23) & 0xff; 
    x = expo - 1; 
    if (x > 0xfd) { 
    /* handle special cases */ 
    x = z << 1; 
    /* zero or small denormal returns infinity of like sign */ 
    if (x <= 0x00400000) { 
     return sign | 0x7f800000; 
    } 
    /* infinity returns zero of like sign */ 
    else if (x == 0xff000000) { 
     return sign; 
    } 
    /* convert SNaNs to QNaNs */ 
    else if (x > 0xff000000) { 
     return z | 0x00400000; 
    } 
    /* large denormal, normalize it */  
    else { 
     y = x < 0x00800000; 
     z = x << y; 
     expo = expo - y; 
    } 
    } 
    y = z & 0x007fffff; 
#if USE_TABLE 
    z = 256 + rcp_tab[y >> 15]; /* approx */ 
#else 
    x = y << 3; 
    z =     0xe39ad7a0; 
    z = mul_32_hi (x, z) + 0x0154bde4; 
    z = mul_32_hi (x, z) + 0xfff87521; 
    z = mul_32_hi (x, z) + 0x00001ffa; 
    z = z >> 4; 
#endif /* USE_TABLE */ 
    y = y | 0x00800000; 
    /* cubically convergent approximation to reciprocal */ 
    x = (unsigned int)-(int)(y * z); /* x = 1 - arg * approx */ 
    x = mul_32_hi (x, x) + x;  /* x = x * x + x */ 
    z = z << 15; 
    z = mul_32_hi (x, z) + z;  /* approx = x * approx + approx */ 
    /* compute result exponent */ 
    expo = 252 - expo; 
    if (expo >= 0) { 
    /* result is a normal */ 
    x = y * z + y; 
    z = (expo << 23) + z; 
    z = z | sign; 
    /* round result */ 
    if ((int)x <= (int)(y >> 1)) { 
     z++; 
    } 
    return z; 
    } 
    /* result is a denormal */ 
    expo = -expo; 
    z = z >> expo; 
    x = y * z + y; 
    z = z | sign; 
    /* round result */ 
    if ((int)x <= (int)(y >> 1)) { 
    z++; 
    } 
    return z; 
} 

函數mul_32_high()是用於返回兩個32位位整數有符號乘法的高32位的特定機器操作的佔位符。替代的特定機器的版本的半便攜式實現

/* 32-bit int, 64-bit long long int */ 
int mul_32_hi (int a, int b) 
{ 
    return (int)(unsigned int)(((unsigned long long)(((long long)a)*b)) >> 32); 
} 

由基於表中,所使用的256項倒數表可以構造如下:

static unsigned char rcp_tab[256]; 
for (int i = 0; i < 256; i++) { 
    rcp_tab[i] = (unsigned char)(((1./(1.+((i+.5)/256.)))*512.+.5)-256.); 
} 
相關問題