[RAND] Linear congruential generator (線性同餘法,LCG)－Edison.X. Blog

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 之倍數。