不好意思，我是覺得這段 code 還是很有問題
另不知道這段 code 是不是就是 LCM, 只是把 constant 設 0 而已 ?
我說一下我的結論：這段 code 只能是近似算法，而不能算是準確算法。
下面是我提的論點。
※ 引述《DJWS (...)》之銘言：
: ※ 引述《DJWS (...)》之銘言：
: : 在網路上看到一段code
: : // r = a * b % m
目前我所知較穩定的作法如 LPH66 所言，
r = ( (a%m) * (b%m) ) % m;
重點是還有一個前提才不會出包： m <= sqrt ( long long int 最大值 )
: : long long int mulmod(long long int a, long long int b, long long int m) {
: : long long int y = (long long int)((double)a*(double)b/m+0.5);
: : long long int r = (a*b-y*m) % m;
: : if (r < 0) r += m;
: : return r;
: : }
: : 請問為什麼這樣寫就不會溢位？
----
要挑的話有二個地方是有點問題。
它可能不會產生溢位，但這和上面連續用到三個 mod 結果會不一樣。
(1) 想一下 a*b > 2^52 的時候會發生什麼事，
在轉 double 相乘的時候會造成捨位誤差。
(2) long long int y = (long long int)((double)a*(double)b/m+0.5);
這段有時候是沒意義的，如 a*b/m 結果大於 2^52 時，這0.5 加不進去。
又如前一篇所說，double 精度只有 52 + 1 位 [註一]，
意思是如果 a*b/m+0.5 其值大於 2^52 時，y 開始就會有捨位誤差
----
結論是，它只是在「趨近」算 a*b % m 之結果而已。
: 感謝各位板友幫忙
: 自己再嘗試回答一次
: a*b 與 y*m 其實都會溢位
: 但是我們知道正確的r一定小於m，一定放得進64bit
重點是在算 y 的時候就有捨位誤差了，到後面會與正確結果愈差愈遠，
這結果是你所能接受的？
< 以下恕刪 >。
[註一]
我不確定 C/C++ 存浮點數是否一定是使用 ieee754 規格，
但之前研究過的 gcc / vc / bcb 的確是從 ieee754 規格，
< compiler 差別在於「(non)-normalize 之最大、最小特殊值」 >，
52 + 1 裡面的+1指的是在 normalize 時會將首數調成 1,
所以精度才會是 2^-53= 2.2E-16。
--
世界上有種，
將 不可能 轉換為 無限可能 的強大力量，
我稱它為 - 信念。
--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 123.195.165.40