//**************************************************************************************//
// This program simulates the performance of the shortened RM(32,26)^2-turbo code case A
// The Trellis-based MAP algorithm is used. This code is an Equal Error Protetion (EEP) code
//**************************************************************************************//

#include <iostream.h>
#include <fstream.h>
#include <string.h>
#include <stdlib.h>
#include <time.h>
#include <alloc.h>
#include <math.h>
#include <iomanip.h>

#define RM_N 32
#define RM_K 26
#define ORDER 3
#define dimfill1 8
#define dimfill2 26
#define Maxiter 5
#define startEBN0 1
#define SNRincr 0.5
#define Numrun 20

int Gx[27][33]; 
int nzero;
void get_data(void);
void add_noise(double sig);
double gauss_rv(double mean, double sigma);
int print_e_r(void);
double prioriprob[33][33];
double Receive[33][33]; 
int DIN[27][33];
double Soft[33][33];
double FinalSoft[33][33];
int numcode;

// RM generator Matrix part
void generateG(int order);
void multiply(int *out, int *v1,int *v2);
double bico(int n, int k);
int factorial(int n);


// for MSGM part
void Find_L_R(int *L, int *R, int *Given_row,int n);
void AddXYtoX(int *X, int *Y);
char Check(int *L, int *R);
void MSGM(int Gen[][33]);
void RMencoder(double *Out, int *In);
void Echelon(int Gen[][33]);
int Left[27];
int Right[27];
int Parity[27];
long total_errors[7];
long total_errorsarea1[7];
long total_errorsarea2[7];
double BERarea1[6],BERarea2[6],BERarea3[6];
double BER[6];
ofstream fout;
int iter;
double total_number_of_data_bits;
double Output_EBN0_dB;
double CodeRate,CodeRate_dB;

class Node {
 public:
  Node();
  int getTime(void);
  int getStage(void);
  void setStage(int s);
  void setTime(int t);
  void setForwardmatric(double fm);
  void setBackwardmatric(double bm); 
  double getForwardmatric(void);
  double getBackwardmatric(void);
 private: 
  int time;
  int stage;
  double forwardmatric;
  double backwardmatric;
};

Node::Node() {time=stage=0; forwardmatric = 0.0; backwardmatric = 0.0;}

int Node::getTime(void) {return time;}
int Node::getStage(void) {return stage;}
double Node::getForwardmatric(void) {return forwardmatric;}
double Node::getBackwardmatric(void) {return backwardmatric;}
void Node::setStage(int s) {stage = s;}
void Node::setTime(int t) {time = t;}
void Node::setForwardmatric(double fm) {forwardmatric = fm;}
void Node::setBackwardmatric(double bm) {backwardmatric = bm;} 

class Branch {
 public:
  Branch();
  int getOriginaltime(void);
  int getOriginal(void);
  int getDestination(void);
  int getBit(void);
  int getOriginalnode(void);
  int getDestinationnode(void);
  void setOriginaltime(int ot);
  void setOriginal(int o);
  void setDestination(int d);
  void setBit(int b);
  void setOriginalnode(int on);
  void setDestinationnode(int dn);
 private:
  int originaltime;
  int original; 
  int originalnode;
  int destinationnode;
  int destination;
  int bit;
};

Branch::Branch() {originaltime=original=destination=bit=0;}

int Branch::getOriginaltime(void) {return originaltime;}
int Branch::getOriginal(void) {return original;}
int Branch::getDestination(void) {return destination;}
int Branch::getBit(void) {return bit;}
int Branch::getOriginalnode(void) {return originalnode;}
int Branch::getDestinationnode(void) {return destinationnode;}
void Branch::setOriginaltime(int ot) {originaltime = ot;}
void Branch::setOriginal(int o) {original = o;}
void Branch::setDestination(int d) {destination = d;}
void Branch::setBit(int b) {bit = b;}
void Branch::setOriginalnode(int on) {originalnode = on;}
void Branch::setDestinationnode(int dn) {destinationnode = dn;}

void Trellis(void);
int findbit(int originaltime,int original,int *index);
void newAddXYtoX(int *X,int *Y,int leng);
int findstate(int *index,int length);
void fillid(int *in,int *out,int s,int a);
void Softout(double *Rx,double Lc, double *Le_in, double *soft_out,int tag); 
int findparallel(int originaltime);
double callogone(double in1,double in2);
double callogzero(double in1,double in2);
double algebra(double in1,double in2);
double Lc;
double LeIn[33];
double Le[33];
double infosoft[27];
int iteration;
int w[33][27];
Node node[1000];
Branch branch[2000];
int statecomplex[33];
int lastbranch, lastnode;
int testindex[27];
int oneindex[27]; 
int numoneindex = 0;
int numindex = 0;
int bin[27];
int startnodeindex[33];
int startbranchindex[33];

int main(void)
{
  short int FieldPoly[9];
  time_t now;
  int i,j,k,e;
  double temp;
  int shift;
  int iEBN0;
  long MaxSample,Sample;
  double EbN0,EbN0_dB,sigma;

  fout.open("BTC_ATM.dat");
  srand((unsigned)time(&now));
  cout << setprecision(3);

  //***********************************************************//
  //                   RM(32,26) Generator
  //***********************************************************//

  generateG(ORDER);

  //*********************************************************//
  // Transform the Generator matrix to be in reduced echelon form  
  //*********************************************************//
  MSGM(Gx);

 //***********************************************************//
  //     Minimal Trellis Construction using Massey Algorithm
  //***********************************************************//

  Trellis();

  fout << setprecision(3);
  nzero = (int) (dimfill1*dimfill2)/RM_K;

  for (iEBN0=0;iEBN0<4;iEBN0++)
    {
      Output_EBN0_dB = startEBN0 + iEBN0*SNRincr;
      CodeRate= (float) (RM_K*RM_K-dimfill1*dimfill2) / (float) (RM_N*RM_K+RM_K*(RM_N-RM_K)-dimfill1*RM_N);
      Lc = 4*CodeRate*pow(10.0,Output_EBN0_dB/10.0);
      CodeRate_dB = 10*log(CodeRate)/log(10.0);
      EbN0_dB = Output_EBN0_dB + CodeRate_dB;
      cout  << "EBNO db = " << Output_EBN0_dB << endl;
      fout  << "EBNO db = " << Output_EBN0_dB << endl;
      EbN0 = exp(log(10.0)*(EbN0_dB/10.0));
      sigma = sqrt(1.0/(2.0*EbN0));
      
      cout << "CodeRate = " << CodeRate << endl;
      fout << "CodeRate = " << CodeRate << endl;
      MaxSample = (1+iEBN0)*Numrun;
      
      for(i=0; i < Maxiter; i++)
	{
	  total_errorsarea1[i] = 0;
	  total_errorsarea2[i] = 0;
	}
	

      double TempSoft1[RM_K][RM_N],TempSoft2[RM_K][RM_N];
      double TempReceive[RM_K][RM_N];
      double Softhor[RM_K][RM_K],Softver[RM_K][RM_K];
      
      for(Sample = 0; Sample < MaxSample; Sample++)
	{
	  // Initialization

	  for(i=0; i < RM_K; i++)
	    for(j=0; j < RM_K; j++)
	      {
		prioriprob[i][j] = 0;
		Softhor[i][j] = 0;
		Softver[i][j] = 0;
	      }
	  

	  //***********************************************************//
	  //            Random Data and Channel Encoder 
	  //***********************************************************//
	 
	  get_data(); 

	  //***********************************************************//
	  //                  AWGN Channel
	  //***********************************************************//
	 
	  add_noise(sigma);
	  
	  //***********************************************************//
	  //Set shortened bits to be -2, they are known at the receiver
	  //***********************************************************//
	  int num = 0;
	  for(i=RM_K-1; i >= 0; i--)
	    for(j=0; j < RM_K; j++)
	      if(num < dimfill1*dimfill2)
		{
		  Receive[i][Left[j]] = -2;
		  num++;
		}

	  for(i=RM_K-1; i >= RM_K-nzero; i--)
	    for(j=0; j < RM_N; j++)
	      Receive[i][j] = -2;

	  
	  for(iteration = 0; iteration < Maxiter; iteration++)
	    {
	      
	      for(i=0; i < RM_K; i++)
		Softout(Receive[i],Lc,prioriprob[i],Soft[i],0);

	      for(i=0; i < RM_K; i++)
		for(j=0; j < RM_K; j++)
		  TempSoft1[i][j] = Softhor[j][i] = Soft[j][i];


	      for(i=0; i < RM_K; i++)
		{
		  for(j=0; j < RM_K; j++)
		    TempReceive[i][Left[j]] = Receive[j][Left[i]];
		  
		  for(j=0; j < RM_N-RM_K; j++)
		    TempReceive[i][Parity[j]] = Receive[RM_K+j][Left[i]];
		}

	      for(i=0; i < RM_K; i++)
		Softout(TempReceive[i],Lc,TempSoft1[i],TempSoft2[i],1);
	   

	      for(i=0; i < RM_K; i++)
		for(j=0; j < RM_K; j++)
		  prioriprob[i][j] = Softver[i][j] = TempSoft2[j][i];
	 

	      for(i=0; i < RM_K; i++)
		for(j=0; j < RM_K; j++)
		  FinalSoft[i][j] = Lc*Receive[i][Left[j]] + Softhor[i][j] + Softver[i][j];

	      e = print_e_r();
	    
	    }// end iter
	
	}// end Sample

      total_number_of_data_bits = (1+iEBN0)*Numrun*(RM_K*RM_K-dimfill1*dimfill2);
      cout << "total_bit = " << total_number_of_data_bits << endl;

    
      for (i=0; i < Maxiter; i++)
	{
	  cout << "total_errorsarea1[ " << i << "] = " << total_errorsarea1[i] << endl;
	  cout << "total_errorsarea2[ " << i << "] = " << total_errorsarea2[i] << endl << endl;
	  
	  
	  fout << "total_errorsarea1[ " << i << "] = " << total_errorsarea1[i] << endl;
	  fout << "total_errorsarea2[ " << i << "] = " << total_errorsarea2[i] << endl << endl;
	  
	}
      
      for (i=0;i<=Maxiter;i++) BERarea1[i]=(float)total_errorsarea1[i]/(float)((1+iEBN0)*Numrun*(RM_K/2)*(RM_K-nzero));
      for (i=0;i<=Maxiter;i++) BERarea2[i]=(float)total_errorsarea2[i]/(float)((1+iEBN0)*Numrun*(RM_K*RM_K-dimfill1*dimfill2-(RM_K/2)*(RM_K-nzero)));
      
      for (i=0;i<=Maxiter;i++) BER[i]=(float)(total_errorsarea1[i]+total_errorsarea2[i])/(float)total_number_of_data_bits;
      
      
      for (i=0; i < Maxiter; i++) cout << "BERarea1[" << i << "] = " << BERarea1[i] << "  ";
      cout << endl;
      for (i=0; i < Maxiter; i++) cout << "BERarea2[" << i << "] = " << BERarea2[i] << "  ";
      cout << endl;  
      
      for (i=0; i < Maxiter; i++) cout << "BER[" << i << "] = " << BER[i] << " ";
      cout << endl;
      
   
      for (i=0; i < Maxiter; i++) fout << "BERarea1[" << i << "] = " << BERarea1[i] << " ";
      fout << endl;
      for (i=0; i < Maxiter; i++) fout << "BERarea2[" << i << "] = " << BERarea2[i] << " ";  
      fout << endl;
      
      for (i=0; i < Maxiter; i++) fout << "BER[" << i << "] = " << BER[i] << " ";
      fout << endl;
    
      cout << endl;
      
    }// end iEBN0
  
  fout.close();
}

void get_data(void)
{
  int i,j;
  int temp;
  double u1;
  int iparity = 0;
  int templeft[RM_N];

  for(i=0; i < RM_N; i++)
    templeft[i] = 0;

  for(i=0; i < RM_K; i++)
    templeft[Left[i]] = 1;
  
  for(i=0; i < RM_N; i++)
    {
      if(templeft[i] == 0)
	{
	  Parity[iparity] = i;
	  iparity++;
	} 
    }

  // Random data
  for(i=0; i < RM_K ; i++)
    {
      for(j=0; j < RM_K; j++)
	{
	  temp = rand();
	  u1 = (float) temp;
	  u1 /= RAND_MAX;
	  
	  if(u1 >= 0.5)
	    DIN[i][j] = 1;
	  else
	    DIN[i][j] = 0;
	}
    }

  // Shortening Pattern A
  int num = 0;
  for(i=RM_K-1; i >= 0; i--)
    for(j=0; j < RM_K; j++)
      if(num < dimfill1*dimfill2)
	{
	  DIN[i][j] = 0;
	  num++;
	}

  // Channel Encoder
  for(i=0; i < RM_K; i++)
    RMencoder(Receive[i],DIN[i]);
  
  int Itemp[RM_K][RM_K];
  double Ctemp[RM_K][RM_N];
  
  for(i=0; i < RM_K; i++)
    for(j=0; j < RM_K; j++)
      Itemp[i][j] = DIN[j][i];

  for(i=0; i < RM_K; i++)
    RMencoder(Ctemp[i],Itemp[i]);
  
  for(i=0; i < RM_K; i++)
    for(j=0; j < RM_N-RM_K; j++)
      Receive[RM_K+j][Left[i]] = Ctemp[i][Parity[j]];

}

void add_noise(double sig)
{
  int i, j;
  
  for(i=0; i < RM_K; i++)
    for(j=0; j < RM_N; j++)
      Receive[i][j] = Receive[i][j] + gauss_rv(0.0,sig);
  
  for(i=RM_K; i < RM_N; i++)
    for(j=0; j < RM_K; j++)
      Receive[i][Left[j]] = Receive[i][Left[j]] + gauss_rv(0.0,sig);
  
}

double gauss_rv(double mean, double sigma)
{
  int i1,i2;
  double temp,u1,u2,g1;
  
  i1 = rand();
  u1 = (float)i1;