1 /* Author: G. Jungman
2 */
3 /* Implementation for Sobol generator.
4 * See
5 * [Bratley+Fox, TOMS 14, 88 (1988)]
6 * [Antonov+Saleev, USSR Comput. Maths. Math. Phys. 19, 252 (1980)]
7 */
8 #include "gsl__config.h"
9 #include "gsl_qrng.h"
10
11
12 /* maximum allowed space dimension */
13 #define SOBOL_MAX_DIMENSION 40
14
15 /* bit count; assumes sizeof(int) >= 32-bit */
16 #define SOBOL_BIT_COUNT 30
17
18 /* prototypes for generator type functions */
19 static size_t sobol_state_size(unsigned int dimension);
20 static int sobol_init(void * state, unsigned int dimension);
21 static int sobol_get(void * state, unsigned int dimension, double * v);
22
23 /* global Sobol generator type object */
24 static const gsl_qrng_type sobol_type =
25 {
26 "sobol",
27 SOBOL_MAX_DIMENSION,
28 sobol_state_size,
29 sobol_init,
30 sobol_get
31 };
32 const gsl_qrng_type * gsl_qrng_sobol = &sobol_type;
33
34
35 /* primitive polynomials in binary encoding
36 */
37 static const int primitive_polynomials[SOBOL_MAX_DIMENSION] =
38 {
39 1, 3, 7, 11, 13, 19, 25, 37, 59, 47,
40 61, 55, 41, 67, 97, 91, 109, 103, 115, 131,
41 193, 137, 145, 143, 241, 157, 185, 167, 229, 171,
42 213, 191, 253, 203, 211, 239, 247, 285, 369, 299
43 };
44
45 /* degrees of the primitive polynomials */
46 static const int degree_table[SOBOL_MAX_DIMENSION] =
47 {
48 0, 1, 2, 3, 3, 4, 4, 5, 5, 5,
49 5, 5, 5, 6, 6, 6, 6, 6, 6, 7,
50 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
51 7, 7, 7, 7, 7, 7, 7, 8, 8, 8
52 };
53
54
55 /* initial values for direction tables, following
56 * Bratley+Fox, taken from [Sobol+Levitan, preprint 1976]
57 */
58 static const int v_init[8][SOBOL_MAX_DIMENSION] =
59 {
60 {
61 0, 1, 1, 1, 1, 1, 1, 1, 1, 1,
62 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
63 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
64 1, 1, 1, 1, 1, 1, 1, 1, 1, 1
65 },
66 {
67 0, 0, 1, 3, 1, 3, 1, 3, 3, 1,
68 3, 1, 3, 1, 3, 1, 1, 3, 1, 3,
69 1, 3, 1, 3, 3, 1, 3, 1, 3, 1,
70 3, 1, 1, 3, 1, 3, 1, 3, 1, 3
71 },
72 {
73 0, 0, 0, 7, 5, 1, 3, 3, 7, 5,
74 5, 7, 7, 1, 3, 3, 7, 5, 1, 1,
75 5, 3, 3, 1, 7, 5, 1, 3, 3, 7,
76 5, 1, 1, 5, 7, 7, 5, 1, 3, 3
77 },
78 {
79 0, 0, 0, 0, 0, 1, 7, 9, 13, 11,
80 1, 3, 7, 9, 5, 13, 13, 11, 3, 15,
81 5, 3, 15, 7, 9, 13, 9, 1, 11, 7,
82 5, 15, 1, 15, 11, 5, 3, 1, 7, 9
83 },
84 {
85 0, 0, 0, 0, 0, 0, 0, 9, 3, 27,
86 15, 29, 21, 23, 19, 11, 25, 7, 13, 17,
87 1, 25, 29, 3, 31, 11, 5, 23, 27, 19,
88 21, 5, 1, 17, 13, 7, 15, 9, 31, 9
89 },
90 {
91 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
92 0, 0, 0, 37, 33, 7, 5, 11, 39, 63,
93 27, 17, 15, 23, 29, 3, 21, 13, 31, 25,
94 9, 49, 33, 19, 29, 11, 19, 27, 15, 25
95 },
96 {
97 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
98 0, 0, 0, 0, 0, 0, 0, 0, 0, 13,
99 33, 115, 41, 79, 17, 29, 119, 75, 73, 105,
100 7, 59, 65, 21, 3, 113, 61, 89, 45, 107
101 },
102 {
103 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
104 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
105 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
106 0, 0, 0, 0, 0, 0, 0, 7, 23, 39
107 }
108 };
109
110
111 /* Sobol generator state.
112 * sequence_count = number of calls with this generator
113 * last_numerator_vec = last generated numerator vector
114 * last_denominator_inv = 1/denominator for last numerator vector
115 * v_direction = direction number table
116 */
117 typedef struct
118 {
119 unsigned int sequence_count;
120 double last_denominator_inv;
121 int last_numerator_vec[SOBOL_MAX_DIMENSION];
122 int v_direction[SOBOL_BIT_COUNT][SOBOL_MAX_DIMENSION];
123 } sobol_state_t;
124
125 /* NOTE: I fixed the size for the preliminary implementation.
126 This could/should be relaxed, as an optimization.
127 */
128
sobol_state_size(unsigned int dimension)129 static size_t sobol_state_size(unsigned int dimension)
130 {
131 return sizeof(sobol_state_t);
132 }
133
134
sobol_init(void * state,unsigned int dimension)135 static int sobol_init(void * state, unsigned int dimension)
136 {
137 sobol_state_t * s_state = (sobol_state_t *) state;
138 unsigned int i_dim;
139 int j, k;
140 int ell;
141
142 if(dimension < 1 || dimension > SOBOL_MAX_DIMENSION) {
143 return GSL_EINVAL;
144 }
145
146 /* Initialize direction table in dimension 0. */
147 for(k=0; k<SOBOL_BIT_COUNT; k++) s_state->v_direction[k][0] = 1;
148
149 /* Initialize in remaining dimensions. */
150 for(i_dim=1; i_dim<dimension; i_dim++) {
151
152 const int poly_index = i_dim;
153 const int degree_i = degree_table[poly_index];
154 int includ[8];
155
156 /* Expand the polynomial bit pattern to separate
157 * components of the logical array includ[].
158 */
159 int p_i = primitive_polynomials[poly_index];
160 for(k = degree_i-1; k >= 0; k--) {
161 includ[k] = ((p_i % 2) == 1);
162 p_i /= 2;
163 }
164
165 /* Leading elements for dimension i come from v_init[][]. */
166 for(j=0; j<degree_i; j++) s_state->v_direction[j][i_dim] = v_init[j][i_dim];
167
168 /* Calculate remaining elements for this dimension,
169 * as explained in Bratley+Fox, section 2.
170 */
171 for(j=degree_i; j<SOBOL_BIT_COUNT; j++) {
172 int newv = s_state->v_direction[j-degree_i][i_dim];
173 ell = 1;
174 for(k=0; k<degree_i; k++) {
175 ell *= 2;
176 if(includ[k]) newv ^= (ell * s_state->v_direction[j-k-1][i_dim]);
177 }
178 s_state->v_direction[j][i_dim] = newv;
179 }
180 }
181
182 /* Multiply columns of v by appropriate power of 2. */
183 ell = 1;
184 for(j=SOBOL_BIT_COUNT-1-1; j>=0; j--) {
185 ell *= 2;
186 for(i_dim=0; i_dim<dimension; i_dim++) {
187 s_state->v_direction[j][i_dim] *= ell;
188 }
189 }
190
191 /* 1/(common denominator of the elements in v_direction) */
192 s_state->last_denominator_inv = 1.0 /(2.0 * ell);
193
194 /* final setup */
195 s_state->sequence_count = 0;
196 for(i_dim=0; i_dim<dimension; i_dim++) s_state->last_numerator_vec[i_dim] = 0;
197
198 return GSL_SUCCESS;
199 }
200
201
sobol_get(void * state,unsigned int dimension,double * v)202 static int sobol_get(void * state, unsigned int dimension, double * v)
203 {
204 sobol_state_t * s_state = (sobol_state_t *) state;
205
206 unsigned int i_dimension;
207
208 /* Find the position of the least-significant zero in count. */
209 int ell = 0;
210 int c = s_state->sequence_count;
211 while(1) {
212 ++ell;
213 if((c % 2) == 1) c /= 2;
214 else break;
215 }
216
217 /* Check for exhaustion. */
218 if(ell > SOBOL_BIT_COUNT) return GSL_EFAILED; /* FIXME: good return code here */
219
220 for(i_dimension=0; i_dimension<dimension; i_dimension++) {
221 const int direction_i = s_state->v_direction[ell-1][i_dimension];
222 const int old_numerator_i = s_state->last_numerator_vec[i_dimension];
223 const int new_numerator_i = old_numerator_i ^ direction_i;
224 s_state->last_numerator_vec[i_dimension] = new_numerator_i;
225 v[i_dimension] = new_numerator_i * s_state->last_denominator_inv;
226 }
227
228 s_state->sequence_count++;
229
230 return GSL_SUCCESS;
231 }
232