Parallelized Sobol Quasi Random Number Generator using Pthreads

/*
********************************************************************
Quasi Random Number Generator - serial version
Abishek Ramdas
********************************************************************
Inputs :
1. Number of patterns (limit 64000) can be varied by varying the LIM paramter
2. Number of dimensions

Outputs:
All patterns are stored in file QRNG_Patterns_sngl_thrd.txt
the patterns are printed along the x axis and the dimensions are along y axis
********************************************************************
*/


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


#define LIM 64000 //limitng the number of patterns generated to max of 64000
#define MAXBIT 30
#define MAXDIM 6
#define DEBUG_INIT 0
#define DEBUG_SEED 0
#define DEBUG_SOBOL 0
#define DEBUG_SOBOL_E 0

float random_numbers[MAXDIM][LIM];//dynamic declaration of 2D array is not possible

int main(int argc, char *argv[])
{
  int loop,init=1;
  unsigned long in_G,in_temp,no_of_vectors;
  int dim;//number of values to be generated
  int *n = &dim;
  float x[6];
  char filename[]="QRNG_Patterns_sngl_thrd.txt";//file contains patterns for a single threaded execution
  FILE *fp;
  fp=fopen(filename,"w");

  //  no_of_vectors = 20;

  int j,k,l;
  unsigned long i,im,ipp;
  static float fac;
  //  static unsigned long in,ix[MAXDIM+1],*iu[MAXBIT+1];
  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("Enter the number of vectors required:");
  scanf("%ld",&no_of_vectors);
  printf("Enter the dimensions:");
  scanf("%d",n);
  printf("single threaded version\n");
  if(init==1)
    {
      //      printf("initializing... \n");
      for(k=1;k<=MAXDIM;k++) ix[k]=0;
      in=0;
      if (iv[1] != 1) return;
      fac=1.0/(1L << MAXBIT);
      if(DEBUG_INIT)
    {
      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_INIT)
        {
          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_INIT)
        {
          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_INIT)
        {
          printf("iu[%d][%d] = %ld\n",j,k,iu[j][k]);
          printf("(maxbit-%d) = %d\n",j,(MAXBIT-j));
        }
        }
      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;
        }
    }
    }
 
  //Generate SEED for a particular index of random number
  in = -1;//in - index of (n-1)random number
  im = in++;
  in_G=in;
  in_temp = in;
  in_temp >>= 1;
  in_G ^= in_temp;
  if(DEBUG_SEED)
    {
      printf("Random # index : %ld\n",in);
      printf("Random # index in Gray in_G : %ld\n",in_G);
      printf("Generating random numbers from index: %ld\n",in+1);
    }
  for(k=1;k<=IMIN(*n,MAXDIM);k++)
    {
      in_G=in;
      in_temp = in;
      in_temp >>= 1;
      in_G ^= in_temp;
      if(DEBUG_SEED)
    {
      printf("ix[%d] =",k);
    }

      for(loop=0;loop*MAXDIM<(MAXBIT*MAXDIM+1);loop++)
    {
      //if LSB in in_G =0 then dont xor else xor
      if((in_G&1))
        {
          ix[k]^=iv[k+(loop*MAXDIM)];
          in_G >>= 1;
           if(DEBUG_SEED)
        {
          printf("iv[%d+%d]^",k,loop*MAXDIM);
        }
        }
      else
        {
          in_G >>= 1;
        }
    }
      if(DEBUG_SEED)
    {
      printf("\nseed value for rand #%ld : ix[%d]=%ld",in,k,ix[k]);
      printf("\n");
    }
    }

  //SEED is generated
  //Generate the random numbers next

  for(i=0;i<no_of_vectors;i++)//i is the number of random numbers to be generated
    {
      im=in++;
      if(DEBUG_SOBOL_E)
    {
      printf("in=%ld\n",in);
      printf("im=%ld\n",im);
    }
      for(j=1;j<=MAXBIT;j++)
    {
      if (!(im & 1)) break; //identify the location of the right most 0 in im. it is indicated by j
      im >>= 1;
      if(DEBUG_SOBOL_E)
        {
          printf("im2=%ld\n",im);
        }
    }
      if (j > MAXBIT) error(1,errno,"MAXBIT too small in sobseq");
      im=(j-1)*MAXDIM;
      if(DEBUG_SOBOL_E)
    {
      printf("j : %d\n",j);
      printf("im3=%ld\n",im);
      printf("IMIN %d\n",IMIN(*n,MAXDIM));
    }
     
      for(k=1;k<=IMIN(*n,MAXDIM);k++)
    {
      if(DEBUG_SOBOL)
        printf("ix[%d] = %ld\n",k,ix[k]);

      ix[k] ^= iv[im+k];
      x[k]=ix[k]*fac;
      random_numbers[k-1][i]=x[k];
      if(DEBUG_SOBOL)
        {
          printf("ix[%d]^= iv[%ld] => %ld ^= %ld => %f\t",k,im+k,ix[k],iv[im+k],x[k]);
      //          printf("fac : %1.30f\n",fac);
      //          printf("iv[%ld] : %ld\n",im+k,iv[im+k]);
        }

      if(!DEBUG_SOBOL)
        {
          printf("%f  ",x[k]);
        }
    }
      printf("\n");
     
    }
    printf("\nFinal Result\n");
  for(k=0;k<IMIN(*n,MAXDIM);k++)
    {
      for(i=0;i<no_of_vectors;i++)
    {
      //  printf("%f  ",random_numbers[k][i]);
      fprintf(fp,"%f  ",random_numbers[k][i]);//%f2spaces
    }
      fprintf(fp,"\n");
      //      printf("\n");
    }
  fclose(fp);
  printf("Random Numbers are stored in ./%s\n",filename);
  return 0;
}