求助：如何修改Fortran标记泊松簇过程随机模拟器的随机数种子？

嘿，你猜的完全没错！问题根源就是伪随机数生成器的种子被固定死了——你贴的两个随机数生成子程序里，初始种子都是用data语句硬编码的（比如第一个里的data ir/574289/，第二个里的data ix,iy,iz /1992,1111,1151/），每次程序运行都会从同一个起点生成随机数，结果自然重复。

先看你提供的两个随机数子程序代码：

subroutine pseud0(r)
c generation of pseudo-random numbers
c data ir/584287/
data ir/574289/
ir=ir*48828125
if(ir) 10,20,20
10 ir=(ir+2147483647)+1
20 r=float(ir)*0.4656613e-9
return
end

subroutine pseudo(random)
c wichmann+hill (1982) Appl. Statist 31
data ix,iy,iz /1992,1111,1151/
ix=171*mod(ix,177)-2*(ix/177)
iy=172*mod(iy,176)-35*(iy/176)
iz=170*mod(iz,178)-63*(iz/178)
if (ix.lt.0) ix=ix+30269
if (iy.lt.0) iy=iy+30307
if (iz.lt.0) iz=iz+30323
random=mod(float(ix)/30269.0+float(iy)/30307.0+ &
float(iz)/30323.0,1.0)
return
end

下面给你针对两个子程序的具体修改方案：

方案1：修改pseud0子程序，支持自定义/自动种子 #

这个是简单的线性同余生成器，我们需要让它的初始种子可以外部传入，或者用系统时间自动初始化：

1. 先修改子程序，加入可选的种子参数，并显式保存状态变量（避免每次调用重置）：

   subroutine pseud0(r, seed)
   c generation of pseudo-random numbers
   real :: r
   integer, optional :: seed
   integer :: ir
   save ir  ! 必须保存ir的状态，不然每次调用都会回到初始值
   
   ! 第一次调用时初始化种子
   if (present(seed)) then
       ir = seed
   else
       ! 可选：如果没传种子，这里可以用系统时间自动生成，或者留一个默认值（但默认值别固定死）
       ! 临时测试的话可以让用户手动传不同的整数当种子
   endif
   
   ir = ir * 48828125
   if(ir) 10,20,20
   

10 ir = (ir + 2147483647) + 1
20 r = float(ir) * 0.4656613e-9
return
end

2. 在主程序里，第一次调用`pseud0`前，传入一个随时间变化的种子（比如系统时钟计数）：
```fortran
program main
integer :: seed
real :: rand_num

! 获取系统时钟计数当种子，保证每次运行都不同
call system_clock(count=seed)
! 初始化随机数生成器
call pseud0(rand_num, seed)
! 后续正常调用生成随机数，不用再传种子
call pseud0(rand_num)
! ... 你的模拟代码 ...
end program main

方案2：修改pseudo（Wichmann-Hill）子程序 #

这个是经典的组合生成器，同样需要把硬编码的种子改成可配置的：

1. 修改子程序，支持传入三个种子参数，并保存状态：

   subroutine pseudo(random, seed1, seed2, seed3)
   c wichmann+hill (1982) Appl. Statist 31
   real :: random
   integer, optional :: seed1, seed2, seed3
   integer :: ix, iy, iz
   save ix, iy, iz  ! 保存三个状态变量，关键！
   
   ! 初始化种子
   if (present(seed1) .and. present(seed2) .and. present(seed3)) then
       ix = seed1
       iy = seed2
       iz = seed3
   else
       ! 可选：自动生成种子，或者留默认（但别固定死）
   endif
   
   ix = 171 * mod(ix,177) - 2 * (ix/177)
   iy = 172 * mod(iy,176) - 35 * (iy/176)
   iz = 170 * mod(iz,178) - 63 * (iz/178)
   if (ix.lt.0) ix = ix + 30269
   if (iy.lt.0) iy = iy + 30307
   if (iz.lt.0) iz = iz + 30323
   random = mod(float(ix)/30269.0 + float(iy)/30307.0 + &
             float(iz)/30323.0, 1.0)
   return
   end
   

2. 主程序里用系统时钟生成三个不同的种子：

   program main
   integer :: seed, s1, s2, s3
   real :: rand_num
   
   call system_clock(count=seed)
   ! 把种子拆成三个符合每个分量范围的数（30269、30307、30323是各自的模数）
   s1 = mod(seed, 30269)
   s2 = mod(seed * 2, 30307)
   s3 = mod(seed * 3, 30323)
   ! 初始化生成器
   call pseudo(rand_num, s1, s2, s3)
   ! 后续正常调用
   call pseudo(rand_num)
   ! ... 你的模拟代码 ...
   end program main
   

额外提醒 #

• 一定要加save语句：旧Fortran编译器可能会默认保留data初始化的变量，但显式加save能确保状态在子程序调用之间不会被重置，这是保证随机数序列连续的关键。
• 如果用的是非常老的Fortran 77编译器，system_clock可能不支持，这时可以用编译器特定的函数（比如Intel Fortran的etime，GNU Fortran的date_and_time）来获取时间当种子，或者临时手动每次修改种子值测试。

内容的提问来源于stack exchange，提问作者Roland