Re: [問題] 估計pi...
※ 引述《td00414955 (banana)》之銘言：
既然有網友想討論，小弟亂數只是初學，便學「隔空抓葯」、拋磚引玉，
期 神之手 出手相助
先回一下 TsinTa，原po邏輯是真的有誤，而下面兩行並非一定完全相等。
u1=double(s(i))/2147483647;
u2=double(s(i))/2147484000;
在這二行 s(i) 是一樣的沒錯，會相等的機率是 353 / 2147483647,
所以執行起來幾乎是怎麼看都相等。
下面敘述為我的推斷
: int s(int a)
: {
: int s0=100000000,s1;
: if(a==1) return (16807*s0)%2147483647;
: else
: {
: s1=(16807*s(a-1))%2147483647;
: if(s1<0) s1=s1+2147483647;
: return s1;
: }
:}
首先這段碼的 if(a==1) 函意應該是相當於 srand() 設亂數種子的動作，
這段碼推測應是用 同餘法 的方式，
也就是 Yn= (A * Yn-1 + B) MOD C
而上下兩段合起來看，可以知道原po應該是要設
Y0 = s0 = 100000000, 另 A=16807, B=0, C=INT_MAX=2147483647
於是這段程式碼真的是「非常有改善空間」
1. if(a==1) return (16807*s0)%2147483647;
16807 * s0 = 1680700000000 > INT_MAX，所以一定爆
2. s1=(16807*s(a-1))%2147483647;
我猜你這裡要做的應該是 s1 = 16807*s0 % 2147483647;
它並不是用 recursive 做的，只和上一數值有所相依而已。
但會爆的機率仍然太大 (s0: 0~INT_MAX-1, s1*s0 > INT_MAX 就爆了)
即使你發現數值很奇怪，但仍然使用了 if(s1<0) s1=s1+2147483647;
硬把爆掉的數值調回來
3. 我想這也是為什麼 RAND_MAX 大多為 32767 的原因，如果要再大上去的話，
其中一個大缺點是亂數產生慢，因最後所有的計算不是用大數，
就是用 double/float 去運算，包含了加法和乘法及模數，
也就是說，這裡不能用 % , 只能用 fmod。
綜合以上三點說明，我想以下程式碼應是你想要的，可優化修改空間大。
int s()
{
static bool first_call = true;
static double s0=0, s1=0;
if(first_call) {
s1 = ceil( fmod( A*S0+B, C) + 0.5); // fmod後四捨五入
first_call=false;
}
else {
s1 = ceil( fmod( A*s0+B, C) + 0.5); // fmod後四捨五入
}
s0=s1;
return int(s1);
}
如果這份 project/ hw 的重點只是蒙地卡羅法求 pi 而不是寫亂數產生器,
那麼比較建議是用內建的亂數產生器去做就好。
[Lemma1]
手邊是用 vc, 編過程式碼測過亂數產生器週期， vc/dev-c 雖 RAND_MAX = 32767，
但亂數週期的確為 INT_MAX (此後開始重覆)，所以若用蒙地卡羅法，
測試次數別超過 INT_MAX 次。
[Lemma2]
亂數產生器測試週期，最後我是記前 n(5~10) 個亂數值，再一直取下去，
最後發現連續 n 個亂數值都一樣時，就代表重覆。除非那個亂數產生器的週期小於 n
又另當別論 (所以最好再寫一個大數演算法的遞增去算週期比較準)。
[Lemma3]
真的不是必要的話，建議不要自己寫亂數產生器，寫出來是一回事；
怎麼測試產生器的品質 (均勻度、週期等) 又是另一回事，通常真的是數學系在搞的。
如果真的要搞大範圍、均勻度、速度等各方面協調都好的話
http://en.wikipedia.org/wiki/Mersenne_twister
試試這個吧！！
--
YouLoveMe() ? LetItBe() : LetMeFree();
--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 180.177.72.67
※ 編輯: tropical72 來自: 180.177.72.67 (04/09 02:48)
推
04/09 02:56, , 1F
04/09 02:56, 1F
→
04/09 09:10, , 2F
04/09 09:10, 2F
→
04/09 10:35, , 3F
04/09 10:35, 3F
→
04/09 13:45, , 4F
04/09 13:45, 4F
→
04/09 21:24, , 5F
04/09 21:24, 5F
→
04/09 21:42, , 6F
04/09 21:42, 6F
→
04/09 21:47, , 7F
04/09 21:47, 7F
討論串 (同標題文章)