之所以把计算机生成的数称之为伪随机数,是因为它并不是真正的“随机的”,它的周期总是有限的,经过足够的循环次数后,它会重复出现。
关于随机数产生的方法,有大量的研究。而被编译器普遍采用的,是线性同余算法(LCG)。
二. Fortran90 函数 Random_Seed 和 Random_Number
在Fortran90之前,各家编译器扩展有不同的随机数生成函数。例如 seed 和 random,ran 或 rand 之类的。
导致不同编译器的写法不同,通用性不强。
而在 Fortran90 以后,语法规范引入了两个标准的函数用来产生随机数。它们就是 random_seed 和 random_number( 通常这两个函数需配合使用 )
1. random_seed
刚才说了,计算机产生的是伪随机数,只能根据一个初始值,计算出一个序列来,看起来,这个序列毫无规律。但实际上,他依然依赖于特定的计算,而这样的计算结果,是唯一存在的。
于是,要产生每次不一样的随机数。必须有不同的初始值,这就是 “种子”
random_seed 用来获取或设置新的种子。不同的编译器,都会有一套“默认”的种子。我们可以在一开始获取它,也可以用自己的种子来覆盖它。
注意,这个函数在整个程序里,通常只调用一次!
call random_seed( size , get , put )
这个函数一共三个参数(一般一次只能选择一个参数使用,而不能三个一起使用):
size 是输出参数,会返回种子的大小。在获取或设置新的种子时,输入输出的种子数组,必须与 size 的大小一样。
比如: call random_seed( size = n ) ,执行后,n 的值就是种子的大小。
get 是输出参数,用来返回当前的种子。必须是 integer 类型的数组(数组大小取决于 size),通常很少使用它。
put 是输入参数,用来设置新的种子。
所以,常规的设置种子的办法是:
1. call random_seed( n ) 获取种子的大小
2. allocate() 分配种子
3. 对种子设置一定的值,例如利用时间函数设置种子。(这样不同的时间可以获得不同的随机数)
4. call random_seed( put ) 设置新的种子。
然而,在 IVF 编译器上,这一切都比较容易。因为它规定,只要random_seed不加入任何参数,则自动用时间设置种子。所以这个操作可以简化为 call random_seed()
2. random_number
这个函数是用来产生随机数的,当我们需要产生时,调用它就可以了。整个程序可以不限次数地调用它。
call random_number( x )
这里的 x 必须是 real 类型,可以是单变量,也可以是数组。调用后,x 的值变为当前的(伪)随机数。
三. 在GFortran和IVF上的区别
在 IVF / CVF 编译器上,用这两个函数产生随机数(每次都不相同)的方法非常简单。
Program www_fcode_cn Implicit None real :: x(3) = 0.0 , y(2) = 0.0 , a = 0.0 integer :: i call random_seed() !// 只调用一次,读者可尝试去掉此行并多次运行 Do i = 1 , size(x) call random_number( x(i) ) !// 可以循环调用 End Do call random_number(y) !// 也可以一次生成一个数组 call random_number(a) !// 单变量也可以 write(*,*) a write(*,*) x write(*,*) y End Program www_fcode_cn
读者可以尝试去掉random_seed,并多次运行,会发现生成的随机数每次都一样。这是因为编译器默认的种子是不变的。
而在 GFortran 或其他编译器上,这个过程可能就要复杂一些。以下代码演示了如何使用时间和PID产生随机数的过程。(读者也可以自行设计自己的随机数及种子产生方法)
Module Random_Mod Implicit None !// 使用内部函数获得int64的Kind值 Integer , parameter :: int64 = Selected_Int_Kind( 10 ) contains Subroutine Init_Random_Seed() integer :: ised , i , pid integer(int64) :: t integer , allocatable :: sed(:) call random_seed( size = ised ) !// 获得种子大小 allocate( sed(ised) ) !// 分配种子 call system_clock(t) !// 获得时间 pid = getpid() !// 获得处理器ID t = ieor(t, int(pid, kind(t))) !// 用 pid 和日期做XOR运算 do i = 1, ised sed(i) = lcg(t) !// 用线性同余计算种子 end do call random_seed( put=sed ) !// 给定种子 End Subroutine Init_Random_Seed Function lcg(s) !// 线性同余算法 Linear congruential generator integer :: lcg integer(int64) :: s if (s == 0) then s = 104729 else s = mod(s, 4294967296_int64) end if s = mod(s * 279470273_int64, 4294967291_int64) lcg = int(mod(s, int(huge(0), int64)), kind(0)) End Function lcg End Module Random_Mod Program www_fcode_cn use Random_Mod Real :: x(8) call Init_Random_Seed() call random_number(x) write(*,*) x End Program www_fcode_cn
读者如果觉得麻烦,或者读懂这段代码比较困难,可以直接把 Module Random_Mod 拷贝过去直接使用。只修改主程序既可。(注意 Init_Random_Seed 也只调用一次)
四. 产生区间上的随机分布和高斯(正态分布)
前面我们讨论的函数,都只能生成 0 只 1 之间的随机数。
那么如果我们需要产生 a 到 b 区间的随机数怎么办呢?此时,我们可以用到线性变换 ax + b
call random_number(x) !// 随机数种子部分忽略不写 a = 1.0 !// 产生从 1.0 b = 100.0 !// 到 100.0 之间的随机数 x = abs(b-a) * x + min(a,b) write(*,*) x
而如果希望得到高斯分布的随机数,可参考本站文章:/code_prof-33-1.html