1. 基本觀念

線性同餘法基本上只有一個公式,X(n+1) = ( a * X(n) + b ) mod c

但並非所有 a, b, c 都適用。以 a=5, b=0, c=10 為例,假定 X(0)=3 代入

X(1) = ( 5 * 3 + 0 ) mod 10 = (15 + 0 ) mod 10 = 5
X(2) = ( 5 * 3 + 0 ) mod 10 = (15 + 0 ) mod 10 = 5

這完成沒辦法產生亂數,0~9 裡面,只有 5 會不停重覆產生。

 

2. 常數選擇

 重點便是在 a, b, c 應選何值、怎麼選為佳。有些規則可供參考,下列三項規則取自 TAOCP

X(n+1) = ( a * X(n) + b ) mod c

(1) 乘數 a 與模數 c 互質。
(2) 對於整除 c 之每個質數 p,a-1 抑可整除 p。
(3) 若 c 為 4 之倍數, b 也為 4 之倍數。

在一些文獻記載中,有提到幾個 a, b, c 之選取,可達到全週期之產生器

a = 16807 , b=0, c=2147483647
a = 630360016, b=0, c=2147483647

上述的 b 都為 0,而在文獻上,LCG 確實還分二種,一種是 b=0的,叫 Lehmer random number generator,或稱 Multipicative LCG;另一種是 b 不為 0 的,叫 Prime Moduls Multiplicative LCG,PMMLCG。


3. 程式碼實做

 實做程式碼參考如下

 

/*
 * filename : LCG.c
 * author   : edisonx / edison.shih
 * compiler : Visual C++ 2008
 *
 */

#include <stdio.h>
#include <time.h>
long long seed;

void LCGSrand(int x){ seed = (long long )x; }

int LCGRand()
{
    const long long a=16807LL; // or a =  630360016LL
    const long long b=0LL;
    const long long c=2147483647LL;
    seed = (a*seed + b ) % c;
    return (int)seed;
}

int main()
{
    int a, b, c, d;
    int x1=0, x2=0, x3=0, x4=0;
    unsigned period=0;
    LCGSrand(time(NULL));
    a=LCGRand(), b=LCGRand(), c=LCGRand(), d=LCGRand();
    do{
        x1 = x2;
        x2 = x3;
        x3 = x4;
        x4=LCGRand();
        ++period;
    }while(x1!=a && x2!=b && x3!=c && x4!=d);
    printf("period: %u\n", period);
    return 0;
}

 


4. RAND_MAX 與週期探討

 這裡有一些問題可探討。

1. 為何要用 long long ?
2. 為何 Visual C++ 系列,RAND_MAX 是 32767,但實測週期 (可能官方有文件紀錄) 是 INT_MAX-1?

第一個問題是溢位問題。當 seed "不小" (若 a=16807,即大於 127773) 時,便會造成 ov 問題。

第二個問題應要拆成二部份來看。Visual C++ 在實做時有用一點修正技巧,它確實也是用 LCG 實做 rand(),但最後在回傳時,只取其中15 bits 而已,故範圍是 0~32767。而這 15 bits,不表示一定是取最低的 15 bits。在 這個 wiki 網頁上,甚至可以知道一些 compiler 在做 rand,若採用 LCG 時之常數選擇,其中裡面便提到了,Visual C++ 取的是

a=214013, b=2531011, c=2^31,傳回 bit30~bit16

模擬該網頁的常數寫法,和 Visual C++ 2008 之 rand() 結果比較,事實上結果並不完全相等 (前幾項的確是很準,後面就幾乎不準),這部份應是 VC 那裡有再用一些手法套入。

這也是為什麼,調用 rand() ,第一個重覆的數字間隔代數,並不等於 RAND_MAX 之因。即

a = rand();
while(rand()!=a) ++times;

在 vc 下, times 不會等於 RAND_MAX ,一些書上會說,要先將0~c-1 產生完一次之後,才會有重覆的情況發生,但若取部份 bit,這情況並不適用,由此例可知。

 


 

5. 找其他 a, b, c 常數

上述程式碼中,運算一直是用 long long int 在運算,速度慢!有辦法可以進行較快速的運算嗎?

有一個簡單的方法,但必須犧牲一個重要的條件:週期會變短 。再看一次公式:

X(n+1) = ( a * X(n) + b ) mod c

若將 c 選在小於 2^16 以下,便可使用 int/unsigned int 完成,但相對的全週期 (如果有的話) 就只剩下 c-1 而已。根據上述的三項規則

X(n+1) = ( a * X(n) + b ) mod c

(1) 乘數 a 與模數 c 互質。
(2) 對於整除 c 之每個質數 p,a-1 抑可整除 p。
(3) 若 c 為 4 之倍數, b 也為 4 之倍數。

額外提,我認為 c 直接選 梅森質數 可能會比較方便一點,如 8191(13bits)、131071(17bits)、524287(19bits)、2147483647 (31bits)、2305843009213693951 (61bits),要找梅林質數的話,建議用「盧卡斯-萊默檢驗法」。

那有使 LCG 週期長、速度又快的方法嗎?我想這部份應直接去看 compiler 出來 asm code 長怎樣,推斷已「不簡單」。


arrow
arrow
    全站熱搜

    edisonx 發表在 痞客邦 留言(0) 人氣()