線性同餘法詳解

一、基礎概念

線性同餘法（Linear congruential generator，簡稱LCG）是一種偽隨機數生成器。其具體思想是以一種線性的方式對一個起始值進行遞推計算，以得到一個序列化的偽隨機數序列。LCG可以形式化地表示為：Xn+1 = (aXn + c) mod m，其中Xn+1是下一個偽隨機數，Xn是當前偽隨機數，a、c、m是先驗設定好的參數，mod表示模運算。

LCG能夠快速地生成大量的偽隨機數，因此被廣泛地應用於模擬、密碼學等領域，同時也是許多編程語言中的隨機數生成器的核心算法。

二、算法細節

在LCG中，參數a、c和m需要精心設計，否則可能導致隨機數生成效果不佳。以下是一些需要注意的細節：

1.參數m

參數m應該為一個較大的質數，以確保隨機數循環周期的長度較長，保證生成的偽隨機數足夠分散。同時，m也需要足夠小，否則會影響算法執行的效率。

const int m = 2147483647; // 2^31 - 1 

2.參數a

參數a需要滿足以下條件：

1）a和m互質，以保證偽隨機數的分散性。

2）a-1能夠被m的所有質因子整除，以避免偽隨機數循環周期過短。

const int a = 16807; // 隨意選取的常數

3.參數c

參數c一般取0或1，不會對偽隨機數的質量產生影響。如果使用C++標準庫中的rand函數，其使用的是a=1103515245, c=12345, m=2147483647的LCG。

const int c = 0; // 常數c一般取0或1

三、優化方法

LCG是一種簡單但不一定高效的偽隨機數生成器，其生成的偽隨機數分佈不夠均勻，容易被攻擊者破解。為了彌補LCG的不足，可以採取以下優化方法：

1.超前迭代法

超前迭代法是一種通過多次遞歸來產生更高質量隨機數的方法，它可以顯著提高LCG的效率和隨機性。在超前迭代法中，可以使用兩個LCG，其中一個作為「控制器」，控制另一個的參數。這種方法可以產生出高質量的隨機數序列，但需要佔用相當大的空間。

// 線性同餘法優化：超前迭代法
class Random {
public:
    Random() {
        // 控制器LCG
        c_r = c0_r;
        a_r = a0_r;
        m_r = m0_r;
        x_r = time(nullptr);
        for (int i = 0; i < 10; i++) {
            next_r();
        }
    }
    int next_r() {
        int x = x_r;
        // LCG1作為控制器
        x_r = (a0_r * x_r + c0_r) % m0_r;
        // LCG2作為普通LGC
        return (a_r * x + c_r) % m_r;
    }
private:
    const int c0_r = 1;
    const int a0_r = 16807;
    const int m0_r = 2147483647;
    int c_r, a_r, m_r, x_r;
};

2.多項式同餘法

多項式同餘法是一種比LCG更高效的偽隨機數生成算法。它使用多項式函數來代替線性函數，進一步提升了隨機數生成效率和質量。不過，這種方法需要較大的空間，在實際應用時需要權衡空間和時間的開銷。

// 多項式同餘法
class Random {
public:
    Random() {
        for (int i = 0; i < degree_r; i++) {
            a_r[i] = rand();
        }
    }
    int next_r() {
        int x = a_r[0];
        for (int i = 1; i < degree_r; i++) {
            x = (a_r[i] + x * seed_r) % m_r;
        }
        for (int i = 0; i < degree_r-1; i++) {
            a_r[i] = a_r[i+1];
        }
        a_r[degree_r-1] = x;
        return x;
    }
private:
    const int m_r = 2147483647;
    const int seed_r = 212885;
    const int degree_r = 10;
    int a_r[degree_r];
};

四、應用實例

以下是一個應用LCG生成模擬數據的實例。在這個例子中，使用LCG和超前迭代法的組合產生了10000個隨機數，用於繪製一個符合正態分佈的數據圖表。

#include <iostream>
#include <vector>
#include <cmath>
#include <fstream>
#include <string>
#include <ctime>

class Random {
public:
    Random() {
        c_r = c0_r;
        a_r = a0_r;
        m_r = m0_r;
        x_r = time(nullptr);
        for (int i = 0; i < 10; i++) {
            next_r();
        }
    }
    int next_r() {
        int x = x_r;
        x_r = (a0_r * x_r + c0_r) % m0_r;
        return (a_r * x + c_r) % m_r;
    }
private:
    const int c0_r = 1;
    const int a0_r = 16807;
    const int m0_r = 2147483647;
    int c_r, a_r, m_r, x_r;
};

int main() {
    Random random;
    std::vector<double> nums;
    for (int i = 0; i < 10000; i++) {
        nums.push_back(random.next_r() / static_cast<double>(random.m0_r));
    }

    // 使用Box-Muller變換將均勻分佈轉換為正態分佈
    std::vector<double> nums2;
    for (int i = 0; i < 10000; i += 2) {
        double u1 = nums[i];
        double u2 = nums[i + 1];
        double z1 = sqrt(-2 * log(u1)) * cos(2 * M_PI * u2);
        double z2 = sqrt(-2 * log(u1)) * sin(2 * M_PI * u2);
        nums2.push_back(z1);
        nums2.push_back(z2);
    }

    // 統計數據
    double sum = 0, avg = 0, stdev = 0;
    for (auto num : nums2) {
        sum += num;
    }
    avg = sum / nums2.size();
    for (auto num : nums2) {
        stdev += pow(num - avg, 2);
    }
    stdev = sqrt(stdev / (nums2.size() - 1));

    // 輸出結果到文件
    std::ofstream out("output.txt");
    for (auto num : nums2) {
        out << num << std::endl;
    }
    out.close();

    std::cout << "avg = " << avg << std::endl;
    std::cout << "stdev = " << stdev << std::endl;

    return 0;
}