作者LPH66 (-858993460)
看板C_and_CPP
標題Re: [問題] long long int 的餘數運算
時間Wed Feb 15 23:16:03 2012
※ 引述《tropical72 (藍影)》之銘言：
: ※ 引述《DJWS (...)》之銘言：
: : 感謝各位板友幫忙
: : 自己再嘗試回答一次
: : a*b 與 y*m 其實都會溢位
: : 但是我們知道正確的r一定小於m，一定放得進64bit
: 重點是在算 y 的時候就有捨位誤差了，到後面會與正確結果愈差愈遠，
: 這結果是你所能接受的？
: < 以下恕刪 >。
這樣一提我才發現果然這樣子還是會出事...
剛剛試了一下
傳 a = b = 4611686018427387902 (2^62 - 2)
m = 4611686018427387903 (2^62 - 1)
進去
理論上因為 a = b = m-1 所以 a*b mod m = 1
但實際跑程式的結果卻是 5
裡面的計算是這樣子的:
如上文 由於 double 的精確度關係 y 這商只有高位 53 bit 是正確的
於是當再乘上一個很大的 m 時 y*m 有效的那 53 bit 就完全不在低 64 bit 的範圍了
如此一來 a*b - y*m 的值會混入了不準確但也溢出 64-bit 的位數 就不再正確了
以此例來說 y 應該要是 m-2 (2^62-3)
但是因為只精準到高 53 bit 於是"算"出來的 y 值是 2^62
因此 y*m 的"值" 2^124 - 2^62 只有那 2^124 可信 -2^62 不在有效數字範圍裡
但偏偏因為溢位的關係丟掉的是 2^124
所以 a*b - y*m 的答案混入了 2^70 ~ 2^65 這些不準確的 bit 因此就不正確了
由此看來 m^2 < 2^(64+53) => m < 2^58.5 ≒ 4.076*10^17 應該是保證安全的大小
這樣即使 y 有誤差 但 a*b 和 y*m 溢出去的部份保證是一樣的
因此相減可以抵消 a*b - y*m 的值才保證是正確的
捨去誤差的部份只要全部都有被 a*b - y*m 掌握住那就沒什麼問題
--
其實我上面故意略過一個問題沒說
那就是
y 還是可能溢位!!
問題發生在當 y 接近 2^63 時 因為精確度只有 53 bit 的關係 它會被近似成 2^63
問題是有號 64-bit 整數是放不下 2^63 的
因此由浮點數轉回整數時它會被轉成 -2^63 這正是溢位
那麼 y 溢位的後果也就可想而知了...
--
'Oh, Harry, don't you
see?' Hermione breathed. 'If she could have done
one thing to make
absolutely sure that every single person in this school
will read your interview, it was
banning it!'
---'Harry Potter and the order of the phoenix', P513
--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 140.112.230.62
推 ericinttu:總之就是會爆炸 02/15 23:23
推 loveme00835:樓上 XD 02/16 01:05
推 tropical72:推導說明漂亮!! 推一個 02/16 14:22
推 DJWS:感謝!! 正在努力理解當中... orz 02/16 20:35
→ DJWS:另外回答tropical172上一篇的問題 這個程式碼是快速數論轉換 02/16 20:36
→ DJWS:的其中一段 他用的模數固定是50000000001507329LL 02/16 20:37
推 DJWS:to LPH66: 我想應該是 m^2 < (52+63)bit 有一個bit是負號 02/16 21:55
→ DJWS:to LPH66: 我想應該是 m^2 < (53+63)bit 有一個bit是負號 02/16 21:56
→ LPH66:嗯, 我在打那條時也有在想應該要扣掉負號才對.... 02/16 22:42
→ LPH66:5*10^16嗎...正好比這個 m 的極限值少一位 02/16 22:44
→ LPH66:如果固定是它的話那就沒問題 02/16 22:44
推 xatier:一樓精譬XD 02/16 23:29
推 DJWS:感謝LPH66 :D 02/17 20:00