首页 >

FAQ 之 Fortran随机数

作者:fcode  日期:04-13
来源:lvlimei @ Fcode 研讨团队
一. 计算机产生的伪随机数
之所以把计算机生成的数称之为伪随机数,是因为它并不是真正的“随机的”,它的周期总是有限的,经过足够的循环次数后,它会重复出现。
关于随机数产生的方法,有大量的研究。而被编译器普遍采用的,是线性同余算法(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
常规|工具|专业|读物|
代码|教学|算法|
首页 >
FortranCoder手机版-导航