一个关于dlog算不出来的问题 - Fortran Coder

[Fortran] 纯文本查看 复制代码

program main

    implicit real(kind=8)(a-h,o-z)

    

    pi=3.1416926d0

    vm_ini=2.0d0

    abc1=dsqrt(-dlog(rf(0))) 

    abc2=2.0d0*pi*rf(0)   

    p=abc1*dsin(abc2)*vm_ini

    end

    

    function rf(idum)           

    implicit real(kind=8)(a-h,o-z)

    save ma,inext,inextp

    parameter (mbig=1000000000,mseed=161803398,mz=0,fac=1.e-9)

    dimension ma(55)

    data iff/0/

    if(idum<0.or.iff==0) then

        iff=1

        mj=mseed-iabs(idum)

        mj=mod(mj,mbig)

        ma(55)=mj

        mk=1

        do 50 i=1,54

            ii=mod(21*i,55)

            ma(ii)=mk

            mk=mj-mk

            if(mk<mz) mk=mk+mbig

            mj=ma(ii)

50      continue

        do 100 k=1,4

            do 60 i=1,55

                ma(i)=ma(i)-ma(1+mod(i+30,55))

                if(ma(i)<mz) ma(i)=ma(i)+mbig

60          continue

100     continue

        inext=0

        inextp=31

    endif

200 inext=inext+1

    if(inext==56) inext=1

    inextp=inextp+1

    if(inextp==56) inextp=1

    mj=ma(inext)-ma(inextp)

    if(mj<mz) mj=mj+mbig

    ma(inext)=mj

    rf=mj*fac

    if(rf>1.e-8.and.rf<0.999999999) return

    goto 200

    pause

    end