推 jchin : 在x=0或1, 分子分母取微分? 07/19 06:02
→ yhliu : 當遇到分母為 0 或絕對值低於某值，檢查分子是否亦 07/19 07:57
→ yhliu : 接近 0, 如是，則 x 代以另一近似值計算之；否則ERR 07/19 07:59
→ yhliu : 分子分母微分法的好處是此例只要 x 靠近 0 或 1 都 07/19 08:14
→ yhliu : 能得到正確結果。而上面我提出的方法好處是適用一般 07/19 08:15
→ yhliu : 0/0 不定式，而不需改函數式，缺點一當然在此例就是 07/19 08:17
→ yhliu : x 靠近 0 或 1 會遇到除以 0 或溢值問題，另外就是 07/19 08:18
→ yhliu : 結果也許不精確，例如第一個函數在 x<=0.1e9 時開始 07/19 08:21
→ yhliu : 不正確。 07/19 08:21
謝謝j大跟y大的idea, 感覺不管微分還是近似值計算, 都還是要決定一個threshold?
即大於某個threshold採取原式, 小於某個threshold採取近似式?
原本是不想要有這個"tunable parameter", 但是目前看起來只能這樣了QQ
好像跟偽逆矩陣一樣, 多小的eigenvalue視作0, 否則視作倒數
推 j0958322080 : 試試 bug number algo. 07/19 13:41
→ j0958322080 : bug-->big 07/19 13:41
謝謝j大關鍵字
→ simonjen : 你可以自己寫一個運算用字串去存 07/20 08:45
推 LPH66 : 樓上那就是上面提的 big number (大數) 演算法 07/21 17:11
我查big number algorithm沒看到特指哪一種算法, 以及s大說的"用字串去存"指的是什麼
有Reference可以提供嗎? 謝謝!
→ PPguest : 我估狗搜尋一下,看起來是指用其他方法如字串,陣列來 07/21 22:17
→ PPguest : 存例如100位有效數字的數 07/21 22:19
推 LPH66 : 大數的基本概念就是這樣沒錯, 用大陣列去模擬小學時 07/22 17:30
→ LPH66 : 做基本加減乘除的演算法 07/22 17:30
→ LPH66 : 例如 GMP (GNU 多重精度運算庫) 就是一個用 C 語言 07/22 17:32
→ LPH66 : 寫的這類型函式庫 07/22 17:32
喔喔 意思就是多開buffer去存更多資料就是了?
推 j0958322080 : 不一定要字串，你也可以開個 int64_t 再注意進位 07/22 20:01
推 wohtp : f(x)在 x = 0 附近改寫一下 07/22 22:13
→ wohtp : f(x) = -1/(x-1) + (exp(x ln 2) -1)/x * 1/(x-1) 07/22 22:16
→ wohtp : 主要問題是 (exp(x ln 2) -1)/x 這項 07/22 22:17
→ wohtp : 我不知道你用什麼，但是大部分函式庫應該都有expm1 07/22 22:18
→ wohtp : 改用這個就不會有誤差問題 07/22 22:20
嗨w大, 第一次看到expm1, 然後查一下還有logp1來增加精度的case, 我再參考, 謝謝!
→ wohtp : f(x)拆成兩截，x<0.5的時候用這個，另一頭改用y = 1 07/22 22:22
→ wohtp : -x 當變數，然後expm1比照辦理 07/22 22:22
了解~
推 wohtp : 傷浮點精度的是大數相減剩小數，不是0/0 07/22 22:26
f(x)/g(x)時, 假設數學上是lim_{x→0}f(x) = lim_{x→0}g(x) = 0
且lim_{x→0} f(x)/g(x) = 1
我就是怕x→0時, f跟g這兩個函數在library上的精度不同, 或是函數本身遞增速率不同
導致某個很小的x代入g(x)時, 如果電腦已經因為精度問題return true 0
但是f(x)還沒, 那f(x)/g(x)就是+inf
反之如果是f(x)先return true 0, 那f(x)/g(x)就是+0
我就是怕以上這種情況才會問這篇, 還是說這個並不會發生?
→ PPguest : 原來有expm1(筆記) 07/22 22:36
→ PPguest : 想請問一下,若想知道實際誤差的狀況,準確的計算值要 07/22 22:37
→ PPguest : 如何得到？ 07/22 22:37
推 wohtp : 這個問題不就是「怎麼讓數值運算更精確」 07/22 22:45
→ wohtp : 開100位浮點數去算啊，不然還能怎樣？ 07/22 22:46
推 wohtp : 喔對了，回到原po的問題，x = 0 or 1的時候要直接讓 07/22 22:52
→ wohtp : 函數給極限值，所以其實需要四個case 07/22 22:52
我真實要應用的case是(1,0), (s,log_2(s)), (p,log_2(p)), (2,1)
然後1, s, p這三個可能都相等, 可能都相異, 可能兩個相等
然後我去觀察圖形時哪個case, 都還是三次多項式
因此要拆的段數, 要算的極限值好多...
推 wohtp : 所以你的函數最複雜就log？那你怕什麼？ 07/23 01:37
→ wohtp : 這些常見的函數實作精度大概都逼近machine precisio 07/23 01:37
→ wohtp : n。 07/23 01:37
→ wohtp : 你給的例子裡面，精度直落的唯一原因是大數相減，完 07/23 01:37
→ wohtp : 全是程式員的鍋。 07/23 01:37
推 wohtp : 不然就像上面建議的，開bigfloat硬幹吧。你的時間應 07/23 01:47
→ wohtp : 該比cpu時間貴……吧？ 07/23 01:47
我應用的例子是不會有大數相減, 單純就是一堆0+/0+
如果照w大你所說的這些函數的精度幾乎沒問題, 那我就放心了
另外問個浮點數"可能出問題"的問題, 今天假設我避開大數相減...等等會大幅喪失精度
的事情, 設計出一個C函數"float f(float x)"
然後我確定其數學理論式中, x!=0的話f(x)就沒問題
但是我還是怕浮點運算中即便x!=0, f(x)也可能會出問題(如這篇討論的問題)
我除了暴力跑IEEE754的所有denormal number然後print出來檢查, 是不是沒其他方式了
換言之, 有沒有怎樣的數學式的coding style可以完全避開遇到這些精度問題?
問工程的朋友目前沒有general的答案, 都要看函數長相
推 LPH66 : 大減大消光這有個名字叫 catastrophic cancellation 07/23 02:42
→ LPH66 : 一般都是只能建議改動式子不要有讓相近數相減的地方 07/23 02:42
→ LPH66 : 但變形方式如你所問到的都得看函數長相才知道 07/23 02:43
→ LPH66 : 你的 0+/0+ 應該是踩到浮點數指數下限所以下溢了 07/23 02:45
→ LPH66 : (denormal number 出現的原因就是緩和接近下限時的 07/23 02:45
→ LPH66 : 下溢狀況, 但它只對部份運算有些許緩和作用而已) 07/23 02:46
還有catastrophic cancellation這專有名詞, 謝謝提供
然後02:42說的"不要有讓相近數相減的地方", 是不論大數或是小數相減都應該避免?
大數相減誤差大(浮點在大數精度很低), 小數相減又有下溢問題...
另外重複的問題是, 如果我自認為某個函數實作我都盡量避免掉這些問題
但是我想要把信心度提高到100%, 就是掃整個浮點數的範圍?
推 wohtp : 我說大數小數是相對的啊。要說得更準確些是「相近數 07/23 09:44
→ wohtp : 相減剩下比原數小得多的差」。 07/23 09:44
→ wohtp : 拿你的f(x)為例，2^x的精度不會有問題，但是2^x - 1 07/23 09:44
→ wohtp : 在零附近就爛掉了。出現這個組合是你自己的鍋，不是 07/23 09:44
→ wohtp : library精度不行。 07/23 09:44
這個我在吸收一下, 感覺我沒搞懂 "單一數的精度問題" VS "相減結果的精度問題"
numpy文件舉例exp(10^-10)-1跟expm1(10^-10)的差異確實差很多
但是我自己舉例: float x = (1.[x_1~x_23])_2 * 2^0, where x_1~x_22=0, x_23=1
float y = (1.[y_1~y_23])_2 * 2^0, where y_1~y_23=0
其中()_2代表用二進位表示, []裡面裝著0 or 1的數
可以知道y就是數學上的1.0, x是float中比1.0大的下一個數, 也就是最接近且大於y的
而觀察z=x-y= (1.[z_1~z_23])_2 * 2^(-23), where z_1~z_23=0
可以知道相減後的結果仍是float, 並且毫無誤差
目前有點混亂, 不知道是我舉例錯誤還是剛好舉到特例, 還是根本跟expm1的案例不同類
而且如果自己刻比double更大的精度, exp(10^-10)-1的精度也會上升直逼expm1(10^-10)
所以感覺綠色的VS不太能完全獨立? 彼此有關係
推 wohtp : 然後，要保證實作數值方法的品質，你當然對原函數的 07/23 09:58
→ wohtp : 性質要有足夠理解。解析方法手算再每個點都挑掉是唯 07/23 09:58
→ wohtp : 一做法，就算像你說掃整個浮點數範圍，意義其實也只 07/23 09:58
→ wohtp : 是輔助判別有麻煩的值，最後你還是要自己把函數改寫 07/23 09:58
→ wohtp : 成沒有精度問題的形式。 07/23 09:58
那可以歸納成 "盡量所有的浮點coding都避開相近數相減" 就好了?
比如我在實作浮點二次多項式P(x)=x^2-x, 我都直接打x*x-x
今天如果我沒有要拿P(0+)去除以一個會接近0的分母, 我維持x*x-x就好
但是如果有要除法, 我就要改寫x*x-x了?
推 PPguest : 原po擔心的0+,就像LPH66講的,應是怕函數算一算結果 07/23 10:50
→ PPguest : 超出浮點數能表示的範圍 07/23 10:50
→ PPguest : 雙精度大約比5e-324小一半就會變0 07/23 10:54
→ PPguest : 感覺就是注意實際運算中不會發生這種情況,如果會發 07/23 11:04
→ PPguest : 生,就代表雙精度不夠用 07/23 11:04
→ PPguest : 例如f(x)=x^2, x有需要用到例如1e-170那麼小嗎？ 07/23 11:29
→ PPguest : 若否,那沒事;若有需要,那雙精度就不夠用 07/23 11:32
嗨P大, 如果單純是f(x)=x^2並且沒有f(0+)/0+
那麼f(0+)的精度我確實也覺得不重要(在我的應用中), f(0+)=1e-170 or 0.0都沒差
只是目前我對於w大的例子" expm1 " & "2^x-1"有感覺
w大又提"2^x-1"的問題並不是2^x的精度問題, 但是精度提升不是會更精準嗎
而自己舉的float x, y, z的相減例子又覺得沒問題
目前有點混亂, 是不是要先固定某個精度後(比如float)才能談相減的誤差
有種...某個例子會這樣看起來是某個問題, 但是另外一個例子看起來又不是那個問題XD
→ PPguest : 主要是 07/22 22:26 後面你說的那種情況 07/23 16:21
→ PPguest : f(x)=x^2, g(x)=x, x=1e-170 07/23 16:22
→ PPguest : f(x)/g(x)=0, g(x)/f(x)=Inf 07/23 16:22
了解~
推 LPH66 : Catastrophic cancellation 的問題只要你有一個中間 07/23 18:00
→ LPH66 : 使用的精度, 那就有可能發生 07/23 18:00
推 LPH66 : 例如 x=1e-10, 你直接求 2^x 再減 1 的話 07/23 18:05
→ LPH66 : 2^x 得 1.00000000006931 (拿十五位十進位有效舉例) 07/23 18:06
→ LPH66 : 減去 1 之後得到 6.931e-11 這個數只有四位有效 07/23 18:06
→ LPH66 : 而如果用其他方式去直接求得 2^x-1 的話 07/23 18:07
→ LPH66 : 或許你能得到 6.93147180583968e-11 這樣十五位有效 07/23 18:08
→ LPH66 : 這個"其他方式"就是運用數學關係嘗試找出沒有這種減 07/23 18:09
→ LPH66 : 的式子 (以 2^x-1 來說可能是泰勒展開後消掉 1) 07/23 18:09
→ LPH66 : 上面這個例子你就算中間用了三十位有效 07/23 18:21
→ LPH66 : 減 1 對消後的有效位數就是不足三十位 07/23 18:21
→ LPH66 : 這就是它被叫做 catastrophic (災難性) 對消的原因 07/23 18:22
L大這個舉例我了解了, 我如此解釋:
(1) 2^(1e-10)=1.00000000006931, 這並不是精度不足
這個值就是在固定精度(假設系統是十五位十進位的系統)下所求出來的最準確值了
理論值就是1.00000000006931abcdefg...
這也是為什麼w大說不是library 2^x的精度問題
(2) 2^(1e-10) - 1的理論值是0.00000000006931abcdefg
而這個理論值在系統是十五位十進位是會寫成6.931abcdefg...*10^(-11)
也就是說, 這個6.931abcdefg...*10^(-11)這個值才是2^(1e-10) - 1
在這個系統的最準確值
但是今天如果照(1)的算法, 我們卻得出6.931*10^(-11)這個不準確的值
這樣理解應該是正確的吧?
---------------------------------------
另外請問各位一個實作問題, 我一直以來都在嵌入式系統寫定點(int32)程式碼
最近第一次有浮點可以用的IC, 寫code的速度飆快, 因為之前算出的數學式就直接刻上去
不用像定點還要控制數值範圍卡overflow之類的
結果最近這幾個需要在意精度的浮點公式讓我有種浮點比定點更難計算與控制誤差的感覺
所以想請板上浮點高手給些建議怎樣的浮點coding style能少踩一點坑
比如之前說的少用相近數相減
→ PPguest : 原po 07/23 09:44 之後紅色那段的例子 07/23 22:45
→ PPguest : 基本上都是10進位可以準確換算成2進位的情況 07/23 22:47
→ PPguest : 例如0.5, 0.25, 0.125。像0.1用二進位的角度來看,在 07/23 22:49
→ PPguest : 浮點數的系統只能用近似值來表達 07/23 22:49
這個我理解, 所以我09:44那段特地取剛好可以表達的浮點數x,y來減少討論的變因
不過這樣看起來會有catastrophic cancellation就是不會發生在這種剛好可以表達的
→ PPguest : 很多數用二進位的角度來看,都是無窮小數,然後搭配 07/23 22:53
→ PPguest : 原po對L大舉例的理解,應該就沒問題了 07/23 22:54
了解~
→ PPguest : 用比double更大的精度, exp(10^-10)-1 大概會看到更 07/23 22:58
→ PPguest : 多小數,但大概有一定比例的小數其實是不準的 07/23 22:59
如果發生這種事可以先歸類在函數實作不精準吼?
前提應該是基於在某個精度系統下, 函數的實作準度是準到尾數的bits數
→ PPguest : 不過小數準的位數應該會變多 07/23 23:01
→ PPguest : "是不是要先固定某個精度......才能談相減的誤差" Y 07/23 23:05
了解~
推 PPguest : 多項式用"nested"的方式求值除了降低計算量,好像也 07/23 23:10
→ PPguest : 算是減少誤差 07/23 23:11
推 PPguest : 07/23 10:50 那段，原po實際應用的情況大概不用擔心 07/23 23:17
→ PPguest : 會發生.s和p在1,2之間,如果不同,再怎麼近大概是 07/23 23:20
→ PPguest : 2e-16等級的差距,要超出浮點數可表示的範圍應該不太 07/23 23:22
→ PPguest : 可能 07/23 23:22
謝謝idea, 放心不少~
→ PPguest : 07/23 22:59"有一定比例的小數...不準" 我用詞不好 07/24 11:15
→ PPguest : 單純是相近的數相減造成不精準 07/24 11:16
→ PPguest : 如果原本有15位有效數字,相減後剩5位有效,後面10位 07/24 11:19
→ PPguest : 不準.用比double更大的精度,計算前例如有30位有效數 07/24 11:20
→ PPguest : 字,相減後還剩20位有效,後面10位不準 07/24 11:21
→ PPguest : 用expm1,可能輸入是15位有效數字,輸出是14位有效, 07/24 11:23
→ PPguest : 1位不準;用比double更大的精度,輸入30位有效數字, 07/24 11:24
→ PPguest : 輸出29位有效,1位不準.大概是這種感覺 07/24 11:24
這個舉例理解~謝謝P大~
推 LPH66 : 初次接觸浮點數運算的話冼鏡光老師這篇可以看看 07/24 14:33
→ LPH66 : blog.dcview.com/article.php?a=Az0HYgNrBDU%3D 07/24 14:33
→ LPH66 : 包含了我們這裡提過的跟一些其他也很重要的運算概念 07/24 14:34
謝謝L大的reference~
※ 編輯: znmkhxrw (59.102.225.191 臺灣), 07/25/2023 01:53:08