File indexing completed on 2023-10-25 10:06:57
0001
0002
0003 #include <stdlib.h>
0004
0005 #include "Random.h"
0006
0007 #ifndef NULL
0008 #define NULL 0
0009 #endif
0010
0011
0012 #define MDIG 32
0013
0014
0015 #define ONE 1
0016
0017 static const int m1 = (ONE << (MDIG - 2)) + ((ONE << (MDIG - 2)) - ONE);
0018 static const int m2 = ONE << MDIG / 2;
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028 static double dm1;
0029
0030
0031
0032 static void initialize(Random R, int seed);
0033
0034 Random new_Random_seed(int seed) {
0035 Random R = (Random)malloc(sizeof(Random_struct));
0036
0037 initialize(R, seed);
0038 R->left = 0.0;
0039 R->right = 1.0;
0040 R->width = 1.0;
0041 R->haveRange = 0 ;
0042
0043 return R;
0044 }
0045
0046 Random new_Random(int seed, double left, double right) {
0047 Random R = (Random)malloc(sizeof(Random_struct));
0048
0049 initialize(R, seed);
0050 R->left = left;
0051 R->right = right;
0052 R->width = right - left;
0053 R->haveRange = 1;
0054
0055 return R;
0056 }
0057
0058 void Random_delete(Random R) { free(R); }
0059
0060
0061
0062 double Random_nextDouble(Random R) {
0063 int k;
0064
0065 int I = R->i;
0066 int J = R->j;
0067 int *m = R->m;
0068
0069 k = m[I] - m[J];
0070 if (k < 0)
0071 k += m1;
0072 R->m[J] = k;
0073
0074 if (I == 0)
0075 I = 16;
0076 else
0077 I--;
0078 R->i = I;
0079
0080 if (J == 0)
0081 J = 16;
0082 else
0083 J--;
0084 R->j = J;
0085
0086 if (R->haveRange)
0087 return R->left + dm1 * (double)k * R->width;
0088 else
0089 return dm1 * (double)k;
0090 }
0091
0092
0093
0094
0095
0096 static void initialize(Random R, int seed) {
0097 int jseed, k0, k1, j0, j1, iloop;
0098
0099 dm1 = 1.0 / (double)m1;
0100
0101 R->seed = seed;
0102
0103 if (seed < 0)
0104 seed = -seed;
0105 jseed = (seed < m1 ? seed : m1);
0106 if (jseed % 2 == 0)
0107 --jseed;
0108 k0 = 9069 % m2;
0109 k1 = 9069 / m2;
0110 j0 = jseed % m2;
0111 j1 = jseed / m2;
0112 for (iloop = 0; iloop < 17; ++iloop) {
0113 jseed = j0 * k0;
0114 j1 = (jseed / m2 + j0 * k1 + j1 * k0) % (m2 / 2);
0115 j0 = jseed % m2;
0116 R->m[iloop] = j0 + m2 * j1;
0117 }
0118 R->i = 4;
0119 R->j = 16;
0120 }
0121
0122 double *RandomVector(int N, Random R) {
0123 int i;
0124 double *x = (double *)malloc(sizeof(double) * N);
0125
0126 for (i = 0; i < N; i++)
0127 x[i] = Random_nextDouble(R);
0128
0129 return x;
0130 }
0131
0132 double **RandomMatrix(int M, int N, Random R) {
0133 int i;
0134 int j;
0135
0136
0137
0138 double **A = (double **)malloc(sizeof(double *) * M);
0139
0140 if (A == NULL)
0141 return NULL;
0142
0143 for (i = 0; i < M; i++) {
0144 A[i] = (double *)malloc(sizeof(double) * N);
0145 if (A[i] == NULL) {
0146 free(A);
0147 return NULL;
0148 }
0149 for (j = 0; j < N; j++)
0150 A[i][j] = Random_nextDouble(R);
0151 }
0152 return A;
0153 }