rand.c

来自「用于LDPC编码译码的仿真实现。包括随机生成校验矩阵、由校验矩阵产生生成矩阵、编」· C语言 代码 · 共 388 行

/* RAND.C - Random number generation module. */

/* Copyright (c) 1995, 1996, 2000, 2001 by Radford M. Neal 
 *
 * Permission is granted for anyone to copy, use, or modify this program 
 * for purposes of research or education, provided this copyright notice 
 * is retained, and note is made of any changes that have been made. 
 *
 * This program is distributed without any warranty, express or implied.
 * As this program was written for research purposes only, it has not been
 * tested to the degree that would be advisable in any important application.
 * All use of this program is entirely at the user's own risk.
 */


/* NOTE:  See rand.html for documentation on these procedures. */


#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

#include "rand.h"
#include "open.h"


/* This module uses the 'drand48' pseudo-random number generator found
   on most Unix systems, the output of which is combined with a file
   of real random numbers.

   Many of the methods used in this module may be found in the following
   reference:

      Devroye, L. (1986) Non-Uniform Random Variate Generation, 
        New York: Springer-Verlag.

   The methods used here are not necessarily the fastest available.  They're
   selected to be reasonably fast while also being easy to write.
*/


/* CONSTANT PI.  Defined here if not in <math.h>. */

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif


/* TABLES OF REAL RANDOM NUMBERS.  A file of 100000 real random numbers
   (NOT pseudo-random) is used in conjunction with pseudo-random numbers
   for extra insurance.  These are employed in the form of five tables
   of 5000 32-bit integers.  

   The file must be located at the path given by RAND_FILE, which should
   be defined on the "cc" command line. */

#define Table_size 5000			/* Number of words in each table */

static int rn[N_tables][Table_size];	/* Random number tables */


/* STATE OF RANDOM NUMBER GENERATOR. */

static int initialized = 0;		/* Has module been initialized? */

static rand_state state0;		/* Default state structure */

static rand_state *state;		/* Pointer to current state */


/* INITIALIZE MODULE.  Sets things up using the default state structure,
   set as if rand_seed had been called with a seed of one. */

static void initialize (void)
{
  int i, j, k, w;
  char b;
  FILE *f;

  if (!initialized)
  {
    f = fopen("RAND_FILE","rb");
    
    if (f==NULL)
    { fprintf(stderr,"Can't open file of random numbers (%s)\n","RAND_FILE");
      exit(1);
    }

    for (i = 0; i<N_tables; i++)
    { for (j = 0; j<Table_size; j++)
      { w = 0;
        for (k = 0; k<4; k++)
        { if (fread(&b,1,1,f)!=1)
          { fprintf(stderr,"Error reading file of random numbers (%s)\n",
                            "RAND_FILE");
            exit(1);
          }
          w = (w<<8) | (b&0xff);
        }
        rn[i][j] = w;
      }
    }

    state = &state0;

    initialized = 1;

    rand_seed(1);
  }
}


/* TABLES OF REAL RANDOM NUMBERS.  A file of 100000 real random numbers
   (NOT pseudo-random) is used in conjunction with pseudo-random numbers
   for extra insurance.  These are employed in the form of five tables
   of 5000 32-bit integers.  

   The file must be located at the path given by RAND_FILE, which should
   be defined on the "cc" command line. */



/* GENERATE RANDOM INTEGER FROM 0, 1, ..., (n-1).  by lin jia ru*/

int rand_int( int n, int seed)
{ 
    int m;

	nextm: if ((m=(int) (n * rand_uniopen(seed)))==n)
			goto nextm;

	return (m);
}


//not repeat random int number from 0 to N-1 by lin jia ru 
void rand_int_num( int * rand_col, int N, int seed)
{	
	int i,j;
	int temp2;
	int *num_temp;
	num_temp = (int *)malloc(N*sizeof(int));

	for (i=0; i<N; i++)
		num_temp[i] = i;		//initiation

	for (i=0; i<N-1; i++)
	{
		temp2 = rand_int(N-i, seed);			//random read pointer
		rand_col[i] = num_temp[temp2];
		for (j=temp2; j<N-i-1; j++) 
			num_temp[j] = num_temp[j+1];
	}
	rand_col[N-1] = num_temp[0];

	free( num_temp);
}


/* GENERATE UNIFORMLY FORM (0,1). by lin jia ru */

static double r[98];
static int iff,ix1,ix2,ix3;
double rand_uniopen (int seed)
//double ran1(int& idum)
{
    double rm1,rm2,t;
	int m1,m2,m3,ia1,ia2,ia3,ic1,ic2,ic3,j;
	 
    m1 = 259200, ia1 = 7141, ic1 = 54773, rm1 = 0.0000038580247,
    m2 = 134456, ia2 = 8121, ic2 = 28411, rm2 = 0.0000074373773,
    m3 = 243000, ia3= 4561, ic3 = 51349;
	
    if ((seed < 0) || (iff == 0))
	{
        iff = 1;
        ix1 =  (ic1 - seed) % m1;
        ix1 = (ia1 * ix1 + ic1) % m1;
        ix2 = ix1 % m2;
        ix1 = (ia1 * ix1 + ic1) % m1;
        ix3 = ix1 % m3;
        for (j = 1; j<=97; j++)
		{
            ix1 = (ia1 * ix1 + ic1)% m1;
            ix2 = (ia2 * ix2 + ic2)% m2;
            r[j] = ((double)(ix1) + (double)(ix2) * rm2) * rm1;
        }
//        seed = 1;
    seed = (ia1 * seed + ic2) % m3;
	}
      ix1 = (ia1 * ix1 + ic1) % m1;
      ix2 = (ia2 * ix2 + ic2) % m2;
      ix3 = (ia3 * ix3 + ic3) % m3;
      j = 1 + (int)((97 * ix3) / m3);
      if ((j > 97) || (j < 1)) 
	  {
		  printf("abnormal exit\n");
	      return 1;
	  }
      t = r[j];
      r[j] = ((float)(ix1) + (float)(ix2) * rm2) * rm1;
	  return t;
}

//generate l Length of block 0 or 1 random numbers and store in *data 
//return(1) if OK! by lin xuehong

void rand01
( int *data,    /* data to write */
  int l_n,       /* Length of block */
  int seed	//initiation number for random
)
{	
	int n;

  /* generate l Length of block 0 or 1 random numbers */
    

	for (n = 0; n<l_n; n++)
	{
		if (rand_uniform(seed)>0.5) data[n] = 1;
		else data[n] = 0;
	}


}
float uniform(int *seed)
{   /*Generate a random number between 0 and 1*/
	float t;
	*seed = 2045*(*seed)+1;
	*seed = *seed - (*seed/1048576)*1048576;
	t = (float)((*seed)/1048576.0);
	return t;
}

void gen_symbol(int *symbol, int symbollength,double p[],int size)
{
	double P,temp;
	int i,j;
	
	int s = rand();
	for(i=1; i<symbollength+1; i++)
	{
		temp = uniform(&s);
		P=0;
		for(j=size-1;j>=0;j--)
		{
			P=P+p[j];
			if ((temp>1-P) && (temp<1-P+p[j]))
			*(symbol+i)=j;
		}
	}
}

/*p=P(0|0)=P(1|1)*/
void markov_symbol(int *symbol, int symbollength,double p)
{
	double temp;
	int i;
	
	int s = rand();

	temp = uniform(&s);
	if (temp>0.5)
		symbol[0]=1;
	else
		symbol[0]=0;

	for(i=1; i<symbollength; i++)
	{
		temp = uniform(&s);
		if (temp<p)
			symbol[i]=symbol[i-1];
		else
			symbol[i]=1^symbol[i-1];
			
	}
	
}
/* GENERATE UNIFORMLY FROM [0,1). by lin jia ru*/

static double r[98];
static int iff,ix1,ix2,ix3;
double rand_uniform (int seed)
{
    double rm1,rm2,t;
	int m1,m2,m3,ia1,ia2,ia3,ic1,ic2,ic3,j;
	 
    m1 = 259200, ia1 = 7141, ic1 = 54773, rm1 = 0.0000038580247,
    m2 = 134456, ia2 = 8121, ic2 = 28411, rm2 = 0.0000074373773,
    m3 = 243000, ia3= 4561, ic3 = 51349;
	
    if ((seed < 0) || (iff == 0))
	{
        iff = 1;
        ix1 =  (ic1 - seed) % m1;
        ix1 = (ia1 * ix1 + ic1) % m1;
        ix2 = ix1 % m2;
        ix1 = (ia1 * ix1 + ic1) % m1;
        ix3 = ix1 % m3;
        for (j = 1; j<=97; j++)
		{
            ix1 = (ia1 * ix1 + ic1)% m1;
            ix2 = (ia2 * ix2 + ic2)% m2;
            r[j] = ((double)(ix1) + (double)(ix2) * rm2) * rm1;
        }
    seed = (ia1 * seed + ic2) % m3;
	}
      ix1 = (ia1 * ix1 + ic1) % m1;
      ix2 = (ia2 * ix2 + ic2) % m2;
      ix3 = (ia3 * ix3 + ic3) % m3;
      j = 1 + (int)((97 * ix3) / m3);
      if ((j > 97) || (j < 1)) 
	  {
		  printf("abnormal exit\n");
	      return 1;
	  }
      t = r[j];
      r[j] = ((double)(ix1) + (double)(ix2) * rm2) * rm1;
	  return t;
}


/* SET CURRENT STATE ACCORDING TO SEED. */

void rand_seed( int seed)
{ 
  int j;

  if (!initialized) initialize();

  state->seed = seed;

  state->state48[0] = seed>>16;
  state->state48[1] = seed&0xffff;
  state->state48[2] = rn[0][(seed&0x7fffffff)%Table_size];

  for (j = 0; j<N_tables; j++) 
  { state->ptr[j] = seed%Table_size;
    seed /= Table_size;
  }
}

/* GAUSSIAN GENERATOR.  Done by using the Box-Muller method, but only one
   of the variates is retained (using both would require saving more state).
   See Devroye, p. 235. 

   As written, should never deliver exactly zero, which may sometimes be
   helpful. */

double rand_gaussian (int seed)
//double gasdev(int& idum)
{
	static int iset;
	static double gset;
	double t,v1,v2,fac,r;
    if (iset == 0.0)
	{
		do
		{
			v1 = 2.0 * rand_uniopen (seed) - 1.0;
			v2 = 2.0 * rand_uniopen (seed) - 1.0;
			r = v1 *v1 + v2 *v2;
		}
		while ((r >= 1.0) || (r == 0.0));
		fac = sqrt(-2.0 * log(r) / r);
		gset = v1 * fac;
		t = v2 * fac;
		return t;
		iset = 1;
	}
    else
	{
		t = gset;
		return t;
		iset = 0;
    }
}

/* LOGISTIC GENERATOR.  Just inverts the CDF. */

double rand_logistic (int seed)
{ double u;
  u = rand_uniopen(seed);
  return log(u/(1-u));
}