htswanalysis/MACS/lib/gsl/gsl-1.11/qrng/sobol.c

   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 <config.h>
   9 #include <gsl/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
 129 static size_t sobol_state_size(unsigned int dimension)
 130 {
 131   return sizeof(sobol_state_t);
 132 }
 133
 134
 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
 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 }