/**********************************************/
/* This file consists of many miscellaneous functions
/**********************************************/

#include "parameter_sets.h"

///////////////////////////////////////////////////////////////////////////////
// Miscellaneous Functions
///////////////////////////////////////////////////////////////////////////////

//Calculate the z
//int CalcZfactorF(int K, int N) 
//{
//	
//	int mb = 16;
//	int kb = 16;
//	double r = ((double) K)/N;
//	int z;
//	if (r <= 0.5) {
//		// lower rate
//		z = ((N-K)%mb) ? ((N-K)/mb+1) : ((N-K)/mb);
//	} else {
//		// higher rate
//		z = (K%kb) ? (K/kb+1) : (K/kb);
//	}
//	return z;
//}
//
const double const16QAM[2][2][2] = {
// MSB = 0
	1.0/3.16227766016838, 3.0/3.16227766016838,
// MSB = 1
	-1.0/3.16227766016838, -3.0/3.16227766016838,
// LSB = 0
	-1.0/3.16227766016838, 1.0/3.16227766016838,
// LSB = 1
	-3.0/3.16227766016838, 3.0/3.16227766016838
};

const double const64QAM[3][2][4] = {
// MSB = 0
	1.0/6.48074069840786, 3.0/6.48074069840786, 5.0/6.48074069840786, 7.0/6.48074069840786, 
// MSB = 1
	-1.0/6.48074069840786, -3.0/6.48074069840786, -5.0/6.48074069840786, -7.0/6.48074069840786,
// CSB = 0
	-3.0/6.48074069840786, -1.0/6.48074069840786, 1.0/6.48074069840786, 3.0/6.48074069840786,
// CSB = 1
	-7.0/6.48074069840786, -5.0/6.48074069840786, 5.0/6.48074069840786, 7.0/6.48074069840786,
// LSB = 0
	-5.0/6.48074069840786, -3.0/6.48074069840786, 3.0/6.48074069840786, 5.0/6.48074069840786,
// LSB = 1
	-7.0/6.48074069840786, -1.0/6.48074069840786, 1.0/6.48074069840786, 7.0/6.48074069840786
};

//const int PuncPos2[15] = {17, 19, 21, 23, 25, 27, 31, 29, 18, 24, 22, 28, 30, 20, 26};

const double MAX = 10000.0;
double logcothsemi(double x) {
	return -1*log(tanh(x/2));
}

int MaxLogDeMap(int start, int length, double symbol, double * llr) {

	int i, j, k;
	double temp, mindis[2];

	if (length == 2) {
		// Demap the inphase or quadrature part of 16QAM
		for (i=0; i<2; i++) {
			// each bit
			mindis[0] = MAX;
			mindis[1] = MAX;
			for (j=0; j<2; j++) {
				// 0 or 1
				for (k=0; k<2; k++) {
					// find the minimum distance
					temp = pow((symbol-const16QAM[i][j][k]), 2);
					if (temp < mindis[j])
						mindis[j] = temp;
				}
			}
			*(llr+start+i) = mindis[1]-mindis[0];
		}
	} else if (length == 3) {
		// Demap the inphase or quadrature part of 64QAM
		for (i=0; i<3; i++) {
			// each bit
			mindis[0] = MAX;
			mindis[1] = MAX;
			for (j=0; j<2; j++) {
				// 0 or 1
				for (k=0; k<4; k++) {
					// find the minimum distance
					temp = pow((symbol-const64QAM[i][j][k]), 2);
					if (temp < mindis[j])
						mindis[j] = temp;
				}
			}
			*(llr+start+i) = mindis[1]-mindis[0];
		}
	} else {
		// Unsupported modulation types
		return 1;
	}
	return 0;
}

//int ResumeBitsF(struct BasicParaS * ctrl, double * input, double * output)
//{
//	//Declaration
//	int diff_sys = (ctrl->numInBits-ctrl->codeK);
//	int diff_par = ctrl->numChk - ctrl->codeM;
//	int i, j, k, m;
//
//	//Resume the system bits
//	if (ctrl->typeDecode == 0)	//SPA
//	{
//		for (i=0; i<diff_sys; i++)
//		{
//			* (output+i) = 0;	// the punctured systematic bits were all zeros
//		}
//		for (i=0; i<ctrl->codeK; i++)
//		{
//			*(output+i+diff_sys) = *(input+i);
//		}
//	}
//	else		// LSPA and SMA
//	{
//		for (i=0; i<diff_sys; i++)
//		{
//			* (output+i) = 10000.0;	// the punctured systematic bits were all zeros
//		}
//		for (i=0; i<ctrl->codeK; i++)
//		{
//			*(output+i+diff_sys) = *(input+i);
//		}
//	}
//	
//	//Resume the parity bits
//	if (ctrl->codeM>=ctrl->codeK)	// r<=1/2
//	{
//		for (i=0; i<ctrl->codeM; i++)
//		{
//			*(output+ctrl->numInBits+i) = *(input+ctrl->codeK+i);
//		}
//		if (ctrl->typeDecode==0)	//SPA
//		{
//			for (i=0; i<diff_par; i++)
//			{
//				*(output+ctrl->numInBits+ctrl->codeM+i)= 0.5;
//			}
//		}
//		else		// LSPA and SMA
//		{
//			for (i=0; i<diff_par; i++)
//			{
//				*(output+ctrl->numInBits+ctrl->codeM+i) = 0;
//			}
//		}
//	}
//	else		//r>1/2
//	{
//			k = diff_par/ctrl->zfactor;
//			m = diff_par%ctrl->zfactor;
//			// set all the parity as -1
//			for (i=0; i<ctrl->numChk; i++) {
//				*(output+ctrl->numInBits+i) = -1;
//			}
//			// label the bits should be punctured
//			if (ctrl->typeDecode == 0) { // SPA
//				for (i=0; i<k; i++) {
//					for (j=0; j<ctrl->zfactor; j++) {
//						*(output+PuncPos2[i]*ctrl->zfactor+j) = 0.5;
//					}
//				}
//				for (i=0; i<m; i++) {
//					*(output+PuncPos2[k]*ctrl->zfactor+i) = 0.5;
//				}
//			} else { // LSPA and MSA
//				for (i=0; i<k; i++) {
//					for (j=0; j<ctrl->zfactor; j++) {
//						*(output+PuncPos2[i]*ctrl->zfactor+j) = 0;
//					}
//				}
//				for (i=0; i<m; i++) {
//					*(output+PuncPos2[k]*ctrl->zfactor+i) = 0;
//				}
//			}
//			// set the living parity bits
//			j = 0;
//			for (i=0; i<ctrl->numChk; i++) {
//				if ((*(output+ctrl->numInBits+i)) == -1) { // only the unlabeled bits would be covered
//					*(output+ctrl->numInBits+i) = *(input+ctrl->codeK+j);
//					j++;
//				}
//			}
//		}	
//
//	return 0;
//}
//
	///////////////////////////////////////////////////////////////////////////////
	//  random number functions

	// uniform random number generator
	//generate a random number between 0 and 1
	double random()
	{
		long z,k;
		static long s1=12345L;
		static long s2=12345473464L;

		k= s1 / 53668L;
		s1 = 40014L * (s1 -  k*53668L) - k*12211L;
		if (s1<0)
		  s1 = s1 + 2147483563L;
		k = s2 / 52774;
		s2 = 40692L * (s2 - k*52774L) - k*3791L;
		if (s2<0)
	      s2 = s2 + 2147483399L;
		z=s1 - s2;
		if (z<1)
		  z = z + 2147483562L;
		return (double) z / (double) 2147483563.0;
	}

	// generate a number between -1 and 1 
	double generate_random_source()
	{
	
		double sign;		// sign of the number 
		if(random()>=0.5) 
			sign=-1.0;
		else
		  sign=1.0;
	
		return sign*random();
	}

	// binary source generator 
	//generate a sequence of uniformely distributed bits

	short generate_binary_source()
	{
		if(random()<.5) return 0;
		else return 1;
	}


	// Noise generator
	//parameters : sigma2 = variance (global variable)
	//generate a gaussian distributed sequence with zero mean

	double AWGN(double sigma2)
	{
		static int iset=0;
		static double gset;
		double fac,r,v1,v2;

		/* perform Box Muller transformation */
	//	if (iset==0)
		{	
		  do
		  {
		    v1=2.0*random()-1.0;
		    v2=2.0*random()-1.0;
		    r=v1*v1+v2*v2;
		  }
	   	while (r>=1.0 || r == 0.0);
				fac=sqrt(-2.0*sigma2*log(r)/r);
			 gset = v1 * fac;
		 iset = 1;
		return (double)( v2*fac);
		}
	/*
		else
		{
		  iset = 0;
		  return (float)(gset);
		}
	*/
	}