c   This program performs calculations
c   for a proposed ASTM duration of load standard.
c
c   Steve Verrill, 5/6/97
c   Steve Verrill, 10/22/97 --- modified to reflect
c   my new understanding of the test protocol
c   Modified 9/1/00 to reflect latest test protocol.
c
c
c


      real red(5),cv(5),probf(5,5)
      real lnfrac
      real sigma(5)

      real samp1(200),samp2(200),samp3(200)
      real sort1(200),sort2(200),sort3(200)

      real strrat2(200),strrat3(200)

      real av(5),bv(5)

      integer iorigord(200)

      integer numf(5,5)
      
      character*80 results

      character*80 argv,intwrite

      call getarg(15,results)

      open(unit=7,file=results,access="append",status="unknown")

      call getarg(16,results)

      write(6,9900) results
 9900 format("The results from this run have been appended to",
     x       " the file ",a80,"
To see how to download", x " this file,", x ' read this description of anonymous ftp', x '.

','Alternatively, you should be able to do a', x ' "Save As" or "Print" in your browser.


') write(7,9700) 9700 format(////,1x,"**************************************",//,1x, x "Output from astm900.f",//,1x, x "**************************************",/) c c reductions in properties c call getarg(1,argv) write(intwrite,7) argv 7 format(a80) read(intwrite,*) numred call getarg(2,argv) write(intwrite,7) argv read(intwrite,*) (red(i),i = 1,numred) c av(j) and bv(j) must be calculated from the fact that the lines go c through point (ln(2160 = 24x90 hrs),1.0 - red(j)) and point c (ln(.05830 hr --- 3.5 minutes), 1.0) testlen = log(2160.0) tdiff = log(2160.0) - log(.0583) xbegin = log(.0583) do 8 i = 1,numred red(i) = red(i)/100.0 bv(i) = -red(i)/tdiff av(i) = 1.0 - bv(i)*xbegin 8 continue c c cvs c call getarg(3,argv) write(intwrite,7) argv read(intwrite,*) numcv call getarg(4,argv) write(intwrite,7) argv read(intwrite,*) (cv(i),i = 1,numcv) do 55 i = 1,numcv cv(i) = abs(cv(i)/100.0) c if ((cv(i) .eq. 0.0) .or. (cv(i) .ge. 1.0)) then c write(6,54) c 54 format("
A coefficient of variation value", c x " did not lie between 0 and 100. Please correct", c x " the form.
") c stop c endif 55 continue c c sample sizes c call getarg(5,argv) write(intwrite,7) argv read(intwrite,*) num1 write(6,61) num1 write(7,62) num1 61 format('
There were ',i3,' unaged observations', x '.
') 62 format(/,1x,'There were',i3,' unaged observations', x '.',/) call getarg(6,argv) write(intwrite,7) argv read(intwrite,*) num2 write(6,63) num2 write(7,64) num2 63 format('
There were ',i3,' batch2 observations', x '.
') 64 format(/,1x,'There were',i3,' batch2 observations', x '.',/) call getarg(7,argv) write(intwrite,7) argv read(intwrite,*) idos write(6,65) idos write(7,66) idos 65 format('
The ID of the order statistic', x' from the unaged sample that was used to', x' determine the base load was ',i3,'.
') 66 format(/,1x,'The ID of the order statistic', x' from the unaged sample that was used',/,1x, x'to determine the base load was ',i3,'.',/) call getarg(8,argv) write(intwrite,7) argv read(intwrite,*) frac write(6,67) frac write(7,68) frac 67 format('
The fraction of the base load', x' that was applied', x' to batches two and (possibly) three was ',f4.2,'.
') 68 format(/,1x,'The fraction of the base load', x' that was applied', x' to batches two and (possibly) three was ',f4.2,'.',/) call getarg(9,argv) write(intwrite,7) argv read(intwrite,*) numbok2 write(6,9967) numbok2 write(7,9968) numbok2 9967 format('
The number of failures', x' permitted', x' in the second batch was ',i3,'.
') 9968 format(/,1x,'The number of failures', x' permitted in the second', x' batch was ',i3,'.',/) call getarg(10,argv) write(intwrite,7) argv read(intwrite,*) num3 write(6,9963) num3 write(7,9964) num3 9963 format('
There were (possibly) ',i3,' batch3 observations', x '.
') 9964 format(/,1x,'There were (possibly) ',i3,' batch3 observations', x '.',/) call getarg(11,argv) write(intwrite,7) argv read(intwrite,*) numbok3 write(6,8867) numbok3 write(7,8868) numbok3 8867 format('
The number of failures', x' permitted', x' in the third batch was ',i3,'.
') 8868 format(/,1x,'The number of failures', x' permitted in the third', x' batch was ',i3,'.',/) call getarg(12,argv) write(intwrite,7) argv read(intwrite,*) rho rhocon = sqrt(1.0 - rho*rho) write(6,7767) rho write(7,7768) rho 7767 format('
The correlation in strength between', x' matched specimens was ',f4.2,'.
') 7768 format(/,1x,'The correlation in strength between', x' matched specimens was ',f4.2,'.',/) call getarg(13,argv) write(intwrite,7) argv read(intwrite,*) numtrial write(6,69) numtrial write(7,70) numtrial 69 format('
The number of simulation', x' trials requested was ',i6,'.
') 70 format(/,1x,'The number of simulation', x' trials requested was ',i6,'.',/) call getarg(14,argv) write(intwrite,7) argv read(intwrite,*) istart write(6,71) istart write(7,72) istart 71 format('
The start value for the random', x' number generator was ',i8,'.
') 72 format(/,1x,'The start value for the random', x' number generator was ',i8,'.',/) c c perform the simulation c do 110 i = 1,5 do 100 j = 1,5 numf(i,j) = 0 probf(i,j) = 0.0 100 continue 110 continue dum = rnor(istart) do 1000 iii = 1,numtrial c c generate the first sample c do 120 j = 1,num1 samp1(j) = rnor(0) 120 continue call rsort(samp1,sort1,iorigord,num1) c c generate the second sample c c Note that num2 is set equal to num1 in the Perl program c that calls this FORTRAN program. do 130 j = 1,num2 samp2(j) = rho*samp1(j) + rhocon*rnor(0) 130 continue call rsort(samp2,sort2,iorigord,num2) c c Draw the third batch. c Here I assume that the third batch is loaded c to the same level as the second batch. This c is not entirely clear from the standard. c c Problem: If there are 25, say, specimens in the c third batch (to bring 28 up to 53), then they c can't be matched 1 to 1 with the 28 specimens used c to determine the load. ... c Assuming that the extra 25 are matched c with 25 specimens other than the first 28, c then you could determine the load from testing c the 25 matched specimens to failure. Does one c use the 1st order statistic from the 25? But c this corresponds to a different order statistic c than the first out of 28 or the second out c of 53 .... ???? c If do all 53 at the start use the 2nd order statistic c from the 53 to determine the load and then only test 28 c for 90 days ... ???? Which 28? c do 132 j = 1,num3 samp3(j) = rnor(0) 132 continue call rsort(samp3,sort3,iorigord,num3) do 900 i = 1,numcv fac = 1.0 + cv(i)*sort1(idos) if (fac .lt. 0.0) then xload = 0.0 else xload = frac*fac endif do 140 j = 1,num2 den = 1.0 + cv(i)*sort2(j) if (xload .eq. 0.0) then strrat2(j) = 0.0 elseif (den .le. 0.0) then strrat2(j) = 1.0 else strrat2(j) = xload/den endif 140 continue do 142 j = 1,num3 den = 1.0 + cv(i)*sort3(j) if (xload .eq. 0.0) then strrat3(j) = 0.0 elseif (den .le. 0.0) then strrat3(j) = 1.0 else strrat3(j) = xload/den endif 142 continue do 800 j = 1,numred a = av(j) b = bv(j) numb = 0 do 700 k = 1,num2 xmu = (strrat2(k) - av(j))/bv(j) c c The standard deviation of the ln(time to failure) distribution c must be based on data --- variability of lifetimes c given a fixed stress ratio. From the plot of the Wood data on c page 3 of Chuck Gerhards' "Effect of Duration and Rate of Loading c on Strength of Wood and Wood-based Materials", I obtained an estimate c of sigma (from ranges) as 10.4 mm where a log base 10 (hours) unit is 15.2 mm. c Thus the sd in log base 10 (hours) units was 10.4/15.2 = .684 which leads c to a sd of 1.57 in log base e (hours) units. c xlnt = xmu + 1.57*rnor(0) if (xlnt .lt. testlen) then numb = numb + 1 endif 700 continue if (numb .gt. numbok2+1) then numf(i,j) = numf(i,j) + 1 elseif (numb .eq. numbok2+1 .and. num3 .eq. 0) then numf(i,j) = numf(i,j) + 1 elseif (numb .eq. numbok2+1 .and. num3 .gt. 0) then numb = 0 do 710 k = 1,num3 xmu = (strrat3(k) - av(j))/bv(j) c c The standard deviation of the ln(time to failure) distribution c must be based on data --- variability of lifetimes c given a fixed stress ratio. From the plot of the Wood data on c page 3 of Chuck Gerhards' "Effect of Duration and Rate of Loading c on Strength of Wood and Wood-based Materials", I obtained an estimate c of sigma (from ranges) as 10.4 mm where a log base 10 (hours) unit is 15.2 mm. c Thus the sd in log base 10 (hours) units was 10.4/15.2 = .684 which leads c to a sd of 1.57 in log base e (hours) units. c xlnt = xmu + 1.57*rnor(0) if (xlnt .lt. testlen) then numb = numb + 1 endif 710 continue if (numb .gt. numbok3) x numf(i,j) = numf(i,j) + 1 endif 800 continue 900 continue 1000 continue xnumt = numtrial do 1007 i = 1,numcv do 1005 j = 1,numred probf(i,j) = numf(i,j)/xnumt 1005 continue 1007 continue write(6,1010) write(7,1012) (red(i), i = 1,numred) 1010 format('



Normal Strength', x' Distribution

', x'Reductions are on top. Coefficients', x' of variation are along the left side. Probabilities', x' of failing the test are in the table.

') 1012 format(///,1x,'Normal Strength Distribution',//, x1x,'Reductions are on top. Coefficients of', x' variation',/,1x,'are along the left side. ', x'Probabilities of failing the test',/,1x,'are in the table.', x//, x5x,5(3x,f4.2),/) do 1016 i = 1,numcv write(7,1014) cv(i),(probf(i,j),j = 1,numred) 1014 format(1x,f4.2,5(2x,f5.3)) 1016 continue write(6,1020) 1020 format("") write(6,1030) numred 1030 format("", x"") do 1034 i = 1,numred write(6,1032) red(i) 1032 format("") 1034 continue write(6,1036) 1036 format("") do 1050 i = 1,numcv write(6,1042) cv(i),(probf(i,j), j = 1,numred) 1042 format("",("")) write(6,1044) 1044 format("") 1050 continue write(6,1052) 1052 format("
Coefficient
of
variation
Reduction
",f4.2,"
",f5.3,"
") c c The lognormal case c lnfrac = log(frac) do 3020 i = 1,numcv sigma(i) = sqrt(log(cv(i)*cv(i) + 1.0)) 3020 continue c c perform the simulation c do 3110 i = 1,5 do 3100 j = 1,5 numf(i,j) = 0 probf(i,j) = 0.0 3100 continue 3110 continue do 4000 iii = 1,numtrial c c generate the first sample c do 3120 j = 1,num1 samp1(j) = rnor(0) 3120 continue call rsort(samp1,sort1,iorigord,num1) c c generate the second sample c do 3130 j = 1,num2 samp2(j) = rho*samp1(j) + rhocon*rnor(0) 3130 continue call rsort(samp2,sort2,iorigord,num2) c c Draw the third batch. c Here I assume that the third batch is loaded c to the same level as the second batch. This c is not entirely clear from the standard. c c Problem: If there are 25, say, specimens in the c third batch (to bring 28 up to 53), then they c can't be matched 1 to 1 with the 28 specimens used c to determine the load. ... c Assuming that the extra 25 are matched c with 25 specimens other than the first 28, c then you could determine the load from testing c the 25 matched specimens to failure. Does one c use the 1st order statistic from the 25? But c this corresponds to a different order statistic c than the first out of 28 or the second out c of 53 .... ???? c If do all 53 at the start use the 2nd order statistic c from the 53 to determine the load and then only test 28 c for 90 days ... ???? Which 28? c do 3132 j = 1,num3 samp3(j) = rnor(0) 3132 continue call rsort(samp3,sort3,iorigord,num3) do 3900 i = 1,numcv xload = exp(lnfrac + sigma(i)*sort1(idos)) do 3140 j = 1,num2 strrat2(j) = xload/exp(sigma(i)*sort2(j)) 3140 continue do 3142 j = 1,num3 strrat3(j) = xload/exp(sigma(i)*sort3(j)) 3142 continue do 3800 j = 1,numred a = av(j) b = bv(j) numb = 0 do 3700 k = 1,num2 xmu = (strrat2(k) - av(j))/bv(j) c c The standard deviation of the ln(time to failure) distribution c must be based on data --- variability of lifetimes c given a fixed stress ratio. From the plot of the Wood data on c page 3 of Chuck Gerhards' "Effect of Duration and Rate of Loading c on Strength of Wood and Wood-based Materials", I obtained an estimate c of sigma (from ranges) as 10.4 mm where a log base 10 (hours) unit is 15.2 mm. c Thus the sd in log base 10 (hours) units was 10.4/15.2 = .684 which leads c to a sd of 1.57 in log base e (hours) units. c xlnt = xmu + 1.57*rnor(0) if (xlnt .lt. testlen) then numb = numb + 1 endif 3700 continue if (numb .gt. numbok2+1) then numf(i,j) = numf(i,j) + 1 elseif (numb .eq. numbok2+1 .and. num3 .eq. 0) then numf(i,j) = numf(i,j) + 1 elseif (numb .eq. numbok2+1 .and. num3 .gt. 0) then numb = 0 do 3710 k = 1,num3 xmu = (strrat3(k) - av(j))/bv(j) c c The standard deviation of the ln(time to failure) distribution c must be based on data --- variability of lifetimes c given a fixed stress ratio. From the plot of the Wood data on c page 3 of Chuck Gerhards' "Effect of Duration and Rate of Loading c on Strength of Wood and Wood-based Materials", I obtained an estimate c of sigma (from ranges) as 10.4 mm where a log base 10 (hours) unit is 15.2 mm. c Thus the sd in log base 10 (hours) units was 10.4/15.2 = .684 which leads c to a sd of 1.57 in log base e (hours) units. c xlnt = xmu + 1.57*rnor(0) if (xlnt .lt. testlen) then numb = numb + 1 endif 3710 continue if (numb .gt. numbok3) x numf(i,j) = numf(i,j) + 1 endif 3800 continue 3900 continue 4000 continue xnumt = numtrial do 4007 i = 1,numcv do 4005 j = 1,numred probf(i,j) = numf(i,j)/xnumt 4005 continue 4007 continue write(6,4010) write(7,4012) (red(i), i = 1,numred) 4010 format('



Lognormal Strength', x' Distribution

', x'Reductions are on top. Coefficients', x' of variation are along the left side. Probabilities', x' of failing the test are in the table.

') 4012 format(///,1x,'Lognormal Strength Distribution',//, x1x,'Reductions are on top. Coefficients of', x' variation',/,1x,'are along the left side. ', x'Probabilities of failing the test',/,1x,'are in the table.', x//, x5x,5(3x,f4.2),/) do 4016 i = 1,numcv write(7,4014) cv(i),(probf(i,j),j = 1,numred) 4014 format(1x,f4.2,5(2x,f5.3)) 4016 continue write(6,4020) 4020 format("") write(6,4030) numred 4030 format("", x"") do 4034 i = 1,numred write(6,4032) red(i) 4032 format("") 4034 continue write(6,4036) 4036 format("") do 4050 i = 1,numcv write(6,4042) cv(i),(probf(i,j), j = 1,numred) 4042 format("",("")) write(6,4044) 4044 format("") 4050 continue write(6,4052) 4052 format("
Coefficient
of
variation
Reduction
",f4.2,"
",f5.3,"
") write(6,6060) 6060 format('



For further information,', x' please contact Steve Verrill at ', x'', x'sverrill@fs.fed.us', x' or 608-231-9375.
Last modified on 9/7/00.') stop end REAL FUNCTION RNOR(JD) C***BEGIN PROLOGUE RNOR C***DATE WRITTEN 810915 C***REVISION DATE 830805 C***CATEGORY NO. L6A14 C***KEYWORDS RANDOM NUMBERS, UNIFORM RANDOM NUMBERS C***AUTHOR KAHANER, DAVID, SCIENTIFIC COMPUTING DIVISION, NBS C MARSAGLIA, GEORGE, COMPUTER SCIENCE DEPT., WASH STATE UNIV C C***PURPOSE GENERATES QUASI NORMAL RANDOM NUMBERS, WITH MEAN ZERO AND C UNIT STANDARD DEVIATION, AND CAN BE USED WITH ANY COMPUTER C WITH INTEGERS AT LEAST AS LARGE AS 32767. C***DESCRIPTION C C RNOR generates quasi normal random numbers with zero mean and C unit standard deviation. C It can be used with any computer with integers at least as C large as 32767. C C C Use C First time.... C Z = RNOR(JD) C Here JD is any n o n - z e r o integer. C This causes initialization of the program C and the first random number to be returned as Z. C Subsequent times... C Z = RNOR(0) C Causes the next random number to be returned as Z. C C..................................................................... C C Note: Users who wish to transport this program to other C computers should read the following .... C C Machine dependencies... C MDIG = A lower bound on the number of binary digits available C for representing integers, including the sign bit. C This must be at least 16, but can be increased in C line with remark A below. C C Remarks... C A. This program can be used in two ways: C (1) To obtain repeatable results on different computers, C set 'MDIG' to the smallest of its values on each, or, C (2) To allow the longest sequence of random numbers to be C generated without cycling (repeating) set 'MDIG' to the C largest possible value. C B. The sequence of numbers generated depends on the initial C input 'JD' as well as the value of 'MDIG'. C If MDIG=16 one should find that C the first evaluation C Z=RNOR(87) gives Z=-.40079207... C The second evaluation C Z=RNOR(0) gives Z=-1.8728870... C The third evaluation C Z=RNOR(0) gives Z=1.8216004... C The fourth evaluation C Z=RNOR(0) gives Z=.69410355... C The thousandth evaluation C Z=RNOR(0) gives Z=.96782424... C C***REFERENCES MARSAGLIA & TSANG, "A FAST, EASILY IMPLEMENTED C METHOD FOR SAMPLING FROM DECREASING OR C SYMMETRIC UNIMODAL DENSITY FUNCTIONS", TO BE C PUBLISHED IN SIAM J SISC 1983. C***ROUTINES CALLED I1MACH,XERROR C***END PROLOGUE RNOR REAL V(65),W(65) INTEGER M(17) SAVE I1,J1,M,M1,M2,RMAX DATA AA,B,C,RMAX/12.37586,.4878992,12.67706,3.0518509E-5/ DATA C1,C2,PC,XN/.9689279,1.301198,.1958303E-1,2.776994/ DATA V/ .3409450, .4573146, .5397793, .6062427, .6631691 +, .7136975, .7596125, .8020356, .8417227, .8792102, .9148948 +, .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.4834526, 1.5100121, 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.9155830, 1.9492166 +,1.9839239, 2.0198430, 2.0571356, 2.0959930, 2.1366450, 2.1793713 +,2.2245175, 2.2725185, 2.3239338, 2.3795007, 2.4402218, 2.5075117 +,2.5834658, 2.6713916, 2.7769943, 2.7769943, 2.7769943, 2.7769943/ DATA W/ .10405134E-04, .13956560E-04, .16473259E-04, + .18501623E-04, .20238931E-04, .21780983E-04, .23182241E-04, + .24476931E-04, .25688121E-04, .26832186E-04, .27921226E-04, + .28964480E-04, .29969191E-04, .30941168E-04, .31885160E-04, + .32805121E-04, .33704388E-04, .34585827E-04, .35451919E-04, + .36304851E-04, .37146564E-04, .37978808E-04, .38803170E-04, + .39621114E-04, .40433997E-04, .41243096E-04, .42049621E-04, + .42854734E-04, .43659562E-04, .44465208E-04, .45272764E-04, + .46083321E-04, .46897980E-04, .47717864E-04, .48544128E-04, + .49377973E-04, .50220656E-04, .51073504E-04, .51937936E-04, + .52815471E-04, .53707761E-04, .54616606E-04, .55543990E-04, + .56492112E-04, .57463436E-04, .58460740E-04, .59487185E-04, + .60546402E-04, .61642600E-04, .62780711E-04, .63966581E-04, + .65207221E-04, .66511165E-04, .67888959E-04, .69353880E-04, + .70922996E-04, .72618816E-04, .74471933E-04, .76525519E-04, + .78843526E-04, .81526890E-04, .84749727E-04, + .84749727E-04, .84749727E-04, .84749727E-04/ DATA M(1),M(2),M(3),M(4),M(5),M(6),M(7),M(8),M(9),M(10),M(11), 1 M(12),M(13),M(14),M(15),M(16),M(17) 2 / 30788,23052,2053,19346,10646,19427,23975, 3 19049,10949,19693,29746,26748,2796,23890, 4 29168,31924,16499 / DATA M1,M2,I1,J1 / 32767,256,5,17 / C Fast part... C C C***FIRST EXECUTABLE STATEMENT RNOR IF(JD.NE.0)GO TO 27 10 CONTINUE I=M(I1)-M(J1) IF(I .LT. 0) I=I+M1 M(J1)=I I1=I1-1 IF(I1 .EQ. 0) I1=17 J1=J1-1 IF(J1 .EQ. 0) J1=17 J=MOD(I,64)+1 RNOR=I*W(J+1) IF( ( (I/M2)/2 )*2.EQ.(I/M2))RNOR=-RNOR IF(ABS(RNOR).LE.V(J))RETURN C Slow part; AA is a*f(0) X=(ABS(RNOR)-V(J))/(V(J+1)-V(J)) Y=UNI(0) 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*RNOR**2))RETURN C Tail part; 3.855849 is .5*XN**2 c c the three lines after the next empty comment line were added by c Steve Verrill on 10/23/89 c 22 tttttt = uni(0) if (tttttt .eq. 0.0) go to 22 s = xn - alog(tttttt)/xn c 22 S=XN-ALOG(UNI(0))/XN c c the three lines after the next empty comment line were c added by Steve Verrill on 10/23/89 c 23 tttttt = uni(0) if (tttttt .eq. 0.0) go to 23 if (3.855849 + alog(tttttt) - xn*s .gt. -.5*s**2) go to 22 c IF(3.855849+ALOG(UNI(0))-XN*S.GT.-.5*S**2)GO TO 22 RNOR=SIGN(S,RNOR) RETURN 11 RNOR=SIGN(B-B*X,RNOR) RETURN C FILL 27 CONTINUE MDIG=32 C BE SURE THAT MDIG AT LEAST 16... ccc IF(MDIG.LT.16)CALL XERROR('RNOR--MDIG LESS THAN 16',23,1,2) M1 = 2**(MDIG-2) + (2**(MDIG-2)-1) M2 = 2**(MDIG/2) JSEED = MIN0(IABS(JD),M1) IF( MOD(JSEED,2).EQ.0 ) JSEED=JSEED-1 K0 =MOD(9069,M2) K1 = 9069/M2 J0 = MOD(JSEED,M2) J1 = JSEED/M2 DO 2 I=1,17 JSEED = J0*K0 J1 = MOD(JSEED/M2+J0*K1+J1*K0,M2/2) J0 = MOD(JSEED,M2) 2 M(I) = J0+M2*J1 J1=17 I1=5 RMAX = 1./FLOAT(M1) C Seed uniform (0,1) generator. (Just a dummy call) RNOR=UNI(JD) DO 28 I=1,65 28 W(I)=RMAX*V(I) GO TO 10 END ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc REAL FUNCTION UNI(JD) C***BEGIN PROLOGUE UNI C***DATE WRITTEN 810915 C***REVISION DATE 830805 C***CATEGORY NO. L6A21 C***KEYWORDS RANDOM NUMBERS, UNIFORM RANDOM NUMBERS C***AUTHOR BLUE, JAMES, SCIENTIFIC COMPUTING DIVISION, NBS C KAHANER, DAVID, SCIENTIFIC COMPUTING DIVISION, NBS C MARSAGLIA, GEORGE, COMPUTER SCIENCE DEPT., WASH STATE UNIV C C***PURPOSE THIS ROUTINE GENERATES QUASI UNIFORM RANDOM NUMBERS ON [0,1) C AND CAN BE USED ON ANY COMPUTER WITH WHICH ALLOWS INTEGERS C AT LEAST AS LARGE AS 32767. C***DESCRIPTION C C THIS ROUTINE GENERATES QUASI UNIFORM RANDOM NUMBERS ON THE INTERVAL C [0,1). IT CAN BE USED WITH ANY COMPUTER WHICH ALLOWS C INTEGERS AT LEAST AS LARGE AS 32767. C C C USE C FIRST TIME.... C Z = UNI(JD) C HERE JD IS ANY N O N - Z E R O INTEGER. C THIS CAUSES INITIALIZATION OF THE PROGRAM C AND THE FIRST RANDOM NUMBER TO BE RETURNED AS Z. C SUBSEQUENT TIMES... C Z = UNI(0) C CAUSES THE NEXT RANDOM NUMBER TO BE RETURNED AS Z. C C C.................................................................. C NOTE: USERS WHO WISH TO TRANSPORT THIS PROGRAM FROM ONE COMPUTER C TO ANOTHER SHOULD READ THE FOLLOWING INFORMATION..... C C MACHINE DEPENDENCIES... C MDIG = A LOWER BOUND ON THE NUMBER OF BINARY DIGITS AVAILABLE C FOR REPRESENTING INTEGERS, INCLUDING THE SIGN BIT. C THIS VALUE MUST BE AT LEAST 16, BUT MAY BE INCREASED C IN LINE WITH REMARK A BELOW. C C REMARKS... C A. THIS PROGRAM CAN BE USED IN TWO WAYS: C (1) TO OBTAIN REPEATABLE RESULTS ON DIFFERENT COMPUTERS, C SET 'MDIG' TO THE SMALLEST OF ITS VALUES ON EACH, OR, C (2) TO ALLOW THE LONGEST SEQUENCE OF RANDOM NUMBERS TO BE C GENERATED WITHOUT CYCLING (REPEATING) SET 'MDIG' TO THE C LARGEST POSSIBLE VALUE. C B. THE SEQUENCE OF NUMBERS GENERATED DEPENDS ON THE INITIAL C INPUT 'JD' AS WELL AS THE VALUE OF 'MDIG'. C IF MDIG=16 ONE SHOULD FIND THAT C THE FIRST EVALUATION C Z=UNI(305) GIVES Z=.027832881... C THE SECOND EVALUATION C Z=UNI(0) GIVES Z=.56102176... C THE THIRD EVALUATION C Z=UNI(0) GIVES Z=.41456343... C THE THOUSANDTH EVALUATION C Z=UNI(0) GIVES Z=.19797357... C C***REFERENCES MARSAGLIA G., "COMMENTS ON THE PERFECT UNIFORM RANDOM C NUMBER GENERATOR", UNPUBLISHED NOTES, WASH S. U. C***ROUTINES CALLED I1MACH,XERROR C***END PROLOGUE UNI INTEGER M(17) C SAVE I,J,M,M1,M2 C DATA M(1),M(2),M(3),M(4),M(5),M(6),M(7),M(8),M(9),M(10),M(11), 1 M(12),M(13),M(14),M(15),M(16),M(17) 2 / 30788,23052,2053,19346,10646,19427,23975, 3 19049,10949,19693,29746,26748,2796,23890, 4 29168,31924,16499 / DATA M1,M2,I,J / 32767,256,5,17 / C***FIRST EXECUTABLE STATEMENT UNI IF(JD .EQ. 0) GO TO 3 C FILL MDIG=32 C BE SURE THAT MDIG AT LEAST 16... c IF(MDIG.LT.16)CALL XERROR('UNI--MDIG LESS THAN 16',22,1,2) M1= 2**(MDIG-2) + (2**(MDIG-2)-1) M2 = 2**(MDIG/2) JSEED = MIN0(IABS(JD),M1) IF( MOD(JSEED,2).EQ.0 ) JSEED=JSEED-1 K0 =MOD(9069,M2) K1 = 9069/M2 J0 = MOD(JSEED,M2) J1 = JSEED/M2 DO 2 I=1,17 JSEED = J0*K0 J1 = MOD(JSEED/M2+J0*K1+J1*K0,M2/2) J0 = MOD(JSEED,M2) 2 M(I) = J0+M2*J1 I=5 J=17 C BEGIN MAIN LOOP HERE 3 K=M(I)-M(J) IF(K .LT. 0) K=K+M1 M(J)=K I=I-1 IF(I .EQ. 0) I=17 J=J-1 IF(J .EQ. 0) J=17 UNI=FLOAT(K)/FLOAT(M1) RETURN END c This subroutine attempts to perform an n*log(n) sort c rather than an n*n sort. c It sorts reals. c It also keeps track of the original order of the x's. c c c Steve Verrill, 5/2/89 c c c subroutine rsort(x,xsort,iorigord,n) c dimension x(n),xsort(n),iorigord(n) c do 10 i = 1,n iorigord(i) = i 10 continue c num = 2 c if (x(1) .lt. x(2)) then xsort(1) = x(1) iorigord(1) = 1 xsort(2) = x(2) iorigord(2) = 2 else xsort(1) = x(2) iorigord(1) = 2 xsort(2) = x(1) iorigord(2) = 1 endif c do 1000 j = 3,n c c insert the next value c if (x(j) .ge. xsort(num)) go to 200 if (x(j) .le. xsort(1)) go to 160 il = 1 iu = num 145 numb = iu - il if (numb .eq. 1) go to 170 numbd2 = numb/2 itest = il + numbd2 if (x(j) .le. xsort(itest)) go to 150 il = itest go to 145 150 iu = itest go to 145 160 il = 0 170 num = num + 1 do 180 k = num,il + 2,-1 km1 = k - 1 xsort(k) = xsort(km1) iorigord(k) = iorigord(km1) 180 continue 190 ilp1 = il + 1 xsort(ilp1) = x(j) iorigord(ilp1) = j go to 1000 c 200 num = num + 1 xsort(num) = x(j) iorigord(num) = j c 1000 continue c return end