htswanalysis/MACS/lib/gsl/gsl-1.11/qrng/niederreiter-2.c

   1 /* Author: G. Jungman
   2  */
   3 /* Implement Niederreiter base 2 generator.
   4  * See:
   5  *   Bratley, Fox, Niederreiter, ACM Trans. Model. Comp. Sim. 2, 195 (1992)
   6  */
   7 #include <config.h>
   8 #include <gsl/gsl_qrng.h>
   9
  10
  11 #define NIED2_CHARACTERISTIC 2
  12 #define NIED2_BASE 2
  13
  14 #define NIED2_MAX_DIMENSION 12
  15 #define NIED2_MAX_PRIM_DEGREE 5
  16 #define NIED2_MAX_DEGREE 50
  17
  18 #define NIED2_BIT_COUNT 30
  19 #define NIED2_NBITS (NIED2_BIT_COUNT+1)
  20
  21 #define MAXV (NIED2_NBITS + NIED2_MAX_DEGREE)
  22
  23 /* Z_2 field operations */
  24 #define NIED2_ADD(x,y) (((x)+(y))%2)
  25 #define NIED2_MUL(x,y) (((x)*(y))%2)
  26 #define NIED2_SUB(x,y) NIED2_ADD((x),(y))
  27
  28
  29 static size_t nied2_state_size(unsigned int dimension);
  30 static int nied2_init(void * state, unsigned int dimension);
  31 static int nied2_get(void * state, unsigned int dimension, double * v);
  32
  33
  34 static const gsl_qrng_type nied2_type =
  35 {
  36   "niederreiter-base-2",
  37   NIED2_MAX_DIMENSION,
  38   nied2_state_size,
  39   nied2_init,
  40   nied2_get
  41 };
  42
  43 const gsl_qrng_type * gsl_qrng_niederreiter_2 = &nied2_type;
  44
  45 /* primitive polynomials in binary encoding */
  46 static const int primitive_poly[NIED2_MAX_DIMENSION+1][NIED2_MAX_PRIM_DEGREE+1] =
  47 {
  48   { 1, 0, 0, 0, 0, 0 },  /*  1               */
  49   { 0, 1, 0, 0, 0, 0 },  /*  x               */
  50   { 1, 1, 0, 0, 0, 0 },  /*  1 + x           */
  51   { 1, 1, 1, 0, 0, 0 },  /*  1 + x + x^2     */
  52   { 1, 1, 0, 1, 0, 0 },  /*  1 + x + x^3     */
  53   { 1, 0, 1, 1, 0, 0 },  /*  1 + x^2 + x^3   */
  54   { 1, 1, 0, 0, 1, 0 },  /*  1 + x + x^4     */
  55   { 1, 0, 0, 1, 1, 0 },  /*  1 + x^3 + x^4   */
  56   { 1, 1, 1, 1, 1, 0 },  /*  1 + x + x^2 + x^3 + x^4   */
  57   { 1, 0, 1, 0, 0, 1 },  /*  1 + x^2 + x^5             */
  58   { 1, 0, 0, 1, 0, 1 },  /*  1 + x^3 + x^5             */
  59   { 1, 1, 1, 1, 0, 1 },  /*  1 + x + x^2 + x^3 + x^5   */
  60   { 1, 1, 1, 0, 1, 1 }   /*  1 + x + x^2 + x^4 + x^5   */
  61 };
  62
  63 /* degrees of primitive polynomials */
  64 static const int poly_degree[NIED2_MAX_DIMENSION+1] =
  65 {
  66   0, 1, 1, 2, 3, 3, 4, 4, 4, 5, 5, 5, 5
  67 };
  68
  69
  70 typedef struct
  71 {
  72   unsigned int sequence_count;
  73   int cj[NIED2_NBITS][NIED2_MAX_DIMENSION];
  74   int nextq[NIED2_MAX_DIMENSION];
  75 } nied2_state_t;
  76
  77
  78 static size_t nied2_state_size(unsigned int dimension)
  79 {
  80   return sizeof(nied2_state_t);
  81 }
  82
  83
  84 /* Multiply polynomials over Z_2.
  85  * Notice use of a temporary vector,
  86  * side-stepping aliasing issues when
  87  * one of inputs is the same as the output
  88  * [especially important in the original fortran version, I guess].
  89  */
  90 static void poly_multiply(
  91   const int pa[], int pa_degree,
  92   const int pb[], int pb_degree,
  93   int pc[], int  * pc_degree
  94   )
  95 {
  96   int j, k;
  97   int pt[NIED2_MAX_DEGREE+1];
  98   int pt_degree = pa_degree + pb_degree;
  99
 100   for(k=0; k<=pt_degree; k++) {
 101     int term = 0;
 102     for(j=0; j<=k; j++) {
 103       const int conv_term = NIED2_MUL(pa[k-j], pb[j]);
 104       term = NIED2_ADD(term, conv_term);
 105     }
 106     pt[k] = term;
 107   }
 108
 109   for(k=0; k<=pt_degree; k++) {
 110     pc[k] = pt[k];
 111   }
 112   for(k=pt_degree+1; k<=NIED2_MAX_DEGREE; k++) {
 113     pc[k] = 0;
 114   }
 115
 116   *pc_degree = pt_degree;
 117 }
 118
 119
 120 /* Calculate the values of the constants V(J,R) as
 121  * described in BFN section 3.3.
 122  *
 123  *   px = appropriate irreducible polynomial for current dimension
 124  *   pb = polynomial defined in section 2.3 of BFN.
 125  * pb is modified
 126  */
 127 static void calculate_v(
 128   const int px[], int px_degree,
 129   int pb[], int * pb_degree,
 130   int v[], int maxv
 131   )
 132 {
 133   const int nonzero_element = 1;    /* nonzero element of Z_2  */
 134   const int arbitrary_element = 1;  /* arbitray element of Z_2 */
 135
 136   /* The polynomial ph is px**(J-1), which is the value of B on arrival.
 137    * In section 3.3, the values of Hi are defined with a minus sign:
 138    * don't forget this if you use them later !
 139    */
 140   int ph[NIED2_MAX_DEGREE+1];
 141   /* int ph_degree = *pb_degree; */
 142   int bigm = *pb_degree;      /* m from section 3.3 */
 143   int m;                      /* m from section 2.3 */
 144   int r, k, kj;
 145
 146   for(k=0; k<=NIED2_MAX_DEGREE; k++) {
 147     ph[k] = pb[k];
 148   }
 149
 150   /* Now multiply B by PX so B becomes PX**J.
 151    * In section 2.3, the values of Bi are defined with a minus sign :
 152    * don't forget this if you use them later !
 153    */
 154    poly_multiply(px, px_degree, pb, *pb_degree, pb, pb_degree);
 155    m = *pb_degree;
 156
 157   /* Now choose a value of Kj as defined in section 3.3.
 158    * We must have 0 <= Kj < E*J = M.
 159    * The limit condition on Kj does not seem very relevant
 160    * in this program.
 161    */
 162   /* Quoting from BFN: "Our program currently sets each K_q
 163    * equal to eq. This has the effect of setting all unrestricted
 164    * values of v to 1."
 165    * Actually, it sets them to the arbitrary chosen value.
 166    * Whatever.
 167    */
 168   kj = bigm;
 169
 170   /* Now choose values of V in accordance with
 171    * the conditions in section 3.3.
 172    */
 173   for(r=0; r<kj; r++) {
 174     v[r] = 0;
 175   }
 176   v[kj] = 1;
 177
 178
 179   if(kj >= bigm) {
 180     for(r=kj+1; r<m; r++) {
 181       v[r] = arbitrary_element;
 182     }
 183   }
 184   else {
 185     /* This block is never reached. */
 186
 187     int term = NIED2_SUB(0, ph[kj]);
 188
 189     for(r=kj+1; r<bigm; r++) {
 190       v[r] = arbitrary_element;
 191
 192       /* Check the condition of section 3.3,
 193        * remembering that the H's have the opposite sign.  [????????]
 194        */
 195       term = NIED2_SUB(term, NIED2_MUL(ph[r], v[r]));
 196     }
 197
 198     /* Now v[bigm] != term. */
 199     v[bigm] = NIED2_ADD(nonzero_element, term);
 200
 201     for(r=bigm+1; r<m; r++) {
 202       v[r] = arbitrary_element;
 203     }
 204   }
 205
 206   /* Calculate the remaining V's using the recursion of section 2.3,
 207    * remembering that the B's have the opposite sign.
 208    */
 209   for(r=0; r<=maxv-m; r++) {
 210     int term = 0;
 211     for(k=0; k<m; k++) {
 212       term = NIED2_SUB(term, NIED2_MUL(pb[k], v[r+k]));
 213     }
 214     v[r+m] = term;
 215   }
 216 }
 217
 218
 219 static void calculate_cj(nied2_state_t * ns, unsigned int dimension)
 220 {
 221   int ci[NIED2_NBITS][NIED2_NBITS];
 222   int v[MAXV+1];
 223   int r;
 224   unsigned int i_dim;
 225
 226   for(i_dim=0; i_dim<dimension; i_dim++) {
 227
 228     const int poly_index = i_dim + 1;
 229     int j, k;
 230
 231     /* Niederreiter (page 56, after equation (7), defines two
 232      * variables Q and U.  We do not need Q explicitly, but we
 233      * do need U.
 234      */
 235     int u = 0;
 236
 237     /* For each dimension, we need to calculate powers of an
 238      * appropriate irreducible polynomial, see Niederreiter
 239      * page 65, just below equation (19).
 240      * Copy the appropriate irreducible polynomial into PX,
 241      * and its degree into E.  Set polynomial B = PX ** 0 = 1.
 242      * M is the degree of B.  Subsequently B will hold higher
 243      * powers of PX.
 244      */
 245     int pb[NIED2_MAX_DEGREE+1];
 246     int px[NIED2_MAX_DEGREE+1];
 247     int px_degree = poly_degree[poly_index];
 248     int pb_degree = 0;
 249
 250     for(k=0; k<=px_degree; k++) {
 251       px[k] = primitive_poly[poly_index][k];
 252       pb[k] = 0;
 253     }
 254
 255     for (;k<NIED2_MAX_DEGREE+1;k++) {
 256       px[k] = 0;
 257       pb[k] = 0;
 258     }
 259
 260     pb[0] = 1;
 261
 262     for(j=0; j<NIED2_NBITS; j++) {
 263
 264       /* If U = 0, we need to set B to the next power of PX
 265        * and recalculate V.
 266        */
 267       if(u == 0) calculate_v(px, px_degree, pb, &pb_degree, v, MAXV);
 268
 269       /* Now C is obtained from V.  Niederreiter
 270        * obtains A from V (page 65, near the bottom), and then gets
 271        * C from A (page 56, equation (7)).  However this can be done
 272        * in one step.  Here CI(J,R) corresponds to
 273        * Niederreiter's C(I,J,R).
 274        */
 275       for(r=0; r<NIED2_NBITS; r++) {
 276         ci[r][j] = v[r+u];
 277       }
 278
 279       /* Advance Niederreiter's state variables. */
 280       ++u;
 281       if(u == px_degree) u = 0;
 282     }
 283
 284     /* The array CI now holds the values of C(I,J,R) for this value
 285      * of I.  We pack them into array CJ so that CJ(I,R) holds all
 286      * the values of C(I,J,R) for J from 1 to NBITS.
 287      */
 288     for(r=0; r<NIED2_NBITS; r++) {
 289       int term = 0;
 290       for(j=0; j<NIED2_NBITS; j++) {
 291         term = 2*term + ci[r][j];
 292       }
 293       ns->cj[r][i_dim] = term;
 294     }
 295
 296   }
 297 }
 298
 299
 300 static int nied2_init(void * state, unsigned int dimension)
 301 {
 302   nied2_state_t * n_state = (nied2_state_t *) state;
 303   unsigned int i_dim;
 304
 305   if(dimension < 1 || dimension > NIED2_MAX_DIMENSION) return GSL_EINVAL;
 306
 307   calculate_cj(n_state, dimension);
 308
 309   for(i_dim=0; i_dim<dimension; i_dim++) n_state->nextq[i_dim] = 0;
 310   n_state->sequence_count = 0;
 311
 312   return GSL_SUCCESS;
 313 }
 314
 315
 316 static int nied2_get(void * state, unsigned int dimension, double * v)
 317 {
 318   static const double recip = 1.0/(double)(1U << NIED2_NBITS); /* 2^(-nbits) */
 319   nied2_state_t * n_state = (nied2_state_t *) state;
 320   int r;
 321   int c;
 322   unsigned int i_dim;
 323
 324   /* Load the result from the saved state. */
 325   for(i_dim=0; i_dim<dimension; i_dim++) {
 326     v[i_dim] = n_state->nextq[i_dim] * recip;
 327   }
 328
 329   /* Find the position of the least-significant zero in sequence_count.
 330    * This is the bit that changes in the Gray-code representation as
 331    * the count is advanced.
 332    */
 333   r = 0;
 334   c = n_state->sequence_count;
 335   while(1) {
 336     if((c % 2) == 1) {
 337       ++r;
 338       c /= 2;
 339     }
 340     else break;
 341   }
 342
 343   if(r >= NIED2_NBITS) return GSL_EFAILED; /* FIXME: better error code here */
 344
 345   /* Calculate the next state. */
 346   for(i_dim=0; i_dim<dimension; i_dim++) {
 347     n_state->nextq[i_dim] ^= n_state->cj[r][i_dim];
 348   }
 349
 350   n_state->sequence_count++;
 351
 352   return GSL_SUCCESS;
 353 }