| c end PRanf |
| |
| c============================================================================== |
| |
| subroutine Rans( NIn, Seed1, SeedArray ) |
| |
| c------------------------------------------------------------------------------ |
| c |
| c This routine divides the sequence of random numbers into N subsequences, |
| c each with its own seed. The seeds for the independent subsequences are |
| c returned in the seed array. if Seed1 is zero, all zeroes will be returned. |
| c To prevent this, Seed1 is set to [3281, 4041, 595, 2376], which is |
| c statistically the best starting seed. The wheel is then divided into the |
| c N pieces (where N is odd and >= NIn) by dividing its period (2**46) by N. |
| c |
| c Then, seed(i) = seed(i-1) * (a**k mod 2**48), and 1<=k<=N. |
| c |
| c Here, 'a' is the multiplier used by the linear congruential generator whose |
| c wheel we are dividing up. |
| c |
| c The number of streams must be odd; if NIn is even N will be NIn+1, and |
| c n extra stream of random numbers will be available that will not get used. |
| c |
| c It returns an array of seeds, each an array of 4 integers that are used as |
| c the digits of a four-digit modulo-4096 integer. |
| c |
| c------------------------------------------------------------------------------ |
| |
| include 'pranf.inc' |
| |
| integer NIn, Seed1, DefaultValues, SeedArray, N, atothek, K, |
| 1 KBinary, InSeed, OutSeed |
| |
| c****************************************************************************** |
| c**** Data common to the PRanf package. |
| integer Multiplier, DefaultSeed |
| |
| real*8 Divisor |
| |
| dimension Multiplier( IRandNumSize ), DefaultSeed( IRandNumSize ), |
| 1 Divisor( IRandNumSize ) |
| |
| common / PRanfCom / Multiplier, DefaultSeed, Divisor |
| c**** End of PRanf common data |
| c****************************************************************************** |
| |
| dimension DefaultValues( IRandNumSize ), |
| 1 SeedArray( IRandNumSize, MaxStreams ), |
| 2 atothek( IRandNumSize ), |
| 3 K( IRandNumSize ), |
| 4 KBinary( IBinarySize ), |
| 5 InSeed( IRandNumSize ), |
| 6 OutSeed( IRandNumSize ) |
| |
| INTEGER NisOdd |
| |
| c scatter, shared SeedArray |
| |
| c------------------------------------------------------------------------------ |
| |
| |
| c write( *, 10 ) NIn, Seed1 |
| 10 format( ' Rans: NIn = ', i6, ' Seed1 = ', i6, /, |
| 1 '------' ) |
| |
| IZero = 0 |
| |
| c**** Make sure we are generating an odd number of random number sequences. |
| if( NisOdd( NIn ) .eq. 1 ) then |
| N = NIN |
| else |
| N = NIn + 1 |
| endif |
| |
| c**** Set up the initial seed to either a legal input value or its default |
| c**** values. The input integer, if nonzero, is used for the first of the |
| c**** four modulo-4096 digits of the actual initial seed. |
| |
| c**** The First element of the returned SeedArray will be used here, since |
| c**** at least one seed will be returned, and the first seed of the set |
| c**** returned will be the seed for the entire wheel. |
| if( Seed1 .eq. IZero ) then |
| SeedArray( 1, 1 ) = DefaultSeed( 1 ) |
| SeedArray( 2, 1 ) = DefaultSeed( 2 ) |
| SeedArray( 3, 1 ) = DefaultSeed( 3 ) |
| SeedArray( 4, 1 ) = DefaultSeed( 4 ) |
| else |
| SeedArray( 1, 1 ) = abs( Seed1 ) |
| SeedArray( 2, 1 ) = IZero |
| SeedArray( 3, 1 ) = IZero |
| SeedArray( 4, 1 ) = IZero |
| endif |
| |
| c**** 'a' is the multiplier for the Ranf linear congruential generator. |
| if( N .eq. 1 ) then |
| c*******If only one stream of random numbers is needed, do not bother to |
| c*******raise 'a' to the first power. |
| atothek( 1 ) = Multiplier( 1 ) |
| atothek( 2 ) = Multiplier( 2 ) |
| atothek( 3 ) = Multiplier( 3 ) |
| atothek( 4 ) = Multiplier( 4 ) |
| else |
| c*******more than one stream is needed; generate the Kth seed by multiplying |
| c*******the K-1st seed by the multiplier raised to the Nth power. |
| CALL ranfk( N, K ) |
| CALL ranfkbinary( K, KBinary ) |
| CALL ranfatok( Multiplier, KBinary, atothek ) |
| do 100 I = 2, N |
| InSeed( 1 ) = SeedArray( 1, I-1 ) |
| InSeed( 2 ) = SeedArray( 2, I-1 ) |
| InSeed( 3 ) = SeedArray( 3, I-1 ) |
| InSeed( 4 ) = SeedArray( 4, I-1 ) |
| CALL ranfmodmult( InSeed, atothek, OutSeed ) |
| SeedArray( 1, I ) = OutSeed( 1 ) |
| SeedArray( 2, I ) = OutSeed( 2 ) |
| SeedArray( 3, I ) = OutSeed( 3 ) |
| SeedArray( 4, I ) = OutSeed( 4 ) |
| 100 continue |
| endif |
| |
| c write( *, 110 ) |
| 110 format( ' leaving Rans' ) |
| |
| return |
| end |
| |
| |
| c MJC mod to remove statement function |
| INTEGER FUNCTION NISODD( M ) |
| c**** check for argument numbers being odd. |
| NIsOdd = mod( m, 2 ) |
| RETURN |
| END |