基于Python的随机数生成器

• 希望大家测试
• 希望大家测试
• 交流+qq 3052408134

• 随机数 
• 伪随机数
• 准随机数

import numpy as np
import matplotlib.pyplot as plt
import csv
import sys
import time
import copy

#uniform distribution
class RG():
    def __init__(self):
        pass 


    #Square random number generator
    def SRG(self,Len=10,seed=42432):
        x = seed
        assert(type(x)==int),"seed should be integer!!!"
        k = len(str(x))/2
        U = np.zeros(Len)
        for i in range(Len):
            x = int(x**2/(10**k))%(10**(2*k))
            U[i] = x/(10**(2*k))
        return U

    #Product medium random number generator
    def PMG(self,Len=10,seed1=123331,seed2=432133,k=3):
        assert(len(str(int(seed1)))==len(str(int(seed2)))==2*k),"seed1 and seed2 should be 2k bits non-negative integer"
        U = np.zeros(Len)
        X = np.zeros(Len+2)
        X[0] = seed1
        X[1] = seed2
        for i in range(Len):
            X[i+2] = (X[i]*X[i+1]/(10**k))%(10**(2*k))
            U[i] = X[i+2]/10**(2*k)
        return U

    
    #Constant multiplier generator
    def CMG(self,Len=10,seed1=12333091,M=14321133,k=4):
        assert(len(str(int(M)))==len(str(int(seed1)))==2*k),"M and seed1 should be 2k bits non-negative integer"
        x = int(seed1)
        U = np.zeros(Len)
        for i in range(Len):
            x = int(M*x/10**k)%(10**(2*k))
            U[i] = x/10**(2*k)
        return U



    #linear fibonacci generator
    def LFG(self,Len=10,seed1=123331,seed2=432132133,m=231):
        assert(type(seed1)==type(seed2)==type(m)==int),"seed should be integer!!!"        
        U = np.zeros(Len)
        X = np.zeros(Len+2)
        X[0] = seed1
        X[1] = seed2
        for i in range(Len):
            X[i+2] = (X[i]+X[i+1])%m
            U[i] = X[i+2]/m
        return U
    

    #simple linear congruential generators
    def LCG(self,Len=10,a=3,c=1,seed=1,m=8):
        U = np.zeros(Len)
        x = seed
        for i in range(Len):
            x = (a*x+c)%m
            U[i] = x/m
        return U

    #multiple recursive generators
    def MRG(self,Len=10,r1=123456789,r2=987654321,a=428773,m=2**32):
        U = np.zeros(Len)
        for i in range(Len):
            r1,r2 = (a*r1+r2)%m,r1
            U[i] = (a*r1+r2)%m/m
        return U
    
    #add-with-carry
    def ACG(self):
        pass

    #subtract-with-borrow
    def SWBG(self):
        pass

    #multiply-with-carry
    def MWCG(self):
        pass
    
    #matrix congruential generators
    def MatrixCG(self):
        pass

    #混(组)合同余法
    #Mixed congruence method
    def MCMG(self,Len=10,a=233,b=41,seed=61,m=28):
        U = np.zeros(Len)
        x = seed
        for i in range(Len):
            x = a*x+b
            U[i] = (x%m)/m
        return U

    #乘同余法
    #Multiplication  congruence method
    def MuCMG(self,Len=10,a=2343,seed=7561,m=2844534):
        U = np.zeros(Len)
        x = seed
        for i in range(Len):
            x = a*x
            U[i] = (x%m)/m
        return U

    
    #加同余法
    #Additive congruence
    def ACMG(self,Len=10,k=4,Y0=np.array([17512435,22345,33422,68974]),m=489636696231):
        assert(len(Y0)==k),"k is not equal with length Y"
        Y = np.zeros(k+Len)
        U = np.zeros(Len)
        try:
            Y[:k] = Y0
        except:
            print("loaded failure!!!")
            sys.exit()
        for i in range(k,Len+k):
            Y[i] = sum(Y[i-k:i])
            U[i-k] = Y[i]%m/m
        return U


    #inversive congruential generators
    def ICG(self):
        pass

    #feedback shift register generators
    def FSRG(self):
        pass


    #逆变法
    #Inverse Transform method
    def ITM(self,Len=10):
        #print("Warning!using sin(x)^2 [0,pi/2]as cdf")
        cdf = lambda x:np.sin(x)**2
        pdf = lambda x:np.sin(x)*np.cos(x)*2
        n_samples=max(1000,10*Len)
        def Inverse_Transform(Y):
            nonlocal n_samples
            nonlocal cdf
            x = [i/n_samples for i in range(n_samples)]
            y = [cdf(i*np.pi/2) for i in x]
            ABS = [abs(Y-i) for i in y]
            return x[ABS.index(min(ABS))]
        X = self.LCG(Len)
        U = np.array([Inverse_Transform(i) for i in X])
        return U
        
    #Accept-reject method
    def ARM(self,Len=10):
        f = lambda x:1
        M = np.sqrt(2*np.exp(1)/np.pi)
        def ARM_f():
            while True:
                x = self.LCG()[0]
                y = self.LCG()[0]
                if y <= f(x)/M:
                    return x
        U = np.zeros(Len)
        for i in range(Len):
            U[i] = ARM_f()
        return U

    
################################################################
#############Performance Test###################################
################################################################

    #参数检验
    def parameter_test(self,f,Len=1000,PLOT=True,SAVE=True):
        #ER1 = 1/2
        #ER2 = 1/3
        #ES2 = 1/12

        #DR1 = 1/(12N)
        #DR2 = 4/(45N)
        #DS2 = 1/(180N)

        U = f(Len=Len)
        MIN = min(U)
        MAX = max(U)
        MEAN = sum(U)/len(U)
        S = (sum((U-MEAN)**2)/(len(U)-1))**0.5
        RANG = MAX - MIN
        MEDIAN = sorted(U)[int(len(U)/2)]
        SKWNESS = sum(((U-MEAN)/S)**3)/len(U)
        KURTOSIS = sum(((U-MEAN)/S)**4)/len(U)-3

        INFO = {"MIN":MIN,
                "MAX":MAX,
                "MEAN":MEAN,
                "S":S,
                "RANG":RANG,
                "MEDIAN":MEDIAN,
                "SKWNESS":SKWNESS,
                "KURTOSIS":KURTOSIS}
        if SAVE:
            f = open(str(round(time.time()))+".txt","w")
            f.write("num: "+str(Len)+"\n")
            f.write(str(INFO))
            f.write("\n#################")
            f.close()
            print(INFO)

        return INFO
    
    #组合规律检验
    def Poker_test(self,f,Len=1000,PLOT=True,SAVE=False):
        U = f(Len=Len)
        P = np.array([0.02,0.1703,0.4205,0.3195,0.0697])
        M = divmod(len(U),8)[0]*P

        Poker = np.array(list((10*np.round(U-0.04999999999,1)%8).astype(int)))   
        Count = np.array([len(set(Poker[i*8:(1+i)*8]))  for i in range(divmod(len(Poker),8)[0])])
        Poker = np.zeros(5)
        for i in Count:
            if i<=3:
                Poker[0] += 1
            elif i==4:
                Poker[1] += 1
            elif i==5:
                Poker[2] += 1
            elif i==6:
                Poker[3] += 1
            elif i>=7:
                Poker[4] += 1
        print("M",M)
        print("Poker",Poker)
        T = sum((M-Poker)**2/M)
        print("T",T)
        chi_4_95 = 9.4877
        POKER = chi_4_95<T
        return {"POKER":POKER}

        
    def minimal_statistical_test(self,f,Len=1000,PLOT=True,SAVE=False):
        List = f(Len=Len)
        ave = np.zeros(Len)
        for i in range(1,Len):
            ave[i] = sum(List[:i])/i
        
        if PLOT:
            plt.plot(range(Len),ave)
            plt.title("$minimal\,statistical\,test$")
            plt.pause(0.01)
        if PLOT and SAVE:
            plt.savefig("minimal_statistical_test "+str(Len)+".jpg")
            plt.cla()            
            
    def Hist_test(self,f,Len=1000,Bins=50,Density=True,SAVE=False):
        List = f(Len=Len)
        plt.title("$Hist\,photo$")
        plt.hist(List,bins=Bins,density=Density)
        plt.pause(0.01)
        if SAVE:
            plt.savefig("Hist_test "+str(Len)+".jpg")
            plt.cla()  

    def PI_test(self,f,square_a=1,Len=1000,TIME=True,PLOT=True,SAVE=False):
        if TIME:
            start = time.time()
        Location = f(Len=Len*2).reshape((Len,2))
        Outside = list(filter(lambda x:x[0]**2+x[1]**2>square_a**2,Location))
        Inside  = list(filter(lambda x:x[0]**2+x[1]**2<=square_a**2,Location))
        Pi = len(Inside)/Len*4*square_a**2
        print("Pi:",Pi)
        if TIME:
            end = time.time()
            print("used time",end-start)
        if PLOT:
            fig, ax = plt.subplots()
            plt.scatter([i[0] for i in Inside],[i[1] for i in Inside],c="red",s=0.5)
            plt.scatter([i[0] for i in Outside],[i[1] for i in Outside],c="blue",s=0.5)
            plt.title("$PI_{test}\,num:"+str(Len)+"\,Pi:"+str(Pi)+"$")
            circle = plt.Circle((0,0),square_a,color="green",fill=False)
            ax.add_artist(circle)
            plt.pause(0.01)
        if PLOT and SAVE:
            plt.savefig("PI_test "+str(Len)+".jpg")
            plt.cla()  

################################################################
#############Quasirandom Numbers################################
################################################################

    def Decomposition_prime_factor(self,N):
        T = []
        M = copy.deepcopy(N)
        T.append(1)
        T.append(M)
        while M>1:
            for i in range(2,M+1):
                if M%i == 0:
                    M = M//i
                    if i not in T:
                        T.append(i)
                    break
        return T


    def Halton(self,Len=10,m=15,BASE=[3,5]):
        newL = divmod(Len,len(BASE))[0]+1
        U = np.zeros((newL,len(BASE)))
        for i in range(newL):
            for j in range(len(BASE)):
                U[i][j] = self.convert_decimal(self.get_decimal(m,BASE[j]),BASE[j])
                m += 1
        return U.reshape(newL*len(BASE))[:Len]

################################################################
#############Tools##############################################
################################################################
    def Decomposition_prime_factor(self,N):
            T = []
            M = copy.deepcopy(N)
            while M>1:
                for i in range(2,M+1):
                    if M%i == 0:
                        M = M//i
                        if i not in T:
                            T.append(i)
                        break
            return T

    def get_decimal(self,number,base):
        if not(2 <= base <= 36):
            return "Base value must be between 2 and 36."
        elif not isinstance(number, int) or number < 0:
            return "Input value must be a positive integer."
        else:
            if number == 0:
                return "0"
            else:
                output = ""
                while number != 0:
                    digit = number % base
                    if digit < 10:
                        output = str(digit) + output
                    else:
                        output = chr(55 + digit) + output
                    number = number // base
                return "0."+output[::-1]

    def convert_decimal(self,num,base):
        decimal = 0
        fraction = num.split('.')[1]
        length = len(fraction)
        for i in range(length):
            decimal += int(fraction[i]) * (base ** -(i+1))
        integer = int(str(num).split('.')[0])
        i = 0
        while integer:
            decimal += (integer % 10) * (base ** i)
            integer //= 10
            i += 1
        return decimal

################################################################
#############Main Function######################################
################################################################    
if __name__ == "__main__":
    r = RG()
    f = r.Halton
    #r.Hist_test(f)
    #r.minimal_statistical_test(f)
    r.parameter_test(f)
    r.PI_test(f,Len=1000)