Python 实现随机数生成算法（线性同余和梅森旋转）

在Python中生成随机数, 一般会调用random模块, 但random模块内也有自己的算法实现。
如何设计自己的算法呢?本文将介绍两大算法线性同余算法和梅森旋转算法的程序实现。

1.线性同余算法

线性同余算法与梅森旋转算法相比相对简单。生成随机数, 首先要有一个随机数种子seed, 然后根据这个种子, 推导出剩余的数字, 也就是伪随机数。
一种方法是将种子经过一定的变换, 再除以一个数,
得到的余数作为新的种子。这种方法也称线性同余算法。

程序实现

from time import perf_counter

class LinearCongruence:
    def __init__(self,seed,k,b,p):
        self.seed=seed
        self.k=k;self.b=b;self.p=p
    def random(self):
        self.seed = (self.seed*self.k + self.b) % self.p
        return self.seed

def count(lst,max_value,d):
    # 统计lst中每个范围内数值的出现次数,d为组距
    result=[0]*math.ceil(max_value/d) # ceil为向上取整
    for i in lst:
        result[int(i//d)]+=1
    return result

seed = 1.7
k=5.6;b=3.5
p=1;d=0.001

l=[] # 存储生成的随机数
cnt=10000
lc=LinearCongruence(seed,k,b,p)
for j in range(cnt):
    l.append(lc.random())

自制算法的检验

对随机数算法的要求是随机数的重复周期很长, 且在生成的区间内分布均匀。
首先生成包含10000个随机数的列表，绘制随机数从0至1的数值分布直方图，以检验随机数的数值分布是否均匀。

import math
import matplotlib.pyplot as plt
from random import random as random2

class LinearCongruence:
    def __init__(self,seed,k,b,p):
        self.seed=seed
        self.k=k;self.b=b;self.p=p
    def random(self):
        self.seed = (self.seed*self.k + self.b) % self.p
        return self.seed

def count(lst,max_value,d):
    # 统计lst中每个范围内数值的出现次数,d为组距
    result=[0]*math.ceil(max_value/d) # ceil为向上取整
    for i in lst:
        result[int(i//d)]+=1
    return result

plt.rcParams['font.sans-serif'] = ['SimHei'] # 用于正常显示中文
seed = 1.7
k=5.6;b=3.5
p=1;d=0.001

l=[]
cnt=10000
lc=LinearCongruence(seed,k,b,p)
for j in range(cnt):
    l.append(lc.random())

# 绘制l从0至1的数值分布直方图，检验随机数的数值分布是否均匀
l=count(l,p,d)
plt.bar(range(len(l)),l)
plt.title("自制的线性同余算法")
plt.show()

l2=[random2() for i in range(cnt)]
# 绘制l2的数值分布直方图
l2=count(l2,p,d)
plt.bar(range(len(l2)),l2)
plt.title("Python内置random模块")
plt.show()

运行结果如下（注：图中横坐标的0 ~ 1000应为0 ~ 1(抱歉)，纵坐标为在每个0.001的区间（如[0.002, 0.003) ）中随机数的出现次数）：

作为对比，下面是Python内置的random模块：

从中可以发现, 前述自制的随机数生成算法生成的值比较均匀, 比较好地模拟了random模块中的随机数生成。
除此之外，还需要计算随机数的数学期望和方差：

# 生成l和l2
seed = 1.7
k=5.6;b=3.5
p=1;d=0.001
l=[]
cnt=500
lc=LinearCongruence(seed,k,b,p)
for j in range(cnt):
    l.append(lc.random())
l2=[random2() for i in range(cnt)]

e = sum(l)/len(l) # 数学期望
a = e # 平均值
d = sum((x-a)**2 for x in l2) /len(l2)
print("自制随机数的期望和方差:",e,d) # 结果为 0.50298433624 0.080202207010

e2 = sum(l2)/len(l2)
a2 = e2
d2 = sum((x-a2)**2 for x in l2) /len(l2)
print("内置random函数的期望和方差:",e2,d2) # 结果为 0.51050220327 0.08014568868

2.梅森旋转算法

由于线性同余算法重复周期短、低维分布、初始种子敏感性、具有线性关系等缺陷，Python、C++等语言的内置函数使用的算法不是线性同余算法，而是梅森旋转算法。

梅森旋转（Mersenne Twister）算法是一种伪随机数生成算法，以梅森素数命名。梅森旋转算法的主要特点是周期长（219937 - 1），并且在统计测试中表现良好。

梅森旋转算法的原理基于线性同余发生器（LCG）和位运算。它使用一个624维的状态数组来存储随机数生成器的内部状态。在初始化时，需要提供一个种子值，然后通过一系列的位运算和状态数组的扭曲操作来生成随机数。

具体来说，梅森旋转算法的状态数组初始化后，每次需要生成随机数时，会从状态数组中取出一个值，对其进行一系列的位运算操作，然后更新状态数组中的值。当状态数组中的值用尽时，会进行扭曲操作，重新生成一组新的状态数组值。这样不断循环，就能够生成一系列的伪随机数。

程序实现

import math
import matplotlib.pyplot as plt

class MersenneTwister:
    def __init__(self, seed):
        # seed的范围为-2147483648到2147483647
        self.index = 624  # 初始化索引
        self.mt = [0] * 624  # 初始化状态数组
        self.mt[0] = seed  # 使用种子初始化状态数组的第一个元素
        for i in range(1, 624):  # 初始化状态数组的其余元素
            self.mt[i] = 0xFFFFFFFF & (1812433253 * (self.mt[i - 1] ^ (self.mt[i - 1] >> 30)) + i)

    def extract_number(self):
        if self.index >= 624:  # 如果索引超出范围，进行扭曲操作
            self.twist()

        y = self.mt[self.index]  # 取出状态数组中的一个元素
        y = y ^ (y >> 11)  # 对取出的元素进行位运算
        y = y ^ ((y << 7) & 0x9D2C5680)  # 对取出的元素进行位运算
        y = y ^ ((y << 15) & 0xEFC60000)  # 对取出的元素进行位运算
        y = y ^ (y >> 18)  # 对取出的元素进行位运算

        self.index += 1  # 更新索引
        return 0xFFFFFFFF & y  # 返回生成的随机数
    def random(self):
        # 返回0-1的随机数
        return self.extract_number() / 0xFFFFFFFF

    def twist(self):
        for i in range(624):  # 扭曲操作
            y = (self.mt[i] & 0x80000000) + (self.mt[(i + 1) % 624] & 0x7fffffff)
            self.mt[i] = self.mt[(i + 397) % 624] ^ (y >> 1)
            if y % 2 != 0:
                self.mt[i] = self.mt[i] ^ 0x9908B0DF
        self.index = 0  # 重置索引

def count(lst,max_value,d):
    # 统计lst中每个范围内数值的出现次数,d为组距
    result=[0]*math.ceil(max_value/d) # ceil为向上取整
    for i in lst:
        result[int(i//d)]+=1
    return result

def main():
    plt.rcParams['font.sans-serif'] = ['SimHei'] # 用于正常显示中文

    seed = 5489  # 设置随机数种子
    mt = MersenneTwister(seed)
    lst=[];cnt=10000
    for i in range(cnt):
        lst.append(mt.random())
    lst=count(lst,1,0.001)
    plt.bar(range(len(lst)),lst)
    plt.title("梅森旋转算法")
    plt.show()

if __name__ == "__main__":
    main()

再次按照前文方法检验算法，绘制随机数从0至1的数值分布直方图，以检验随机数的数值分布是否均匀。

此外，按前文所述的方法检验，程序生成随机数的期望和方差分别为0.5068390900592408和0.07965702405546546。

3.总结

前面所述的线性同余算法和梅森旋转算法是生成随机数算法的一种，
除此之外，还有随机数表法等算法。随机数表法则是在已有的大随机数表中按一定的规律取数字。
文章到此结束，欢迎点赞、收藏！