/* This program is directly translated from ???'s program written in Fortran
   By Yaohang Li */
#include <stdio.h>

#define NIEDMAXDIM 100
#define NIEDNBITS  31
#define NIEDMAXQ   50

/*
  The following definitions concern the representation of
  polynomials.
*/
#define NIEDMAXDEG 50
#define NIEDDEG    10

/*
  In order to adjust the limit of C array, we change the
  definition of the polynomial array.
  In fortran code, POLY(-1) stores the degree. In our C code
  POLY(MAXDEG+1) stores the degree.
 
                    Yaohang Li Nov 2000

  The parameter MAXDEG gives the highest degree of polynomial
  to be handled.  Polynomials stored as arrays have the
  coefficient of degree n in POLY(N), and the degree of the
  polynomial in POLY().  The parameter DEG is just to remind
  us of this last fact.  A polynomial which is identically 0
  is given degree -1.

  A polynomial can also be stored in an integer I, with
        I = AN*Q**N + ... + A0.
  Routines ITOP and PTOI convert between these two formats.

  MAXE   - We need MAXDIM irreducible polynomials over Z2.
            MAXE is the highest degree among these.
  MAXV   - The maximum possible index used in V.
*/
#define NIEDMAXE 9
#define NIEDMAXV NIEDNBITS+NIEDMAXE

/*
     CJ    - Contains the packed values of Niederreiter's C(I,J,R)
     DIMEN   - The dimension of the sequence to be generated
     COUNT - The index of the current item in the sequence,
             expressed as an array of bits.  COUNT(R) is the same
             as Niederreiter's AR(N) (page 54) except that N is
             implicit.
     NEXTQ - The numerators of the next item in the series.  These
             are like Niederreiter's XI(N) (page 54) except that
             N is implicit, and the NEXTQ are integers.  To obtain
             the values of XI(N), multiply by RECIP (see GOLO2).

   Array CJ of the common block is set up by subroutine CALCC2.
   The other items in the common block are set up by INLO2.
*/
  int niedcj[NIEDMAXDIM][NIEDNBITS];
  int niedcount;
  int niednextq[NIEDMAXDIM];

/*
  The following COMMON block, used by many subroutines,
  gives the order Q of a field, its characteristic P, and its
  addition, multiplication and subtraction tables.
  The parameter MAXQ gives the order of the largest field to
  be handled.
*/

  int niedp, niedq, niedadd[NIEDMAXQ][NIEDMAXQ];
  int niedmul[NIEDMAXQ][NIEDMAXQ], niedsub[NIEDMAXQ][NIEDMAXQ];

int niedmax(int, int);
int niedmin(int, int);
int charac(int);
void calcc2(int);
void calcv(int *, int, int *, int *, int *, int);
void plymul(int *, int, int *, int, int *, int *);
void setfld(int);

int nied(int dimen, double *quasi)
{
/*
  The parameter RECIP is the multiplier which changes the
  integers in NEXTQ into the required real values in QUASI.
*/
  double recip;
  int i,r;

  recip=0.5/(double)(1<<NIEDNBITS-1);

/* 
  Multiply the numerators in NEXTQ by RECIP to get the next
  quasi-random vector
*/
  for (i=0;i<dimen;i++)
    quasi[i]=(double)niednextq[i]*recip;

/*
  Find the position of the right-hand zero in COUNT.  This
  is the bit that changes in the Gray-code representation as
  we go from COUNT to COUNT+1.
*/
  r=0;
  i=niedcount;
  while (i&1!=0)
    {
      r++;
      i>>=1;
    }

/*
  Check that we have not passed 2**NBITS calls on GOLO2
*/
  if (r>=NIEDNBITS)
    {
      fprintf(stderr, "Too many calls\n");
      return(-1);
    }

/*
  Compute the new numerators in vector NEXTQ
*/
  for (i=0;i<dimen;i++)
    niednextq[i]^=niedcj[i][r];

  niedcount++;
}

/*
  This subroutine calculates the values of Niederreiter's
  C(I,J,R) and performs other initialisation necessary
  before calling GOLO2.

  INPUT :
   DIMEN - The dimension of the sequence to be generated.
        {DIMEN is called DIM in the argument of INLO2,
        because DIMEN is subsequently passed via COMMON
        and is called DIMEN there.}

   SKIP  - The number of values to throw away at the beginning
           of the sequence.
*/
int init_nied(int dimen, int skip)
{
  int i, r, gray;

/*
  This assignment just relabels the variable for
  subsequent use.
*/
  if (dimen<0||dimen>NIEDMAXDIM)
    {
      fprintf(stderr,"Error: init_nied() bad dimension\n");
      return(-1);
    }

  calcc2(dimen);

/*
  Translate SKIP into Gray code
*/
  gray=skip^(skip/2);

/*
  Now set up NEXTQ appropriately for this value of the Gray code
*/
  for (i=0;i<dimen;i++)
    niednextq[i]=0;

  r=0;
  while (gray!=0)
    {
      if (gray&1!=0)
        for (i=0;i<dimen;i++)
          niednextq[i]^=niedcj[i][r];
      gray/=2;
      r++;
    }

  niedcount=skip;
  return(0);
}

/*
  This version :  12 February 1992

     *****  For base-2 only.

  See the general comments on implementing Niederreiter's
  low-discrepancy sequences.

  This program calculates the values of the constants C(I,J,R).
  As far as possible, we use Niederreiter's notation.
  For each value of I, we first calculate all the corresponding
  values of C :  these are held in the array CI.  All these
  values are either 0 or 1.  Next we pack the values into the
  array CJ, in such a way that CJ(I,R) holds the values of C
  for the indicated values of I and R and for every value of
  J from 1 to NBITS.  The most significant bit of CJ(I,R)
  (not counting the sign bit) is C(I,1,R) and the least
  significant bit is C(I,NBITS,R).
    When all the values of CJ have been calculated, we return
  this array to the calling program.

 INPUT :
   The array CJ to be initialised, and DIMEN the number of
   dimensions we are using, are transmitted through /COMM2/.

 OUTPUT :
   The array CJ is returned to the calling program.

 USES :
   Subroutine SETFLD is used to set up field arithmetic tables.
   (Although this is a little heavy-handed for the field of
   order 2, it makes for uniformity with the general program.)
   Subroutine CALCV is used to for the auxiliary calculation
   of the values V needed to get the Cs.
*/
void calcc2(int dimen)
{
/*
   This DATA statement supplies the coefficients and the
   degrees of the first 100 irreducible polynomials over Z2.
   They are taken from Lidl and Niederreiter, FINITE FIELDS,
   Cambridge University Press (1984), page 553.
   The order of these polynomials is the same as the order in
   file 'irrtabs.dat' used by the general program.

   In this block PX(I, NIEDDEG) is the degree of the Ith polynomial,
   and PX(I, N) is the coefficient of x**n in the Ith polynomial.
*/

  int irreddeg[NIEDMAXDIM]=
    {
      1, 1, 2, 3, 3, 4, 4, 4, 5, 5, 5, 5, 5, 5,
      6, 6, 6, 6, 6, 6, 6, 6, 6,
      7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
      8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
      8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
      9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9,
      9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9
    };
  int irred[NIEDMAXDIM][NIEDDEG]=
      {
         { 0,1,0,0,0,0,0,0,0,0 },
         { 1,1,0,0,0,0,0,0,0,0 },
         { 1,1,1,0,0,0,0,0,0,0 },
         { 1,1,0,1,0,0,0,0,0,0 },
         { 1,0,1,1,0,0,0,0,0,0 },
         { 1,1,0,0,1,0,0,0,0,0 },
         { 1,0,0,1,1,0,0,0,0,0 },
         { 1,1,1,1,1,0,0,0,0,0 },
         { 1,0,1,0,0,1,0,0,0,0 },
         { 1,0,0,1,0,1,0,0,0,0 },
         { 1,1,1,1,0,1,0,0,0,0 },
         { 1,1,1,0,1,1,0,0,0,0 },
         { 1,1,0,1,1,1,0,0,0,0 },
         { 1,0,1,1,1,1,0,0,0,0 },
         { 1,1,0,0,0,0,1,0,0,0 },
         { 1,0,0,1,0,0,1,0,0,0 }, 
         { 1,1,1,0,1,0,1,0,0,0 },
         { 1,1,0,1,1,0,1,0,0,0 },
         { 1,0,0,0,0,1,1,0,0,0 },
         { 1,1,1,0,0,1,1,0,0,0 }, 
         { 1,0,1,1,0,1,1,0,0,0 },
         { 1,1,0,0,1,1,1,0,0,0 },
         { 1,0,1,0,1,1,1,0,0,0 },
         { 1,1,0,0,0,0,0,1,0,0 }, 
         { 1,0,0,1,0,0,0,1,0,0 },
         { 1,1,1,1,0,0,0,1,0,0 },
         { 1,0,0,0,1,0,0,1,0,0 },
         { 1,0,1,1,1,0,0,1,0,0 },
         { 1,1,1,0,0,1,0,1,0,0 },
         { 1,1,0,1,0,1,0,1,0,0 }, 
         { 1,0,0,1,1,1,0,1,0,0 },
         { 1,1,1,1,1,1,0,1,0,0 },
         { 1,0,0,0,0,0,1,1,0,0 }, 
         { 1,1,0,1,0,0,1,1,0,0 },
         { 1,1,0,0,1,0,1,1,0,0 },
         { 1,0,1,0,1,0,1,1,0,0 },
         { 1,0,1,0,0,1,1,1,0,0 },
         { 1,1,1,1,0,1,1,1,0,0 },
         { 1,0,0,0,1,1,1,1,0,0 },
         { 1,1,1,0,1,1,1,1,0,0 },
         { 1,0,1,1,1,1,1,1,0,0 },
         { 1,1,0,1,1,0,0,0,1,0 },
         { 1,0,1,1,1,0,0,0,1,0 },
         { 1,1,0,1,0,1,0,0,1,0 },
         { 1,0,1,1,0,1,0,0,1,0 },
         { 1,0,0,1,1,1,0,0,1,0 },
         { 1,1,1,1,1,1,0,0,1,0 },
         { 1,0,1,1,0,0,1,0,1,0 },
         { 1,1,1,1,1,0,1,0,1,0 },
         { 1,1,0,0,0,1,1,0,1,0 },
         { 1,0,1,0,0,1,1,0,1,0 },
         { 1,0,0,1,0,1,1,0,1,0 },
         { 1,0,0,0,1,1,1,0,1,0 },
         { 1,1,1,0,1,1,1,0,1,0 },
         { 1,1,0,1,1,1,1,0,1,0 },
         { 1,1,1,0,0,0,0,1,1,0 },
         { 1,1,0,1,0,0,0,1,1,0 },
         { 1,0,1,1,0,0,0,1,1,0 },
         { 1,1,1,1,1,0,0,1,1,0 },
         { 1,1,0,0,0,1,0,1,1,0 },
         { 1,0,0,1,0,1,0,1,1,0 },
         { 1,0,0,0,1,1,0,1,1,0 },
         { 1,0,1,1,1,1,0,1,1,0 },
         { 1,1,0,0,0,0,1,1,1,0 },
         { 1,1,1,1,0,0,1,1,1,0 },
         { 1,1,1,0,1,0,1,1,1,0 },
         { 1,0,1,1,1,0,1,1,1,0 },
         { 1,1,1,0,0,1,1,1,1,0 },
         { 1,1,0,0,1,1,1,1,1,0 },
         { 1,0,1,0,1,1,1,1,1,0 },
         { 1,0,0,1,1,1,1,1,1,0 },
         { 1,1,0,0,0,0,0,0,0,1 },
         { 1,0,0,0,1,0,0,0,0,1 },
         { 1,1,1,0,1,0,0,0,0,1 },
         { 1,1,0,1,1,0,0,0,0,1 },
         { 1,0,0,0,0,1,0,0,0,1 },
         { 1,0,1,1,0,1,0,0,0,1 },
         { 1,1,0,0,1,1,0,0,0,1 },
         { 1,1,0,1,0,0,1,0,0,1 },
         { 1,0,0,1,1,0,1,0,0,1 },
         { 1,1,1,1,1,0,1,0,0,1 },
         { 1,0,1,0,0,1,1,0,0,1 },
         { 1,0,0,1,0,1,1,0,0,1 },
         { 1,1,1,1,0,1,1,0,0,1 },
         { 1,1,1,0,1,1,1,0,0,1 },
         { 1,0,1,1,1,1,1,0,0,1 }, 
         { 1,1,1,0,0,0,0,1,0,1 },
         { 1,0,1,0,1,0,0,1,0,1 },
         { 1,0,0,1,1,0,0,1,0,1 },
         { 1,1,0,0,0,1,0,1,0,1 },
         { 1,0,1,0,0,1,0,1,0,1 },
         { 1,1,1,1,0,1,0,1,0,1 },
         { 1,1,1,0,1,1,0,1,0,1 },
         { 1,0,1,1,1,1,0,1,0,1 },
         { 1,1,1,1,0,0,1,1,0,1 },
         { 1,0,0,0,1,0,1,1,0,1 },
         { 1,1,0,1,1,0,1,1,0,1 },
         { 1,0,1,0,1,1,1,1,0,1 },
         { 1,0,0,1,1,1,1,1,0,1 },
         { 1,0,0,0,0,0,0,0,1,1 }
    };

  int px[NIEDMAXDEG+1], b[NIEDMAXDEG+1];
  int v[NIEDMAXV+1], ci[NIEDNBITS+1][NIEDNBITS];
  int pxdeg, bdeg; /* the degree of px and b */
  int e, i, j, r, u, term;

  /* prepare to work in Z2 */
  setfld(2);

/*
  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.
*/
  for (i=0;i<dimen;i++)
    {
      e=irreddeg[i];
      for (j=0;j<=e;j++)
        px[j]=irred[i][j];
      pxdeg=irreddeg[i];
      bdeg=0;
      b[0]=1;

      /*
        Niederreiter (page 56, after equation (7), defines two
        variables Q and U.  We do not need Q explicitly, but we
        do need U.
      */
      u=0;

      for (j=1;j<=NIEDNBITS;j++)
        {
          /*
            If U = 0, we need to set B to the next power of PX
            and recalculate V.  This is done by subroutine CALCV.
          */
          if (u==0)
            calcv(px,pxdeg,b,&bdeg,v,NIEDMAXV);

          /*
            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<NIEDNBITS;r++)
            ci[j][r]=v[r+u];
          u++;
          if (u==e)
            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<NIEDNBITS;r++)
        {
          term=0;
          for (j=1;j<=NIEDNBITS;j++)
            term=term*2+ci[j][r];
          niedcj[i][r]=term;
        }
    }
} 

/*
   This version :  12 February 1991

   See the general comments on implementing Niederreiter's
   low-discrepancy sequences.

   This program calculates the values of the constants V(J,R) as
   described in BFN section 3.3.  It is called from either CALCC or
   CALCC2.  The values transmitted through common /FIELD/ determine
   which field we are working in.

 INPUT :
   PX is the appropriate irreducible polynomial for the dimension
     currently being considered.  The degree of PX will be called E.
   B is the polynomial defined in section 2.3 of BFN.  On entry to
     the subroutine, the degree of B implicitly defines the parameter
     J of section 3.3, by degree(B) = E*(J-1).
   MAXV gives the dimension of the array V.
   On entry, we suppose that the common block /FIELD/ has been set
     up correctly (using SETFLD).

 OUTPUT :
   On output B has been multiplied by PX, so its degree is now E*J.
   V contains the values required.

 USES :
   The subroutine PLYMUL is used to multiply polynomials.
*/
void calcv(int *px, int pxdeg, int *b, int *bdeg, int *v, int maxv)
{
  int e, i, j, kj, m, bigm, r, term;
  int h[NIEDMAXDEG+1];
  int hdeg; /* degree for h */

  int arbit, nonzer;

/*
   We use ARBIT() to indicate where the user can place
   an arbitrary element of the field of order Q, while NONZER
   shows where he must put an arbitrary non-zero element
   of the same field.  For the code,
   this means 0 <= ARBIT < Q and 0 < NONZER < Q.  Within these
   limits, the user can do what he likes.  ARBIT is declared as
   a function as a reminder that a different arbitrary value may
   be returned each time ARBIT is referenced.

    BIGM is the M used in section 3.3.
    It differs from the [little] m used in section 2.3,
    denoted here by M.
*/
  arbit=1;

  nonzer=1;

  e=pxdeg;

/*
   The poly H 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 !
*/

  for (i=0;i<=*bdeg;i++)
    h[i]=b[i];
  hdeg=*bdeg;
  bigm=hdeg;

/*
   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 !
*/

  plymul(px,pxdeg,b,*bdeg,b,bdeg);
  m=*bdeg;

/* We don't use j explicitly anywhare, but here it is just in case */
  j=m/e;

/*
  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.
*/
  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)
    {
      term=niedsub[0][h[kj]];
      for (r=kj+1;r<bigm;r++)
        {
          v[r]=arbit;
          /*
            Check the condition of section 3.3,
            remembering that the H's have the opposite sign.
          */
          term=niedsub[term][niedmul[h[r]][v[r]]];
        }
      /* Now V(BIGM) is anything but TERM */
      v[bigm]=niedadd[nonzer][term];
      for (r=bigm+1;r<m;r++)
        v[r]=arbit;
    }
  else
    /* this is the case kj greater than bigm */
    for (r=kj+1;r<m;r++)
      v[r]=arbit;

/*
   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++)
    {
      term=0;
      for (i=0;i<m;i++)
        term=niedsub[term][niedmul[b[i]][v[r+i]]];
      v[r+m]=term;
    }
}

/*
   Multiplies polynomial PA by PB putting the result in PC.
   Coefficients are elements of the field of order Q.
*/
void plymul(int *pa, int padeg, int *pb, int pbdeg, int *pc, int *pcdeg)
{
  int i, j, dega, degb, degc, term;
  int pt[NIEDMAXDEG+1];

  dega=padeg;
  degb=pbdeg;

  if (dega==-1||degb==-1)
    degc=-1;
  else
    degc=dega+degb;

  if (degc>NIEDMAXDEG)
    {
      fprintf(stderr,"Error: PLYMUL: Degree of product exceeds NIEDMAXDEG\n");
      exit(-1);
    }

  for (i=0;i<=degc;i++)
    {
      term=0;
      for (j=niedmax(0,i-dega);j<=niedmin(degb,i);j++)
        term=niedadd[term][niedmul[pa[i-j]][pb[j]]];
      pt[i]=term;
    }

  *pcdeg=degc;
  for (i=0;i<=degc;i++)
    pc[i]=pt[i];
  for (i=degc+1;i<=NIEDMAXDEG;i++)
    pc[i]=0;
}

/*
   This version : 12 December 1991

   This subroutine sets up addition, multiplication, and
   subtraction tables for the finite field of order QIN.
   If necessary, it reads precalculated tables from the file
   'gftabs.dat' using unit 1.  These precalculated tables
   are supposed to have been put there by GFARIT.

      *****  For the base-2 programs, these precalculated
      *****  tables are not needed and, therefore, neither
      *****  is GFARIT.


   Unit 1 is closed both before and after the call of SETFLD.

 USES
   Integer function CHARAC gets the characteristic of a field.
*/
void setfld(int qin)
{
  int i, j, n;

  if (qin<1||qin>NIEDMAXQ)
    {
      fprintf(stderr,"Error: setfld: Bad Value of Q\n");
      exit(-1);   
    }

  niedq=qin;
  niedp=charac(niedq-1);
 
  if (niedp==0)
    {
      fprintf(stderr,"Error: setfld: There is no field of order %d\n", niedq);
      exit(-1);
    }

  /* set up to handle a field of prime order: calculate ADD and MUL */
  if (niedp==niedq)
    for (i=0;i<niedq;i++)
      for (j=0;j<niedq;j++)
        {
          niedadd[i][j]=(i+j)%niedp;
          niedmul[i][j]=(i*j)%niedp;
        }
  else
    {
      /*
        Set up to handle a field of prime-power order :  tables for
        ADD and MUL are in the file 'gftabs.dat'.
        We are currently lack of this table
                          Yaohang Li
      */
      fprintf(stderr,"Error: setfld: Need the table gftabs.dat\n");
      exit(-1);
    } 

  /* Now use the addition table to set the subtraction table */
  for (i=0;i<niedq;i++)
    for (j=0;j<niedq;j++)
      niedsub[niedadd[i][j]][i]=j;
}

/*
   This version :  12 December 1991

   This function gives the characteristic for a field of
   order QIN.  If no such field exists, or if QIN is out of
   the range we can handle, returns 0.
*/
int charac(int qin)
{
  int ch[NIEDMAXQ]=
    {
       0,  2,  3,  2,  5,  0,  7,  2,  3,  0,
      11,  0, 13,  0,  0,  2, 17,  0, 19,  0,
       0,  0, 23,  0,  5,  0,  3,  0, 29,  0,
      31,  2,  0,  0,  0,  0, 37,  0,  0,  0,
      41,  0, 43,  0,  0,  0, 47,  0,  7,  0
    };

  if (qin<1||qin>NIEDMAXQ)
    return(0);
  else
    return(ch[qin]);
}

int niedmax(int x, int y)
{
  if (x>y)
    return(x);
  else
    return(y);
}

int niedmin(int x, int y)
{
  if (x<y)
    return(x);
  else
    return(y);
}

int main()
{
  int i;

  double quasi[NIEDMAXDIM];

  init_nied(5,0);
  for (i=0;i<1000;i++)
    {
      nied(NIEDMAXDIM,quasi);
      printf("%f, %f, %f, %f, %f\n",quasi[0], quasi[1], quasi[2], quasi[3], quasi[4]);
    }
  
  return(0);
}