c モンテカルロ法による円の面積の計算
      write(6,*) 'nstep '
      read(5,*) nstep
      write(6,*) 'iseed ?'
      read(5,*) iseed

      r0=1.0               /* 円の半径*/
      r02=r0**2

      nin=0
      do 10 i=1,nstep
        x=rand(iseed)
        y=rand(iseed)
        r2=x**2+y**2       /* 原点からの距離の2乗 */
        if(r2.le.r02) nin=nin+1   /* 円の中 */
 10   continue
      s=4.0 *(float(nin)/float(nstep))
      write(6,*) 'nstep= ',nstep,' s= ',s
      stop
      end 

      function irand(iseed)
      iseed=iseed*48828125
      if(iseed.le.0) iseed=(iseed+2147483647)+1
      irand=iseed
      return
      end

      function rand(iseed)
      iseed=iseed*48828125
      if(iseed.le.0) iseed=(iseed+2147483647)+1
      rand=float(iseed)/2147483647
      return
      end