Intel CPU 的三角函數精度問題，會對現實生活產生什麼影響？

01-12

英特爾處理器三角函數指令被發現會產生極其不精確的結果_Intel 英特爾_cnBeta.COM

e一般做sin/cos是先把輸入截斷到[-Pi/2, +Pi/2]，然後用泰勒級數之類的多項式級數逼近或者cordic演算法算出來。

intel里文檔里解釋產生問題的原因是intel在把大輸入數截斷到[-Pi/2, +Pi/2]過程，也就是sinf(x)=sinf(x-N * Pi/2)時，把pi的近似值作為常數放到了一個66位的ROM里。

對於大數來說，這個精度是不夠的(想想把pi的幾億倍，小數點後剩下的有效數字還有多少)，而sin/cos這種恰恰對小數很敏感。

一般libm里干這個事情的單精度版本e_rem_pio2f() libm/upstream-freebsd/lib/msun/src/e_rem_pio2f.c

裡面都用了3個double存pi,相當於(52+1)*3=159位有效數字。arm上截斷到[-Pi/2, +Pi/2]花的時間不比後面算多項式短。

像這種就有可能造成地球地圖之類的在離坐標原點遠的地方出現鋸齒，或者線歪掉。

因為是pi的值出了問題，所以應該是對大輸入引入了一個N*(Pi-Pi")的偏移量，輸出隨著輸入增大而越來越不準確（絕對誤差）。由於三角函數的周期性，當輸入在Pi的整數倍附近時，相對誤差尤其大（輸出數量級很小，看上去錯誤的ULP更多了）。

問題還是比較嚴重的。

——————————————————————

感謝評論中指出是fsin。單精度雙精度都會錯。誤用的話問題可能比較大，如果在錯誤的範圍調用，或者編譯器直接把三角函數編譯成快速的超越函數指令，都會產生錯誤結果。

按照intel解釋，fsin保證的是歸一化以後1～1.5ulp精度，也就是fsin比較適合作為快速數學函數計算小輸入或者在前述sinf/ sin實現中替代多項式級數計算部分。

有些Gpu里也有快速數學函數，一樣是大輸入不精確。

手冊寫清楚避免誤用就行了。

在cnbeta我自己發的評論，估計還在審核，直接轉過來：

按上面的意思測試了一下，sin(1e16)應該是0.7796880，但是算出的是0.7796799，這麼看的確精度很低。但是實際上sin值是0.7796799的數大約是1e16+0.000013，以long double的精度本來就是無法表示的。可能intel在設計的時候，認為精度不夠的時候最後一位之後的部分可以當成任意值，但是某人認為後面應該當成全是零，而且找到了一個需要後面全是零的應用，而且發現官方文檔的意思與實際表現不同，所以提出來了這個問題。簡單地說：intel的處理方法按照常識並不會造成什麼問題，但是文檔沒寫清楚，而且存在某種常識以外的應用。

這個問題的原因 Shu Zhang 已經解釋的很全了。

說點別的吧，其實這裡說的精度蠻不準確的。

額，首先我們要先了解Accuracy和Precision，一般來說我們所謂的精度是precision，也就是雙精度單精度擴展精度里的精度，這裡的精度是指accuracy，即準確度，關係如下圖。

對於Accuracy而言，一般測試結果表示為ULP，通俗的解釋就是如果1.0和2.0之間有10個數，即ULP為0.1，如果準確度為2ULP，那麼如果一個實際結果為1.25那貌似就因為偏差可能在1.05到1.45之間挑，也就是1.1、1.2、1.3、1.4之類都可能出現的樣子。

然後其實我們電腦進行的很多浮點計算都並不是完全準確的（+-*等是準確的，別的都有誤差存在），以下是intel手冊里關於X87 FPU準確度的一部分。

可以看到，其實並不能完全準確。

對於CPU的FPU如此，對顯卡更是這樣，NVIDIA的資料中也寫到https://developer.nvidia.com/sites/default/files/akamai/cuda/files/NVIDIA-CUDA-Floating-Point.pdf 大多數sin正確，部分1ULP，少數2ULPs之類。AMD貌似2-3ULPs吧，一般CPU比GPU要準確。這裡有比較http://iris.lib.neu.edu/cgi/viewcontent.cgi?article=1057context=elec_comp_theses。

在OpenCL 1.2中關於各個運算的最小準確度要求的更低。一般科學計算也主要是使用X87 FPU的擴展精度進行以平衡準確度和速度。

所以我們其實一直是生活在錯誤里的......只不過intel這個問題的錯誤較大。

另，我不是CS專業的，關於體系結構也不懂，求輕拍。

呃........

因為工作關係測過這個玩意，擔心會影響testcase的正確性，於是翻了翻ubuntu的源代碼。

[svn] Index of /trunk/libc/sysdeps/x86_64/fpu

http://x86.renejeschke.de/html/file_module_x86_id_115.html

結果發現：

assembly - Calling fsincos instruction in LLVM slower than calling libc sin/cos functions?

fsincos用了硬體，sinf和cosf還是彙編。

/* ------------------------誤--------------------------------- 但ubuntu裡面用的eglibc沒有做任何的優化，還是用了fsincos。

如果用ubuntu做科學計算的同學，要注意了，還是乖乖的自己寫庫吧。要不看看GNU那邊的libc是不是優化了。 */

不過幸好，我的數據源在很小的範圍內，還trigger不到這個bug。

工程設計……

摺疊