real x(1000) character*15 filename,fileout c**** converts binary file of 32-bit integers to c**** single precision standard normal variables. print 12345 12345 format(10x,'This program will read a binary file of',/, &10x,'32-bit random integers and write a binary file',/, &10x,'of standard normal random variables, using the',/, &10x,'ziggurat method of Marsaglia and Tsang, ',/, &10x,'SIAM J. Scient. and Stat. Computing, v5, 349-359, 1984',/) print*,' Enter name of source file ( binary, 32-bit integers):' read 20,filename 20 format(a15) open(1,file=filename,form='unformatted',access='direct', & recl=16384) print*,' Enter name of output file:' read 20,fileout open(2,file=fileout,form='unformatted',access='direct', & recl=4000) print*,' Enter n for output of RNOR''S (in thousands):' read*,n print*,' Please wait ...............' jk=0 do 2 i=1,n do 3 k=1,1000 3 x(k)=zignor() jk=jk+1 2 write(2,rec=jk) x close(2) print*,' Finished' print 25, n,fileout 25 format(i10,'000 standard normal random variables',/, &' have been written to the file ',a15) go to 123 11111 print*,' ERROR: not enough integers in source file' 123 end function jtbl() integer b(4096) data j,jk/4097,1/ if(j.gt.4096) then read(1,rec=jk) b j=1 jk=jk+1 endif jtbl=b(j) j=j+1 return end c***Ziggurat rnor. Takes 2.2 microsecs(Microsoft, 5.9 in Lahey) c**************** takes 2.0 if jsr is only rng, 2.2 if jsr+icng c*** mpynor takes 3.1, rnor92 takes 3.8 function zignor() real v(0:64) parameter(aa=12.37586,b=.4878992,c=12.67706,rmax=.4656613e-9, + c1=.9689279,c2=1.301198,pc=.1958303e-1,xn=2.776994) data v/.3409450,.4573146,.5397793,.6062427,.6631691,.7136975, + .7596125, .8020356, .8417227, .8792102, .9148948, 0.9490791, + .9820005, 1.0138492, 1.0447810, 1.0749254, 1.1043917, 1.1332738, + 1.1616530, 1.1896010, 1.2171815, 1.2444516, 1.2714635, 1.2982650, + 1.3249008, 1.3514125, 1.3778399, 1.4042211, 1.4305929, 1.4569915, + 1.4834527, 1.5100122, 1.5367061, 1.5635712, 1.5906454, 1.6179680, + 1.6455802, 1.6735255, 1.7018503, 1.7306045, 1.7598422, 1.7896223, + 1.8200099, 1.8510770, 1.8829044, 1.9155831, 1.9492166, 1.9839239, + 2.0198431, 2.0571356, 2.0959930, 2.1366450, 2.1793713, 2.2245175, + 2.2725186, 2.3239338, 2.3795008, 2.4402218, 2.5075117, 2.5834658, + 2.6713916, 2.7769942, 2.7769942, 2.7769942, 2.7769942/ uni()=.5+jtbl()*.5**32 c--FAST PART i=jtbl() j=and(i,63) zignor=i*rmax*v(j+1) if(abs(zignor) .le. v(j)) return c-------SLOW PART; aa is a*f(0) y=uni() s=x+y if(s .gt. c2) go to 11 if(s .le. c1) return if(y .gt. c-aa*exp(-.5*(b-b*x)**2)) go to 11 if(exp(-.5*v(j+1)**2)+y*pc/v(j+1) .le. exp(-.5*zignor**2)) return c---------------TAIL PART; .3601016 is 1/xn 22 x=.3601016*alog(uni()) if(-2.*alog(uni()) .le. x*x) go to 22 33 zignor=sign(xn-x,zignor) return 11 zignor=sign(b-b*x,zignor) return end