/* Author: G. Jungman

 */

/* Implement Niederreiter base 2 generator.

 * See:

 *   Bratley, Fox, Niederreiter, ACM Trans. Model. Comp. Sim. 2, 195 (1992)

 */

#include <config.h>

#include <gsl/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

    v[r] = 0;

  }

  v[kj] = 1;

  if(kj >= bigm) {

    for(r=kj+1; r

      v[r] = arbitrary_element;

    }

  }

  else {

    /* This block is never reached. */

    int term = NIED2_SUB(0, ph[kj]);

    for(r=kj+1; 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

      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

      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

    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

      px[k] = 0;

      pb[k] = 0;

    }

    pb[0] = 1;

    for(j=0; 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

        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

      int term = 0;

      for(j=0; 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

    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

    n_state->nextq[i_dim] ^= n_state->cj[r][i_dim];

  }

  n_state->sequence_count++;

  return GSL_SUCCESS;

}