這個求指數函數exp()的快速近似方法的原理是什麼?

01-13

下面的代碼可以快速求出指數函數exp()的近似解, 裡面用到了

reinterpret_cast以及兩個神奇數, 它們的原理是什麼?
inline double fast_exp(double y) { double d; *(reinterpret_cast&(d) + 0) = 0; *(reinterpret_cast&(d) + 1) = static_cast&(1512775 * y + 1072632447); return d; }
經測試輸入值在-10.0~10.0時誤差在-4%~2%之間波動:

下面鏈接中還有更一般的求a的b次冪的函數pow(a,b)的近似解:
Optimized Approximative pow() in C / C++

計算近似值最重要的一步是計算出雙精度浮點數的階碼，階碼是這麼個數：

$1023 + lfloor log_2 d floor$

也就是說要計算

$1023 + lfloor log_2 exp(y) floor$

先取指數後取對數，這個我們很熟悉，這就變成乘法了：

$1023 + lfloor y cdot frac{1}{ln 2} floor$

那麼尾數的那部分呢？它實際上是

$2^{{ y cdot frac{1}{ln 2} }} - 1$

減去1，是因為浮點數已經親切地替我們把最高位的1省略掉了。

也就是2的小數部分次冪。

現在我們把整個雙精度浮點數的前32位看成是一個定點數，小數點在階碼的後面，我們會想辦法把這個整數部分和小數部分湊成同一個乘法，也就是：

$R = 1023 + lfloor y cdot frac{1}{ln 2} floor + 2^{{ y cdot frac{1}{ln 2} }} - 1 = 1023 + y cdot frac{1}{ln 2} + 2^{{ y cdot frac{1}{ln 2} }} - 1 - { y cdot frac{1}{ln 2} }$

剩下最後幾項，由於是小數部分組成的，範圍有限，我們用一個常數來代替它。考慮到 ${y cdot frac{1}{ln 2}}$ 的範圍是[0,1)，我們可以用積分平均值來取代這個數：

$int_0^1 (2^x-1-x) dx = left.left(frac{2^x}{ln 2} - x - frac{x^2}{2} ight) ight|^1_0 = -0.0573$

也就是

$R = 1022.9427 + y cdot 1.4427$

用整數來代替浮點數，由於階碼後面有52位的位數，扣掉最後32位不要，還有20位，所以要乘以 $2^{20}$ ，四捨五入也就是：

$R = 1512775 * y + 1072633159$

大致與題目中的公式接近，題目中最後尾數部分的一點微小差別應該是通過實驗確定的最小誤差的值吧。

========================================================

這個方法可以擴展到很多其他運算上，主要就是階碼 + 尾數構成的定點數，與原浮點數的2-base的對數接近相等的原理。反過來，我們將一個浮點數強制轉換成整數，然後減去一個常數，就可以得到一個定點數版本的對數值；那麼許多運算都可以用對數來代替，比如說開平方根，先用這個方法轉換成定點數版本的對數值，然後右移一位，再補回原來的常數然後強制轉換回到浮點數，這就是著名的magic number的原理；比如說任意兩個數的乘冪，可以寫成前一個數的定點數版本的對數值與浮點數相乘，然後轉換回浮點數。不僅對於雙精度數，成立對於單精度數也是一樣的。

原理非常簡單，思路非常巧妙，全靠ieee浮點數構成形式：

我們要求exp(x)，也就是 $e^x$ 並得到一個double的結果。

觀察double浮點數的構成形式：

$Y = M * 2^E$

其中E為指數，M為尾數。

只觀察指數部分，會發現實際上它們均為 $X^Y$ 形式的。

那麼就有意思了，我們要求的是 $X^Y$ 形式的一個數，而由於都是浮點表達，我們的輸入x也是

一個 $X^Y$ 的數，那麼，這之間有沒有什麼可以玩的東西呢？

慢慢來，我們先來看看 $2^Y$ 如何簡便計算:

假設我們有個浮點數 $Y$ ，要計算 $2^Y$

觀察double浮點數的構成形式：

$Y = M * 2^E$

我擦，原來 $2^Y$ 已經藏在結構裡面了啊！就是 $E$ 啊！

所以說，只要我們把 $Y=E$ ,再把 $Y$ 這個數放在正確的bit位置上，我們就構建好這個數了。

那麼，如何放置呢？

我們觀察ieee的double浮點數的構成形式：

第零位bit是符號位

往後11位是指數位

往後剩下的是尾數位

也就是說，我們只需要把 $Y$ 轉成int，然後放入這個浮點數的第63-52位上就行了！

那麼，怎麼放法？

首先，我們將這個double轉為int，得到一個32位的int32，再將這個int32向左移動20位，

這樣，這個int32的前12位就挪到了這個int32的最高12位上，這樣這個int32就和ieee的double浮點數的前63-32位對應上了。

但是別忙，題主的方法中，並不是直接使用移位計算的，而是直接使用乘法，乘上了 $2^{20} = 1048576$ ，我們知道乘2^20和移位能達到同樣的效果，這也就是代碼里的這一步：

1512775 * y

你會說啦，那不對啊，這明明乘的是1512775啊，對，那是因為求的是e的pow而不是2的pow。

e的咱們待會兒再說。

再然後還沒完，後面為什麼還帶一個

+ 1072632447

這是因為，double的指數部分的int並非常見補碼形式的int，而是一個帶有偏移量的int，偏移量為2^12次方的一半也就是1023，大於1023的數值才被認為是正數。因此我們要加上這個偏移量，即加上1023*2^20=1072693248。

那你就又問了，不對啊人家明明加的是1072632447啊，沒錯，多出來的60801是一個補償值，

具體是如何補償的，詳見「A Fast, Compact Approximation of the Exponential Function」的第四節。

綜上，題主的代碼就好理解了：

double d;


 // 先將尾數的後32位抹零。

*(reinterpret_cast&(d) + 0) = 0;

//再計算指數位，移位，加上偏移量和補償值 *(reinterpret_cast&(d) + 1) = static_cast&(1512775 * y + 1072632447);

這麼做之後，大家有沒有發現有什麼問題？有的！那就是尾數並不都是0！

雖然後32位被抹零了，但前面指數加尾數的部分：

並不都為0，因為這一步：

//再計算指數位，移位，加上偏移量和補償值 *(reinterpret_cast&(d) + 1) = static_cast&(1512775 * y + 1072632447);

這一步中的(int)(1512775 * y + 1072632447)這個int32隻保證前面12位是指數沒問題，可後面還跟著20位呢！什麼鬼？

首先，跟著的這20位是加到尾數里的，本身影響就不大，另外人家論文里特意說了，跟著的這20位不但不會降低精度反而對精度有幫助，具體的我實在沒仔細看，有興趣那你只好去看論文了。

至此，一個簡單的powerOfTwo的快速計算就完成了。那麼你就問了，說了半天，exp可怎麼辦呢？

能算2就能算e！觀察此式：

$e^{b}$

有如下等式：

$e^{ln(a^b)} = a^b$

有：

$e^{ln(a)*b} = a^b$

令a=2：

$e^{ln(2)*b} = 2^b$

令x=ln(2)*b,則有：

$e^x = 2^{x/ln(2)}$

也就是說，我們要計算exp(x)，只需要計算2^(x/ln(2))就行了。而計算2^x的方法咱們前文已經講了。

那麼，重新觀察前文中的這一行：

1512775 * y

剛才說了，這裡應該乘1048576，為什麼乘了1512775？因為它算的是e，不是2，所以乘的是1048576 / ln(2) !

至此，這個簡便方法的大致步驟就講完了。

這個演算法的誤差主要來源於主要發生在這裡：

static_cast&(1512775 * y + 1072632447);

這個double轉int時，小數部分應該被計入最終數字中的尾數，但卻被忽略了，所以會產生一個來回晃動的誤差。就像題主那張圖那樣。

IEEE754規定的，浮點數表達方式是S+E+M

S 符號位，尾數的符號位；

E 階，即指數；

M 尾數，即有效小數位數；

其中，雙精度浮點數double的三個參數分別是

符號位 1位，bit63

階 11位，bit62~52

尾數 53位，bit51~0

三個係數生成的浮點數值為

$M imes 2^E$

然後，reinterpret_cast&(d)這一句的意思就是，把這串浮點數直接按照整形規則來解析

也就是相當於

$2^{63} S+2^{52}E+M$

計算機領域中，有相當多的快速演算法，就是用近似計算，把圖一作為x，則常規演算法是f(x)

然後通過強制類型轉換，把浮點數轉為整數，相當於把x轉換為圖二表示的x"，然後將演算法用數學變換，加上一些近似處理，轉換成整數計算的新演算法f(x")。

近似的方法並不複雜，比如泰勒展開之類。作近似轉換後，把係數擬合到一定精度的常量，就可以生成如題目中那樣的演算法，那些常量就是magic number

具體的思路，可以看這個答案0x5f3759df這個快速開方中的常數的數學依據是什麼？ - Belleve 的回答

這是這方面演算法中，最經典的，快速計算平方根倒數演算法的解釋，比我簡單寫的這個答案要更加詳實，並且更加容易看懂

這個我記得是一個purdue university的數學教授寫的paper。這個演算法還是很有用的，很多需要用到exp函數（譬如各種regression），而精度要求有限的時候，這個近似演算法基本上可以提高几十倍的速度。如果exp本身是程序的瓶頸之一，那對性能的提升就很客觀了。以前的工作就遇到過這樣的情形。