|
FAQ 之 Fortran随机数
2015-04-13 19:44:47 来源：lvlimei @ Fcode 研讨团队 评论：0 点击：
本文详细介绍了 Fortran95 语法规范中的两个随机数函数，并分析了在VF编译器和GFortran编译器上使用的不同。适合新手阅读。
一. 计算机产生的伪随机数
之所以把计算机生成的数称之为伪随机数，是因为它并不是真正的“随机的”，它的周期总是有限的，经过足够的循环次数后，它会重复出现。
关于随机数产生的方法，有大量的研究。而被编译器普遍采用的，是线性同余算法（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
相关热词搜索：随机数
上一篇：Intel Visual Fortran OpenGL 初步
下一篇：堆栈（stack）的故事
|