| 1030 | /******************************************************************************/ |
| 1031 | |
| 1032 | double r8_random ( int iseed[4] ) |
| 1033 | |
| 1034 | /******************************************************************************/ |
| 1035 | /* |
| 1036 | Purpose: |
| 1037 | |
| 1038 | R8_RANDOM returns a uniformly distributed random number between 0 and 1. |
| 1039 | |
| 1040 | Discussion: |
| 1041 | |
| 1042 | This routine uses a multiplicative congruential method with modulus |
| 1043 | 2**48 and multiplier 33952834046453 (see G.S.Fishman, |
| 1044 | 'Multiplicative congruential random number generators with modulus |
| 1045 | 2**b: an exhaustive analysis for b = 32 and a partial analysis for |
| 1046 | b = 48', Math. Comp. 189, pp 331-344, 1990). |
| 1047 | |
| 1048 | 48-bit integers are stored in 4 integer array elements with 12 bits |
| 1049 | per element. Hence the routine is portable across machines with |
| 1050 | integers of 32 bits or more. |
| 1051 | |
| 1052 | Parameters: |
| 1053 | |
| 1054 | Input/output, integer ISEED(4). |
| 1055 | On entry, the seed of the random number generator; the array |
| 1056 | elements must be between 0 and 4095, and ISEED(4) must be odd. |
| 1057 | On exit, the seed is updated. |
| 1058 | |
| 1059 | Output, double R8_RANDOM, the next pseudorandom number. |
| 1060 | */ |
| 1061 | { |
| 1062 | int ipw2 = 4096; |
| 1063 | int it1; |
| 1064 | int it2; |
| 1065 | int it3; |
| 1066 | int it4; |
| 1067 | int m1 = 494; |
| 1068 | int m2 = 322; |
| 1069 | int m3 = 2508; |
| 1070 | int m4 = 2549; |
| 1071 | double one = 1.0; |
| 1072 | double r = 1.0 / 4096.0; |
| 1073 | double value; |
| 1074 | /* |
| 1075 | Multiply the seed by the multiplier modulo 2**48. |
| 1076 | */ |
| 1077 | it4 = iseed[3] * m4; |
| 1078 | it3 = it4 / ipw2; |
| 1079 | it4 = it4 - ipw2 * it3; |
| 1080 | it3 = it3 + iseed[2] * m4 + iseed[3] * m3; |
| 1081 | it2 = it3 / ipw2; |
| 1082 | it3 = it3 - ipw2 * it2; |
| 1083 | it2 = it2 + iseed[1] * m4 + iseed[2] * m3 + iseed[3] * m2; |
| 1084 | it1 = it2 / ipw2; |
| 1085 | it2 = it2 - ipw2 * it1; |
| 1086 | it1 = it1 + iseed[0] * m4 + iseed[1] * m3 + iseed[2] * m2 + iseed[3] * m1; |
| 1087 | it1 = ( it1 % ipw2 ); |
| 1088 | /* |
| 1089 | Return updated seed |