random.c

信道编解码中的turbo编解码程序

/*
*******************************************************************************
* ansi c source code
* file name:
*	Random.c
* abstract:
*	to provide a series of functions, which can offer random numbers, including
*	uniform, exponential and normal distribution.
* reference:
*	
* author:
*	Wang Jiaheng	2004-08-24
* revision history:
*	Wang Jiaheng	2004-08-24	original version
*	Wang Jiaheng	2005-02-03	revised version
*******************************************************************************
*/


/*
*******************************************************************************
*                               include files
*******************************************************************************
*/
#include <stdio.h>
#include <string.h>
#include <math.h>


/*
*******************************************************************************
*			           constants and define declarations
*******************************************************************************
*/
/* for function RunRand1 */
#define IA1		16807
#define IM1		2147483647
#define AM1		(1.0/IM1)
#define IQ1		127773
#define IR1		2836
#define NTAB1	32
#define NDIV1	(1 + (IM1 - 1)/NTAB1)
#define EPS1	(1.2e-7)
#define RNMX1	(1.0 - EPS1)

/* for function RunRand2 */
#define IM2_1	2147483563
#define IM2_2	2147483399
#define AM2		(1.0/IM2_1)
#define IMM2_1	(IM2_1 - 1)
#define IA2_1	40014
#define IA2_2	40692
#define IQ2_1	53668
#define IQ2_2	52774
#define IR2_1	12211
#define IR2_2	3791
#define NTAB2	32
#define NDIV2	(1 + IMM2_1/NTAB2)
#define EPS2	(1.2e-7)
#define RNMX2	(1.0 - EPS2)

/* for function RunRand3 */
#define MBIG	1000000000
#define MSEED	161803398
#define MZ		0
#define FAC		(1.0 / MBIG)

/* for function RunRandBit1 */
#define IB1		1
#define IB2		2
#define IB5		16
#define IB18	131072
#define MASK	(IB1 + IB2 + IB5)


/*
*******************************************************************************
* 		                    global object definition
*******************************************************************************
*/
/* local variables for function RunRand1() */
long	rand1_seed = 1;
long	rand1_iy = 0;
long	rand1_iv[NTAB1];

/* local variables for function RunRand2() */
long	rand2_seed = 1;
long	rand2_idum2 = 123456789;
long	rand2_iy = 0;
long	rand2_iv[NTAB2];

/* local variables for function RunRand3() */
long	rand3_seed = 1;
int		rand3_inext, rand3_inextp;
long	rand3_ma[56];
int		rand3_iff = 0;

/* local variables for function RunRand4() */
long	rand4_seed = 1;
long	rand4_ix, rand4_iy, rand4_iz;
int		rand4_ii = 0;

/* local variables for function RunRandGas() */
int		rand_gas_iset = 0;
double	rand_gas_gset;

/* the shift register for function RunRandBit1() */
long	rand_bit1_reg = 7;


/*
*******************************************************************************
*                            function definition
*******************************************************************************
*/
/*
*******************************************************************************
* description:
*	initialize the local variables for function RunRand1()
* input: 
*	seed: the seed for generating random figures
* output:
*	none
* function reference:
* author:
*	Wang Jiaheng, 2005-02-03	created
*******************************************************************************
*/
void SetRand1(long seed)
{
	rand1_seed = seed;
	rand1_iy = 0;
	memset(rand1_iv, 0, NTAB1 * sizeof(long));
}


/*
*******************************************************************************
* description:
*	to produce uniform iid sequence within [0.0 1.0)
* input: 
*	none
* output:
*	return: a uniformly distributed number within [0.0 1.0)
* function reference:
* author:
*	Wang Jiaheng, 2004-08-24	created
*******************************************************************************
*/
double RunRand1(void)
{
	int		j;
	long	k;
	double	temp;

	if ( (rand1_seed <= 0) || (!rand1_iy) )
	{
		if ( -rand1_seed < 1 ) rand1_seed = 1;
		else rand1_seed = -rand1_seed;

		for (j = NTAB1+7; j >= 0; j--)
		{
			k = rand1_seed / IQ1;
			rand1_seed = IA1 * (rand1_seed - k*IQ1) - IR1*k;

			if ( rand1_seed < 0 ) rand1_seed += IM1;
			if ( j < NTAB1 ) rand1_iv[j] = rand1_seed;
		}

		rand1_iy = rand1_iv[0];
	}

	k = rand1_seed / IQ1;
	rand1_seed = IA1 * (rand1_seed - k*IQ1) - IR1*k;
	if ( rand1_seed < 0 ) rand1_seed += IM1;
	j = rand1_iy / NDIV1;
	rand1_iy = rand1_iv[j];
	rand1_iv[j] = rand1_seed;
	temp = AM1 * rand1_iy;
	if ( temp > RNMX1 ) return RNMX1;
	else return temp;
}


/*
*******************************************************************************
* description:
*	initialize the local variables for function RunRand2()
* input: 
*	seed: the seed for generating random figures
* output:
*	none
* function reference:
* author:
*	Wang Jiaheng, 2005-02-03	created
*******************************************************************************
*/
void SetRand2(long seed)
{
	rand2_seed = seed;
	rand2_idum2 = 123456789;
	rand2_iy = 0;
	memset(rand2_iv, 0, NTAB2 * sizeof(long));
}


/*
*******************************************************************************
* description:
*	to produce uniform iid sequence within [0.0 1.0)
* input: 
*	none
* output:
*	return: a uniformly distributed number within [0.0 1.0)
* function reference:
* author:
*	Wang Jiaheng, 2004-08-24	created
*******************************************************************************
*/
double RunRand2(void)
{
	int		j;
	long	k;
	double	temp;

	if ( rand2_seed <= 0 )
	{
		if ( -rand2_seed < 1 ) rand2_seed = 1;
		else rand2_seed = -rand2_seed;
		rand2_idum2 = rand2_seed;

		for (j = NTAB2 + 7; j >= 0; j--)
		{
			k = rand2_seed / IQ2_1;
			rand2_seed = IA2_1 * (rand2_seed - k*IQ2_1) - k*IR2_1;

			if ( rand2_seed < 0 ) rand2_seed += IM2_1;
			if ( j < NTAB2 ) rand2_iv[j] = rand2_seed;
		}

		rand2_iy = rand2_iv[0];
	}

	k = rand2_seed / IQ2_1;
	rand2_seed = IA2_1 * (rand2_seed - k*IQ2_1) - k*IR2_1;
	if ( rand2_seed < 0 ) rand2_seed += IM2_1;
	k = rand2_idum2 / IQ2_2;
	rand2_idum2 = IA2_2 * (rand2_idum2 - k*IQ2_2) - k*IR2_2;
	if ( rand2_idum2 < 0 ) rand2_idum2 += IM2_2;
	j = rand2_iy / NDIV2;
	rand2_iy = rand2_iv[j] - rand2_idum2;
	rand2_iv[j] = rand2_seed;
	if ( rand2_iy < 1 ) rand2_iy += IMM2_1;
	temp = AM2 * rand2_iy;
	if ( temp > RNMX2 ) return RNMX2;
	else return temp;
}


/*
*******************************************************************************
* description:
*	initialize the local variables for function RunRand3()
* input: 
*	seed: the seed for generating random figures
* output:
*	none
* function reference:
* author:
*	Wang Jiaheng, 2005-02-03	created
*******************************************************************************
*/
void SetRand3(long seed)
{
	rand3_seed = seed;
	rand3_inext = rand3_inextp = 0;
	memset(rand3_ma, 0, 56 * sizeof(long));
	rand3_iff = 0;
}


/*
*******************************************************************************
* description:
*	to produce uniform iid sequence within [0.0 1.0)
* input: 
*	none
* output:
*	return: a uniformly distributed number within [0.0 1.0)
* function reference:
* author:
*	Wang Jiaheng, 2004-08-24	created
*******************************************************************************
*/
double RunRand3(void)
{
	long	mj, mk;
	int		i, ii, k;

	if ( (rand3_seed < 0) || (rand3_iff == 0) )
	{
		rand3_iff = 1;
		mj = labs(MSEED - labs(rand3_seed));
		mj %= MBIG;
		rand3_ma[55] = mj;
		mk = 1;

		for (i = 1; i <= 54; i++)
		{
			ii = (21*i) % 55;
			rand3_ma[ii] = mk;
			mk = mj - mk;
			if ( mk < MZ ) mk += MBIG;
			mj = rand3_ma[ii];
		}

		for (k = 1; k <= 4; k++)
		{
			for (i = 1; i <= 55; i++)
			{
				rand3_ma[i] -= rand3_ma[1 + (i+30) % 55];
				if ( rand3_ma[i] < MZ ) rand3_ma[i] += MBIG;
			}
		}

		rand3_inext = 0;
		rand3_inextp = 31;
		rand3_seed = 1;
	}

	if ( ++rand3_inext == 56 ) rand3_inext = 1;
	if ( ++rand3_inextp == 56 ) rand3_inextp = 1;
	mj = rand3_ma[rand3_inext] - rand3_ma[rand3_inextp];
	if ( mj < MZ ) mj += MBIG;
	rand3_ma[rand3_inext] = mj;
	return mj*FAC;
}


/*
*******************************************************************************
* description:
*	initialize the local variables for function RunRand4()
* input: 
*	seed: the seed for generating random figures
* output:
*	none
* function reference:
* author:
*	Wang Jiaheng, 2005-02-03	created
*******************************************************************************
*/
void SetRand4(long seed)
{
	rand4_seed = seed;
	rand4_ix = rand4_iy = rand4_iz = 0;
	rand4_ii = 0;
}


/*
*******************************************************************************
* description:
*	to produce uniform iid sequence within [0.0 1.0)
* input: 
*	none
* output:
*	return: a uniformly distributed number within [0.0 1.0)
* function reference:
* author:
*	Wang Jiaheng, 2004-08-24	created
*******************************************************************************
*/
double RunRand4(void)
{
     double	temp;

     if ( rand4_ii == 0 )
     {
		rand4_ix = rand4_seed;
		rand4_iy = rand4_seed;
		rand4_iz = rand4_seed;
		rand4_ii++ ;
     }
     rand4_ix = (rand4_ix * 249) % 61967;
     rand4_iy = (rand4_iy * 251) % 63443;
     rand4_iz = (rand4_iz * 252) % 63599;
     temp= (rand4_ix/61967.0 + rand4_iy/63443.0 + rand4_iz/63599.0) - (int)(rand4_ix/61967.0 + rand4_iy/63443.0 + rand4_iz/63599.0);
     return temp;
}


/*
*******************************************************************************
* description:
*	to produce exponential iid sequence with the mean and deviation being 1
* input: 
*	none
* output:
*	return: an exponentially distributed number
* function reference:
* author:
*	Wang Jiaheng, 2004-08-24	created
*******************************************************************************
*/
double RunRandExp(void)
{
	double dum;

	do
	{
		dum = RunRand3();
	} while ( dum == 0.0 );
	return -log(dum);
}


/*
*******************************************************************************
* description:
*	initialize the local variables for function RunRandGas() and RunRandGasMd()
* input: 
*	seed: the seed for generating random figures
* output:
*	none
* function reference:
* author:
*	Wang Jiaheng, 2005-02-03	created
*******************************************************************************
*/
void SetRandGas(long seed)
{
	SetRand3(seed);
	rand_gas_iset = 0;
	rand_gas_gset = 0;
}


/*
*******************************************************************************
* description:
*	to produce normal iid sequence with the mean being 0, and deviation being 1
* input: 
*	none
* output:
*	return: a normally distributed number
* function reference:
* author:
*	Wang Jiaheng, 2004-08-24	created
*******************************************************************************
*/
double RunRandGas(void)
{
	double	fac;
	double	rsq;
	double	v1, v2;

	if ( rand3_seed < 0 ) rand_gas_iset = 0;

	if ( rand_gas_iset == 0 )
	{
		do
		{
			v1 = 2.0 * RunRand3() - 1.0;
			v2 = 2.0 * RunRand3() - 1.0;
			rsq = v1*v1 + v2*v2;
		} while ( rsq >= 1 || rsq == 0.0 );

		fac = sqrt(-2.0 * log(rsq) / rsq);

		rand_gas_gset = v1 * fac;
		rand_gas_iset = 1;
		return (v2 * fac);
	}
	else
	{
		rand_gas_iset = 0;
		return rand_gas_gset;
	}
}


/*
*******************************************************************************
* description:
*	to produce normal iid sequence with the mean being 0, and deviation being 1
* input: 
*	mean and variance
* output:
*	return: a normally distributed number
* function reference:
* author:
*	Wang Jiaheng, 2004-08-24	created
*******************************************************************************
*/
double RunRandGasMd(double mean, double var)
{
	double	fac;
	double	rsq;
	double	v1, v2;

	if ( rand3_seed < 0 ) rand_gas_iset = 0;

	if ( rand_gas_iset == 0 )
	{
		do
		{
			v1 = 2.0 * RunRand3() - 1.0;
			v2 = 2.0 * RunRand3() - 1.0;
			rsq = v1*v1 + v2*v2;
		} while ( rsq >= 1 || rsq == 0.0 );
		fac = sqrt(-2.0 * log(rsq) / rsq);

		rand_gas_iset = 1;
		rand_gas_gset = v1 * fac * sqrt(var) + mean;
		return (v2 * fac * sqrt(var) + mean);
	}
	else
	{
		rand_gas_iset = 0;
		return rand_gas_gset;
	}
}


/*
*******************************************************************************
* description:
*	set the initial register state for m sequence
* input: 
*	reg_state: the initial register state
* output:
*	return: invalid input indication
* function reference:
* author:
*	Wang Jiaheng, 2004-08-26	created
*******************************************************************************
*/
int SetRandBit1(long reg_state)
{
	if ( reg_state <= 0 )
	{
		return 0;
	}

	rand_bit1_reg = reg_state;
	return 1;
}


/*
*******************************************************************************
* description:
*	to produce a pseudo-random sequence with the generating polynomial 
*	x18+x5+x2+x1+1.
* input: 
*	none
* output:
*	return: a bit in the m suquence
* function reference:
* author:
*	Wang Jiaheng, 2004-08-26	created
*******************************************************************************
*/
int RunRandBit1(void)
{
	if ( rand_bit1_reg & IB18 )
	{
		rand_bit1_reg = ((rand_bit1_reg ^ MASK) << 1) | IB1;
		return 1;
	}
	else
	{
		rand_bit1_reg = rand_bit1_reg << 1;
		return 0;
	}
}


/*
*******************************************************************************
* description:
*	set the initial register state for m sequence
* input: 
*	reg_state: the initial register state
* output:
*	p_shift_reg: pointer to the shift register for generating random bits
*	return: invalid input indication
* function reference:
* author:
*	Wang Jiaheng, 2005-02-03	created
*******************************************************************************
*/
int SetRandBit2(unsigned long *p_shift_reg, unsigned long reg_state)
{
	if ( reg_state == 0 )
	{
		return 0;
	}

	*p_shift_reg = reg_state;
	return 1;
}


/*
*******************************************************************************
* description:
*	to produce a pseudo-random sequence with the generating polynomial 
*	x18+x5+x2+x1+1.
* input: 
*	p_shift_reg: pointer to the shift register for generating random bits
* output:
*	return: a bit in the m suquence
* function reference:
* author:
*	Wang Jiaheng, 2005-02-03	created
*******************************************************************************
*/
int RunRandBit2(unsigned long *p_shift_reg)
{
	if ( *p_shift_reg & IB18 )
	{
		*p_shift_reg = ((*p_shift_reg ^ MASK) << 1) | IB1;
		return 1;
	}
	else
	{
		*p_shift_reg = *p_shift_reg << 1;
		return 0;
	}
}