Sobol Quasi Random Number Generator single thread version

#include<stdio.h>
#include<stdlib.h>
#include<error.h>
#include<errno.h>
#include"nrutil.h"

#define MAXBIT 30
#define MAXDIM 6
#define DEBUG 1

void sobseq(int *, float []);

int main(int argc, char *argv[])
{
  int n,i;
  float x[100];
  for(i=0;i<100;i++)
    x[i]=0;
  n=-1;
  sobseq(&n,x);
  //  for(i=0;i<100;i++)
  //  printf("%f",x[i]);
  n=4;
  sobseq(&n,x);
  for(i=0;i<100;i++)
    printf("%f\n",x[i]);
 n=5;
  sobseq(&n,x);
  for(i=0;i<100;i++)
    printf("%f\n",x[i]);
  return 0;
}

void sobseq(int *n,float x[])
{
  int j,k,l;
  unsigned long i,im,ipp;
  static float fac;
  static unsigned long in,ix[MAXDIM+1],*iu[MAXBIT+1];
  static unsigned long mdeg[MAXDIM+1]={0,1,2,3,3,4,4};
  static unsigned long ip[MAXDIM+1]={0,0,1,1,2,1,4};
  static unsigned long iv[MAXDIM*MAXBIT+1]={0,1,1,1,1,1,1,3,1,3,3,1,1,5,7,7,3,3,5,15,11,5,15,13,9};
  printf("n = %d\n",*n);
  if(*n<0)
    {
      for(k=1;k<=MAXDIM;k++) ix[k]=0;
      in=0;
      if (iv[1] != 1) return;
      fac=1.0/(1L << MAXBIT);
      if(DEBUG)
    {
      printf("debug:\n");
      printf("n: %d\n",*n);
      printf("fac: %1.30f\n",fac);
      printf("1L: %ld\n",1L);
      printf("end debug\n");
    }
      for (j=1,k=0;j<=MAXBIT;j++,k+=MAXDIM)
    {
      //iu is array of arrays. contains the address of each element of array iv
      //initialization iu so that addresses of iv are loaded in to iu
      iu[j] = &iv[k];
      if(DEBUG)
        {
          printf("j= %d,k = %d:\t%ld\n",j,k,*iu[j]);
        }
    }
      for (k=1;k<=MAXDIM;k++)
    {
      for (j=1;j<=mdeg[k];j++)
        {
          if(DEBUG)
        {
          printf("mdeg = %ld\n",mdeg[k]);
          printf("iu[%d][%d] = %ld\n",j,k,iu[j][k]);
          printf("(maxbit-%d) = %d\n",j,(MAXBIT-j));
        }
          //initialize vi values
          iu[j][k] <<= (MAXBIT-j);
          if(DEBUG)
        {
          printf("iu[%d][%d] = %ld\n",j,k,iu[j][k]);
          printf("(maxbit-%d) = %d\n",j,(MAXBIT-j));
        }
          printf("\n");
        }
      for(j=mdeg[k]+1;j<=MAXBIT;j++)
        {
          ipp=ip[k];
          i=iu[j-mdeg[k]][k];
          i ^= (i >> mdeg[k]);
          for (l=mdeg[k]-1;l>=1;l--)
        {
          if (ipp & 1) i ^= iu[j-l][k];
          ipp >>= 1;
        }
          iu[j][k]=i;
          printf("iu[%d][%d]: %f\n",j,k,iu[j][k]*fac);
        }
    }
    }
  else
    {
      printf("inside else\n ");
      im=in++;
      for(j=1;j<=MAXBIT;j++)
    {
      if (!(im & 1)) break;
      im >>= 1;
    }
      if (j > MAXBIT) error(1,errno,"MAXBIT too small in sobseq");
      im=(j-1)*MAXDIM;
      printf("IMIN %d\n",IMIN(*n,MAXDIM));
      for(k=1;k<=IMIN(*n,MAXDIM);k++)
    {
      ix[k] ^= iv[im+k];
      x[k]=ix[k]*fac;
      if(DEBUG)
        {
          printf("x[%d] : %f\n",k,x[k]);
          printf("fac : %1.30f\n",fac);
          printf("iv[%ld] : %ld\n",im+k,iv[im+k]);
        }
    }
    }
}