module my_type INTEGER, PARAMETER :: IB = SELECTED_INT_KIND(9) !integer*4 INTEGER, PARAMETER :: SP = SELECTED_REAL_KIND(6,37) !real*4 INTEGER, PARAMETER :: DP = SELECTED_REAL_KIND(15,307) !real*8 end module my_type program circle use my_type implicit none integer(ib) :: seed real(sp) :: x1,y1,ran2 real(dp) :: vol,summ,var integer(ib) :: n,i,it integer(ib),dimension(4) :: nit data nit /100000,1000000,10000000,100000000/ do it=1,3 n=nit(it) seed=-5 vol=1. summ=0._sp do i=1,n x1=ran2(seed) y1=ran2(seed) if(x1**2+y1**2 .le. 1.) then summ=summ+0.1 var=var+0.1*0.1 endif enddo summ=summ*vol/n var=vol*sqrt(((var/n)-(summ/n)**2)/n) write(*,10) n,40.*summ,40.*var 10 format(' n =',i10,' pi = ',f10.6,' +/- ',e12.6) enddo end FUNCTION ran2(idum) !Langperiodischer (>2*10^18) Zufallszahlengenerator von L Ecuyer, ergibt !gleichförmig verteilte Zufallszahlen im Bereich (0,1). !Bzgl. Details, siehe NumRec, Kap. 7.1. !Beim ersten Aufruf ist idum mit einer negativen Zahl zu initialisieren !und danch nicht mehr zu ändern. Beispiel ! program main ! use my_type ! implicit none ! integer(ib) :: idum,i ! real(sp) :: x, ran2 ! data idum /-11/ ! do i=1,10000 ! x=ran2(idum) ! enddo ! end ! !RNMX soll die größte Real-Zahl < 1 approximieren, ist bei Verwendung von !double prec zu ändern ! use my_type implicit none INTEGER(ib) :: idum REAL (sp):: ran2 INTEGER(ib), parameter :: IM1=2147483563,IM2=2147483399,IMM1=IM1-1, & & IA1=40014,IA2=40692,IQ1=53668,IQ2=52774, & & IR1=12211,IR2=3791,NTAB=32,NDIV=1+IMM1/NTAB REAL(sp), PARAMETER :: AM=1./IM1,EPS=1.2e-7,RNMX=1.-EPS INTEGER(ib) :: idum2,j,k,iy INTEGER(ib), dimension(ntab) :: iv SAVE iv,iy,idum2 DATA idum2/123456789/, iv/NTAB*0/, iy/0/ if (idum.le.0) then idum=max(-idum,1) idum2=idum do j=NTAB+8,1,-1 k=idum/IQ1 idum=IA1*(idum-k*IQ1)-k*IR1 if (idum.lt.0) idum=idum+IM1 if (j.le.NTAB) iv(j)=idum enddo iy=iv(1) endif k=idum/IQ1 idum=IA1*(idum-k*IQ1)-k*IR1 if (idum.lt.0) idum=idum+IM1 k=idum2/IQ2 idum2=IA2*(idum2-k*IQ2)-k*IR2 if (idum2.lt.0) idum2=idum2+IM2 j=1+iy/NDIV iy=iv(j)-idum2 iv(j)=idum if(iy.lt.1)iy=iy+IMM1 ran2=min(AM*iy,RNMX) return END