Fortran生成正态分布随机数

分享一个fortran生成正态分布随机数的自函数，拿来即用，已经在intel的ifort编译器上测试通过

      subroutine gaussrandom(y1,y2)
c output y1 and y2 obey normal distribution of (0,1)
      real::x1,x2,y1,y2 
      data pii/3.14159265/
      call random_seed()
      call random_number(x1)
      call random_number(x2)
      y1=sqrt(-2*log(x1))*cos(2*pii*x2)
      y2=sqrt(-2*log(x1))*sin(2*pii*x2)
      return 
      end 

将Y1和Y2输出，并测试是否为正态分布。测试结果如下：