123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354 |
- /* Author: G. Jungman
- */
- /* Implement Niederreiter base 2 generator.
- * See:
- * Bratley, Fox, Niederreiter, ACM Trans. Model. Comp. Sim. 2, 195 (1992)
- */
- #include "gsl__config.h"
- #include "gsl_qrng.h"
- #define NIED2_CHARACTERISTIC 2
- #define NIED2_BASE 2
- #define NIED2_MAX_DIMENSION 12
- #define NIED2_MAX_PRIM_DEGREE 5
- #define NIED2_MAX_DEGREE 50
- #define NIED2_BIT_COUNT 30
- #define NIED2_NBITS (NIED2_BIT_COUNT+1)
- #define MAXV (NIED2_NBITS + NIED2_MAX_DEGREE)
- /* Z_2 field operations */
- #define NIED2_ADD(x,y) (((x)+(y))%2)
- #define NIED2_MUL(x,y) (((x)*(y))%2)
- #define NIED2_SUB(x,y) NIED2_ADD((x),(y))
- static size_t nied2_state_size(unsigned int dimension);
- static int nied2_init(void * state, unsigned int dimension);
- static int nied2_get(void * state, unsigned int dimension, double * v);
- static const gsl_qrng_type nied2_type =
- {
- "niederreiter-base-2",
- NIED2_MAX_DIMENSION,
- nied2_state_size,
- nied2_init,
- nied2_get
- };
- const gsl_qrng_type * gsl_qrng_niederreiter_2 = &nied2_type;
- /* primitive polynomials in binary encoding */
- static const int primitive_poly[NIED2_MAX_DIMENSION+1][NIED2_MAX_PRIM_DEGREE+1] =
- {
- { 1, 0, 0, 0, 0, 0 }, /* 1 */
- { 0, 1, 0, 0, 0, 0 }, /* x */
- { 1, 1, 0, 0, 0, 0 }, /* 1 + x */
- { 1, 1, 1, 0, 0, 0 }, /* 1 + x + x^2 */
- { 1, 1, 0, 1, 0, 0 }, /* 1 + x + x^3 */
- { 1, 0, 1, 1, 0, 0 }, /* 1 + x^2 + x^3 */
- { 1, 1, 0, 0, 1, 0 }, /* 1 + x + x^4 */
- { 1, 0, 0, 1, 1, 0 }, /* 1 + x^3 + x^4 */
- { 1, 1, 1, 1, 1, 0 }, /* 1 + x + x^2 + x^3 + x^4 */
- { 1, 0, 1, 0, 0, 1 }, /* 1 + x^2 + x^5 */
- { 1, 0, 0, 1, 0, 1 }, /* 1 + x^3 + x^5 */
- { 1, 1, 1, 1, 0, 1 }, /* 1 + x + x^2 + x^3 + x^5 */
- { 1, 1, 1, 0, 1, 1 } /* 1 + x + x^2 + x^4 + x^5 */
- };
- /* degrees of primitive polynomials */
- static const int poly_degree[NIED2_MAX_DIMENSION+1] =
- {
- 0, 1, 1, 2, 3, 3, 4, 4, 4, 5, 5, 5, 5
- };
- typedef struct
- {
- unsigned int sequence_count;
- int cj[NIED2_NBITS][NIED2_MAX_DIMENSION];
- int nextq[NIED2_MAX_DIMENSION];
- } nied2_state_t;
- static size_t nied2_state_size(unsigned int dimension)
- {
- return sizeof(nied2_state_t);
- }
- /* Multiply polynomials over Z_2.
- * Notice use of a temporary vector,
- * side-stepping aliasing issues when
- * one of inputs is the same as the output
- * [especially important in the original fortran version, I guess].
- */
- static void poly_multiply(
- const int pa[], int pa_degree,
- const int pb[], int pb_degree,
- int pc[], int * pc_degree
- )
- {
- int j, k;
- int pt[NIED2_MAX_DEGREE+1];
- int pt_degree = pa_degree + pb_degree;
- for(k=0; k<=pt_degree; k++) {
- int term = 0;
- for(j=0; j<=k; j++) {
- const int conv_term = NIED2_MUL(pa[k-j], pb[j]);
- term = NIED2_ADD(term, conv_term);
- }
- pt[k] = term;
- }
- for(k=0; k<=pt_degree; k++) {
- pc[k] = pt[k];
- }
- for(k=pt_degree+1; k<=NIED2_MAX_DEGREE; k++) {
- pc[k] = 0;
- }
- *pc_degree = pt_degree;
- }
- /* Calculate the values of the constants V(J,R) as
- * described in BFN section 3.3.
- *
- * px = appropriate irreducible polynomial for current dimension
- * pb = polynomial defined in section 2.3 of BFN.
- * pb is modified
- */
- static void calculate_v(
- const int px[], int px_degree,
- int pb[], int * pb_degree,
- int v[], int maxv
- )
- {
- const int nonzero_element = 1; /* nonzero element of Z_2 */
- const int arbitrary_element = 1; /* arbitray element of Z_2 */
- /* The polynomial ph is px**(J-1), which is the value of B on arrival.
- * In section 3.3, the values of Hi are defined with a minus sign:
- * don't forget this if you use them later !
- */
- int ph[NIED2_MAX_DEGREE+1];
- /* int ph_degree = *pb_degree; */
- int bigm = *pb_degree; /* m from section 3.3 */
- int m; /* m from section 2.3 */
- int r, k, kj;
- for(k=0; k<=NIED2_MAX_DEGREE; k++) {
- ph[k] = pb[k];
- }
- /* Now multiply B by PX so B becomes PX**J.
- * In section 2.3, the values of Bi are defined with a minus sign :
- * don't forget this if you use them later !
- */
- poly_multiply(px, px_degree, pb, *pb_degree, pb, pb_degree);
- m = *pb_degree;
- /* Now choose a value of Kj as defined in section 3.3.
- * We must have 0 <= Kj < E*J = M.
- * The limit condition on Kj does not seem very relevant
- * in this program.
- */
- /* Quoting from BFN: "Our program currently sets each K_q
- * equal to eq. This has the effect of setting all unrestricted
- * values of v to 1."
- * Actually, it sets them to the arbitrary chosen value.
- * Whatever.
- */
- kj = bigm;
- /* Now choose values of V in accordance with
- * the conditions in section 3.3.
- */
- for(r=0; r<kj; r++) {
- v[r] = 0;
- }
- v[kj] = 1;
- if(kj >= bigm) {
- for(r=kj+1; r<m; r++) {
- v[r] = arbitrary_element;
- }
- }
- else {
- /* This block is never reached. */
- int term = NIED2_SUB(0, ph[kj]);
- for(r=kj+1; r<bigm; r++) {
- v[r] = arbitrary_element;
- /* Check the condition of section 3.3,
- * remembering that the H's have the opposite sign. [????????]
- */
- term = NIED2_SUB(term, NIED2_MUL(ph[r], v[r]));
- }
- /* Now v[bigm] != term. */
- v[bigm] = NIED2_ADD(nonzero_element, term);
- for(r=bigm+1; r<m; r++) {
- v[r] = arbitrary_element;
- }
- }
- /* Calculate the remaining V's using the recursion of section 2.3,
- * remembering that the B's have the opposite sign.
- */
- for(r=0; r<=maxv-m; r++) {
- int term = 0;
- for(k=0; k<m; k++) {
- term = NIED2_SUB(term, NIED2_MUL(pb[k], v[r+k]));
- }
- v[r+m] = term;
- }
- }
- static void calculate_cj(nied2_state_t * ns, unsigned int dimension)
- {
- int ci[NIED2_NBITS][NIED2_NBITS];
- int v[MAXV+1];
- int r;
- unsigned int i_dim;
- for(i_dim=0; i_dim<dimension; i_dim++) {
- const int poly_index = i_dim + 1;
- int j, k;
- /* Niederreiter (page 56, after equation (7), defines two
- * variables Q and U. We do not need Q explicitly, but we
- * do need U.
- */
- int u = 0;
- /* For each dimension, we need to calculate powers of an
- * appropriate irreducible polynomial, see Niederreiter
- * page 65, just below equation (19).
- * Copy the appropriate irreducible polynomial into PX,
- * and its degree into E. Set polynomial B = PX ** 0 = 1.
- * M is the degree of B. Subsequently B will hold higher
- * powers of PX.
- */
- int pb[NIED2_MAX_DEGREE+1];
- int px[NIED2_MAX_DEGREE+1];
- int px_degree = poly_degree[poly_index];
- int pb_degree = 0;
- for(k=0; k<=px_degree; k++) {
- px[k] = primitive_poly[poly_index][k];
- pb[k] = 0;
- }
- for (;k<NIED2_MAX_DEGREE+1;k++) {
- px[k] = 0;
- pb[k] = 0;
- }
- pb[0] = 1;
- for(j=0; j<NIED2_NBITS; j++) {
- /* If U = 0, we need to set B to the next power of PX
- * and recalculate V.
- */
- if(u == 0) calculate_v(px, px_degree, pb, &pb_degree, v, MAXV);
- /* Now C is obtained from V. Niederreiter
- * obtains A from V (page 65, near the bottom), and then gets
- * C from A (page 56, equation (7)). However this can be done
- * in one step. Here CI(J,R) corresponds to
- * Niederreiter's C(I,J,R).
- */
- for(r=0; r<NIED2_NBITS; r++) {
- ci[r][j] = v[r+u];
- }
- /* Advance Niederreiter's state variables. */
- ++u;
- if(u == px_degree) u = 0;
- }
- /* The array CI now holds the values of C(I,J,R) for this value
- * of I. We pack them into array CJ so that CJ(I,R) holds all
- * the values of C(I,J,R) for J from 1 to NBITS.
- */
- for(r=0; r<NIED2_NBITS; r++) {
- int term = 0;
- for(j=0; j<NIED2_NBITS; j++) {
- term = 2*term + ci[r][j];
- }
- ns->cj[r][i_dim] = term;
- }
- }
- }
- static int nied2_init(void * state, unsigned int dimension)
- {
- nied2_state_t * n_state = (nied2_state_t *) state;
- unsigned int i_dim;
- if(dimension < 1 || dimension > NIED2_MAX_DIMENSION) return GSL_EINVAL;
- calculate_cj(n_state, dimension);
- for(i_dim=0; i_dim<dimension; i_dim++) n_state->nextq[i_dim] = 0;
- n_state->sequence_count = 0;
- return GSL_SUCCESS;
- }
- static int nied2_get(void * state, unsigned int dimension, double * v)
- {
- static const double recip = 1.0/(double)(1U << NIED2_NBITS); /* 2^(-nbits) */
- nied2_state_t * n_state = (nied2_state_t *) state;
- int r;
- int c;
- unsigned int i_dim;
- /* Load the result from the saved state. */
- for(i_dim=0; i_dim<dimension; i_dim++) {
- v[i_dim] = n_state->nextq[i_dim] * recip;
- }
- /* Find the position of the least-significant zero in sequence_count.
- * This is the bit that changes in the Gray-code representation as
- * the count is advanced.
- */
- r = 0;
- c = n_state->sequence_count;
- while(1) {
- if((c % 2) == 1) {
- ++r;
- c /= 2;
- }
- else break;
- }
- if(r >= NIED2_NBITS) return GSL_EFAILED; /* FIXME: better error code here */
- /* Calculate the next state. */
- for(i_dim=0; i_dim<dimension; i_dim++) {
- n_state->nextq[i_dim] ^= n_state->cj[r][i_dim];
- }
- n_state->sequence_count++;
- return GSL_SUCCESS;
- }
|