/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:19 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "ranpk2.h" #include #include /* PARAMETER translations */ #define X1DP 123456789.0e0 #define X1SP 123456789.0e0 /* end of PARAMETER translations */ /* COMMON translations */ struct t_rancom { double xcurdp; float xcursp; long int mode; } rancom; /* end of COMMON translations */ void /*FUNCTION*/ rn1() { /* Copyright (c) 1996 California Institute of Technology, Pasadena, CA. * ALL RIGHTS RESERVED. * Based on Government Sponsored Research NAS7-03001. *>> 2006-04-16 RANPK2 Krogh Fixed invalid mode message. *>> 2005-12-21 RANPK2 Krogh Changec ivalid to invalid is stops *>> 2005-12-07 RANPK2 Krogh Fixed bug in C conversion. *>> 1997-12-17 RANPK2 Krogh Removed unreferenced labels *>> 1995-11-22 RANPK2 Krogh Removed multiple entries. *>> 1992-03-17 CLL Moved SAVE stmt ahead of DATA stmts. *>> 1992-03-09 CLL Removed "save FIRST" because FIRST is not used. *>> 1992-03-02 CLL Fix error: Set MODE = 1 in data stmt. *>> 1991-11-21 CLL Add MODE with values 1, 2, 3, & 4. *>> 1989-09-11 CLL Multiversion file. RANPK2 or RANPK3 *>> 1987-05-05 RANPK2 Lawson Initial code. * * This program unit, along with RANPK1, supports random number * generation. * * The functionality of this random number package was modeled on the * specifications of the random number routines RANDOM and RANDOMSEED * in the February, 1987 working draft of the Fortran 8x language * standard. This functionality remains similar to the later draft, * Fortran 90, S8.115, June 1990, in which the routine names have * been changed to RANDOM_MUMBER and RANDOM_SEED. This should * facilitate replacement of use of this package by Fortran intrinsic * when and if Fortran 90 compilers come into widespread use. * * The library user may call RANSIZ or RANGET in this prog unit, * or RAN1 or RANPUT in RANPK1 to obtain the functionality of * RANDOM_SEED of Fortran Fixed bug in C conversion.90. This relates to sett * the seed. * Entries RN1 and RNPUT in this prog unit should not be called by * library users. They are intended only to be called from RANPK1. * Entry RN2 returns the value of MODE. This is a convenience in * case one is interested in knowing the value of MODE the package * has selected. * Entries SRANUA and DRANUA (s.p. and d.p. respectively) * generate arrays of pseudorandom numbers from the uniform * distribution in the range [0., 1.]. These may be called by users * and are called by other library routines. * Entries SRANUS and DRANUA (s.p. and d.p. respectively) * generate arrays of pseudorandom numbers from the uniform * distribution in the range [0., 1.] and then transformed to * A + B*U. These are intended for direct use by users. * ------------------------------------------------------------------ * Algorithm for generating pseudorandom numbers * from the uniform distribution. * * The current integer value in the random integer sequence is XCUR, * and the next is defined mathematically by the statement: * * XCUR = mod(AFAC * XCUR, MDIV) * where * MDIV = m = 6_87194_76503 = 2**36 - 233 * and * AFAC = a = 612_662 = (approx.) 0.58 * 2**20 * * XCUR may be any integer value in the range from 1 to m-1, and all * integer values in this range will be generated before the sequence * cycles. * * We call the above computational algorithm for XCUR the "short" * algorithm. There is also a "long" algorithm that produces exactly * the same sequence of values of XCUR: * Q = aint(XCUR/B) * R = XCUR - Q * B * XCUR = AFAC * R - C * Q * do while(XCUR .lt. 0.0) * XCUR = XCUR + MDIV * end do * where B and C are constants related to MDIV and AFAC by * MDIV = B * AFAC + C * We use B = 112165 and C = 243273. The average number of executions * of the statement XCUR = XCUR + MDIV is 1.09 and the maximum number of * executions is 3. * * The largest number that must be handled in the "short" algorithm * is the product of AFAC with the max value of XCUR, i.e., * 612_662 * 6_87194_76502 = 42_10181_19126_68324 ~= 0.58 * 2**56. * Thus this algorithm requires arithmetic exact to at least 56 bits. * * The largest number that must be handled in the "long" algorithm * is the product of C with the max value of aint(XCUR/B), i.e., * 243273 * 612664 ~= 0.14904e12 ~= 0.54 * 2**38. * Thus this algorithm requires arithmetic exact to at least 38 bits. * * To accommodate different compiler/computer systems this program unit * contains code for 3 different ways of computing the new XCUR from the * old XCUR, each producing exactly the same sequence of of XCUR. * * Initially we have MODE = 1. When MODE = 1 the code does tests to * see which of three implementation methods will be used, and sets * MODE = 2, 3, or 4 to indicate the choice. * * Mode 2 will be used in machines such as the Cray that have at * least a 38 bit significand in SP arithmetic. XCUR will be advanced * using the "long" algorithm in SP arithmetic. * * Mode 3 will be used on machines that don't meet the Mode 2 test, * but can maintain at least a 56 bits exactly in computing */ /* mod(AFAC*XCUR, MDIV) in DP arithmetic. This includes VAX, UNISYS, * IBM 30xx, and some IEEE machines that have clever compilers that * keep an extended precision representation of the product AFAC*XCUR * within the math processor for use in the division by MDIV. XCUR will * be advanced using the "short" algorithm in DP arithmetic. * * Mode 4 will be used on machines that don't meet the Mode 2 or 3 * tests, but have at least a 38 bit significand in DP arithmetic. * This includes IEEE machines that have not-so-clever compilers. * XCUR is advanced using the "long" algorithm in DP arithmetic. * --------------------------------------------------------------------- * Properties of the generated sequence. * * This m is one of the prime numbers cited in * Table 1, Page 390, of Knuth, Seminumerical Algorithms, Second * edition, 1981, Addison-Wesley. * The prime factorization of m-1 is * m-1 = p1 * p2 * p3 = 2 * 43801 * 784_451 * The complementary factors are * q(1) = 3_43597_38251, q(2) = 15_68902, and q(3) = 87602. * * The value a is a primitive root of m as is verified by * computing a**q(i) mod m for i = 1,3, and finding these values are * not 1. These values are m-1, 2_49653_21011, and 1_44431_31136. * The fact that a is a primitive root of m assures that the period * of the generator is m-1, i.e. starting with any integer from 1 * through m-1, all integers in this range will be produced. * * The value a has relatively large values of the measures nu and mu * computed for the Spectral Test as described in Knuth, pp. 89-105. * (Log10(nu(i)), i=2,6) = 5.4, 3.6, 2.6, 2.2, 1.8 * (mu(i), i=2,6) = 3.0, 3.05, 3.39, 4.55, 6.01 * This assures that the generated sequence will have relatively low * autocorrelation. * ------------------------------------------------------------------ * Alternative algorithm * * An alternative set of constants that has been used widely in * commercial and public domain software packages is * m = 21474_83647 = 2**31 - 1 * a = 16807 = 7**5 = (approx.) 0.513 * 2**15 * * The largest product that must be handled exactly is approximately * 0.513 * 2**46 which is approximately 0.36E14. This is within the * double-precision capability of most computer systems. * * The sequence can be started with any integer from 1 through m-1 * and will generate all integers in this range. The autocorrelation * properties of the whole sequence will not be as good as with the * larger values for m and a. * ------------------------------------------------------------------ * C. L. Lawson, F. T. Krogh & S. Chiu, JPL, July 1986, Apr 1987. * ------------------------------------------------------------------ * Common Block */ /* These same parameters are also defined below in RANMOD. */ /* ------------------------------------------------------------------ * Entered using CALL RN1 * This entry should not be called by general users. * User should call RAN1 in RANPK1. * */ rancom.xcurdp = X1DP; rancom.xcursp = X1SP; return; } /* end of function */ /* ------------------------------------------------------------------ */ void /*FUNCTION*/ ransiz( long *ksize) { *ksize = 2; return; } /* end of function */ /* ------------------------------------------------------------------ * Entered using CALL RNPUT(KSEED) * This entry should not be called by general users. * User should call RANPUT(KSEED) in RANPK1. * */ /* PARAMETER translations */ #define MDIVDP 68719476503.0e0 #define MDIVSP 68719476503.0e0 #define SCALDP 100000.0e0 #define SCALSP 100000.0e0 /* end of PARAMETER translations */ void /*FUNCTION*/ rnput( long kseed[]) { static LOGICAL32 first = TRUE; /* OFFSET Vectors w/subscript range: 1 to dimension */ long *const Kseed = &kseed[0] - 1; /* end of OFFSET VECTORS */ /* Common Block */ if (first) { first = FALSE; ranmod(); } if (rancom.mode == 2) { rancom.xcursp = SCALSP*(float)( labs( Kseed[1] ) ) + (float)( labs( Kseed[2] ) ); rancom.xcursp = fmaxf( 1.e0, fminf( rancom.xcursp, MDIVSP - 1.0e0 ) ); } else { /* Here for MODE = 3 or 4 */ rancom.xcurdp = SCALDP*(double)( labs( Kseed[1] ) ) + (double)( labs( Kseed[2] ) ); rancom.xcurdp = fmax( 1.e0, fmin( rancom.xcurdp, MDIVDP - 1.0e0 ) ); } return; } /* end of function */ /* ------------------------------------------------------------------ */ void /*FUNCTION*/ ranget( long kseed[]) { static LOGICAL32 first = TRUE; /* OFFSET Vectors w/subscript range: 1 to dimension */ long *const Kseed = &kseed[0] - 1; /* end of OFFSET VECTORS */ /* Common Block */ if (first) { first = FALSE; ranmod(); } if (rancom.mode == 2) { Kseed[1] = (long)( rancom.xcursp/SCALSP ); Kseed[2] = (long)( rancom.xcursp - SCALSP*(float)( Kseed[1] ) ); } else { /* Here for MODE = 3 or 4 * */ Kseed[1] = (long)( rancom.xcurdp/SCALDP ); Kseed[2] = (long)( rancom.xcurdp - SCALDP*(double)( Kseed[1] ) ); } return; } /* end of function */ /* ------------------------------------------------------------------ */ void /*FUNCTION*/ rn2( long *mode1) { static LOGICAL32 first = TRUE; /* Common Block */ if (first) { first = FALSE; ranmod(); } *mode1 = rancom.mode; return; } /* end of function */ /* ------------------------------------------------------------------ */ /* PARAMETER translations */ #define AFACDP 612662.0e0 #define AFACSP 612662.0e0 #define BFACDP 112165.0e0 #define BFACSP 112165.0e0 #define C2DP 243273.0e0 #define C2SP 243273.0e0 /* end of PARAMETER translations */ void /*FUNCTION*/ sranua( float usp[], long n) { long int i; float qsp, rsp; double qdp, rdp; static LOGICAL32 first = TRUE; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const Usp = &usp[0] - 1; /* end of OFFSET VECTORS */ /* Common Block */ if (first) { first = FALSE; ranmod(); } switch (rancom.mode) { case 1: goto L_310; case 2: goto L_320; case 3: goto L_330; case 4: goto L_340; } L_310: printf("In SRANUA inside RANPK2 invalid MODE =%10ld\n", rancom.mode); exit(0); /* Mode 2. */ L_320: ; for (i = 1; i <= n; i++) { qsp = trunc( rancom.xcursp/BFACSP ); rsp = rancom.xcursp - qsp*BFACSP; rancom.xcursp = AFACSP*rsp - C2SP*qsp; L_322: if (rancom.xcursp < 0.0e0) { rancom.xcursp += MDIVSP; goto L_322; } Usp[i] = rancom.xcursp/MDIVSP; } goto L_350; /* Mode 3. */ L_330: ; for (i = 1; i <= n; i++) { rancom.xcurdp = fmod( AFACDP*rancom.xcurdp, MDIVDP ); Usp[i] = (float)( rancom.xcurdp )/MDIVSP; } goto L_350; /* Mode 4. */ L_340: ; for (i = 1; i <= n; i++) { qdp = trunc( rancom.xcurdp/BFACDP ); rdp = rancom.xcurdp - qdp*BFACDP; rancom.xcurdp = AFACDP*rdp - C2DP*qdp; L_342: if (rancom.xcurdp < 0.0e0) { rancom.xcurdp += MDIVDP; goto L_342; } Usp[i] = (float)( rancom.xcurdp )/MDIVSP; } L_350: ; return; } /* end of function */ /* ------------------------------------------------------------------ */ void /*FUNCTION*/ dranua( double udp[], long n) { long int i; float qsp, rsp; double qdp, rdp; static LOGICAL32 first = TRUE; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Udp = &udp[0] - 1; /* end of OFFSET VECTORS */ /* Common Block */ if (first) { first = FALSE; ranmod(); } switch (rancom.mode) { case 1: goto L_410; case 2: goto L_420; case 3: goto L_430; case 4: goto L_440; } L_410: printf("In DRANUA inside RANPK2 invalid MODE =%10ld\n", rancom.mode); exit(0); /* Mode 2. */ L_420: ; for (i = 1; i <= n; i++) { qsp = trunc( rancom.xcursp/BFACSP ); rsp = rancom.xcursp - qsp*BFACSP; rancom.xcursp = AFACSP*rsp - C2SP*qsp; L_422: if (rancom.xcursp < 0.0e0) { rancom.xcursp += MDIVSP; goto L_422; } Udp[i] = (double)( rancom.xcursp )/MDIVDP; } goto L_450; /* Mode 3. */ L_430: ; for (i = 1; i <= n; i++) { rancom.xcurdp = fmod( AFACDP*rancom.xcurdp, MDIVDP ); Udp[i] = rancom.xcurdp/MDIVDP; } goto L_450; /* Mode 4. */ L_440: ; for (i = 1; i <= n; i++) { qdp = trunc( rancom.xcurdp/BFACDP ); rdp = rancom.xcurdp - qdp*BFACDP; rancom.xcurdp = AFACDP*rdp - C2DP*qdp; L_442: if (rancom.xcurdp < 0.0e0) { rancom.xcurdp += MDIVDP; goto L_442; } Udp[i] = rancom.xcurdp/MDIVDP; } L_450: ; return; } /* end of function */ /* ------------------------------------------------------------------ */ void /*FUNCTION*/ sranus( float usp[], long n, float asp, float bsp) { long int i; float qsp, rsp; double qdp, rdp; static LOGICAL32 first = TRUE; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const Usp = &usp[0] - 1; /* end of OFFSET VECTORS */ /* Common Block */ if (first) { first = FALSE; ranmod(); } switch (rancom.mode) { case 1: goto L_510; case 2: goto L_520; case 3: goto L_530; case 4: goto L_540; } L_510: puts( "In file RANPK2, subroutine SRANUS -- Invalid value for MODE" ); exit(0); /* Mode 2. */ L_520: ; for (i = 1; i <= n; i++) { qsp = trunc( rancom.xcursp/BFACSP ); rsp = rancom.xcursp - qsp*BFACSP; rancom.xcursp = AFACSP*rsp - C2SP*qsp; L_522: if (rancom.xcursp < 0.0e0) { rancom.xcursp += MDIVSP; goto L_522; } Usp[i] = asp + bsp*rancom.xcursp/MDIVSP; } goto L_550; /* Mode 3. */ L_530: ; for (i = 1; i <= n; i++) { rancom.xcurdp = fmod( AFACDP*rancom.xcurdp, MDIVDP ); Usp[i] = asp + bsp*(float)( rancom.xcurdp )/MDIVSP; } goto L_550; /* Mode 4. */ L_540: ; for (i = 1; i <= n; i++) { qdp = trunc( rancom.xcurdp/BFACDP ); rdp = rancom.xcurdp - qdp*BFACDP; rancom.xcurdp = AFACDP*rdp - C2DP*qdp; L_542: if (rancom.xcurdp < 0.0e0) { rancom.xcurdp += MDIVDP; goto L_542; } Usp[i] = asp + bsp*(float)( rancom.xcurdp )/MDIVSP; } L_550: ; return; } /* end of function */ /* ------------------------------------------------------------------ */ void /*FUNCTION*/ dranus( double udp[], long n, double adp, double bdp) { long int i; float qsp, rsp; double qdp, rdp; static LOGICAL32 first = TRUE; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Udp = &udp[0] - 1; /* end of OFFSET VECTORS */ /* Common Block */ if (first) { first = FALSE; ranmod(); } switch (rancom.mode) { case 1: goto L_610; case 2: goto L_620; case 3: goto L_630; case 4: goto L_640; } L_610: puts( "In file RANPK2, subroutine DRANUS -- Invalid value for MODE" ); exit(0); /* Mode 2. */ L_620: ; for (i = 1; i <= n; i++) { qsp = trunc( rancom.xcursp/BFACSP ); rsp = rancom.xcursp - qsp*BFACSP; rancom.xcursp = AFACSP*rsp - C2SP*qsp; L_622: if (rancom.xcursp < 0.0e0) { rancom.xcursp += MDIVSP; goto L_622; } Udp[i] = adp + bdp*(double)( rancom.xcursp )/MDIVDP; } goto L_650; /* Mode 3. */ L_630: ; for (i = 1; i <= n; i++) { rancom.xcurdp = fmod( AFACDP*rancom.xcurdp, MDIVDP ); Udp[i] = adp + bdp*rancom.xcurdp/MDIVDP; } goto L_650; /* Mode 4. */ L_640: ; for (i = 1; i <= n; i++) { qdp = trunc( rancom.xcurdp/BFACDP ); rdp = rancom.xcurdp - qdp*BFACDP; rancom.xcurdp = AFACDP*rdp - C2DP*qdp; L_642: if (rancom.xcurdp < 0.0e0) { rancom.xcurdp += MDIVDP; goto L_642; } Udp[i] = adp + bdp*rancom.xcurdp/MDIVDP; } L_650: ; return; } /* end of function */ /* ------------------------------------------------------------------ */ void /*FUNCTION*/ ranmod() { long int i, j, k; float qsp, rsp, xtstsp; double diff, qdp, rdp, xtstdp; static LOGICAL32 done = FALSE; static double test[1-(0)+1][3-(0)+1]={24997965550.0e0,68719476502.0e0, 68718863841.0e0,36962132774.0e0,43721510953.0e0,1.0e0,612662.0e0, 31757343729.0e0}; /* Do tests to decide whether to set the MODE = 2, 3, or 4. * Outcome will depend on precision of SP and DP floating point * arithmetic on the host system. * Common Block */ if (done) return; done = TRUE; for (rancom.mode = 2; rancom.mode <= 4; rancom.mode++) { for (j = 0; j <= 1; j++) { xtstdp = test[j][0]; xtstsp = (float)( xtstdp ); for (i = 1; i <= 3; i++) { switch (rancom.mode - 1) { case 1: goto L_820; case 2: goto L_830; case 3: goto L_840; } /* Test of MODE 2 */ L_820: ; qsp = trunc( xtstsp/BFACSP ); rsp = xtstsp - qsp*BFACSP; xtstsp = AFACSP*rsp - C2SP*qsp; for (k = 1; k <= 3; k++) { if (xtstsp >= 0.0e0) goto L_825; xtstsp += MDIVSP; } L_825: ; diff = (double)( xtstsp ) - test[j][i]; goto L_850; /* Test of MODE 3 */ L_830: ; xtstdp = fmod( AFACDP*xtstdp, MDIVDP ); diff = xtstdp - test[j][i]; goto L_850; /* Test of MODE 4 */ L_840: ; qdp = trunc( xtstdp/BFACDP ); rdp = xtstdp - qdp*BFACDP; xtstdp = AFACDP*rdp - C2DP*qdp; for (k = 1; k <= 3; k++) { if (xtstdp >= 0.0e0) goto L_845; xtstdp += MDIVDP; } L_845: ; diff = xtstdp - test[j][i]; L_850: ; /* print'(1x,a,3i3,g11.3)', 'RANPK2.. MODE,J,I,DIFF=', * * MODE,J,I,DIFF !****** For Testing ****** */ if (diff != 0.0e0) goto L_880; /* Following line ends I loop. */ } /* Following line ends J loop. */ } /* Here the computations using the current value of MODE have * passed all tests, so we accept this value of MODE. */ rancom.xcurdp = X1DP; rancom.xcursp = X1SP; return; /* print*,'From RANPK2.. MODE =',MODE !***** For Testing ***** * Following line ends MODE loop. */ L_880: ; } /* The computations were unsuccessful for all values of MODE. * This means this random number package will not work on the * current host system. ****** Fatal Error Stop ****** * */ ermsg( "RANPK2", 1, 2, "This random no. code will not work on this computer system." , '.' ); return; } /* end of function */