else if ( SNR[w] == 1.0) 
   loop = 1e+1;
   else if ( SNR[w] == 1.5)
   loop = 2e+1;
   else if ( SNR[w] == 2.0) 
   loop = 6e+1;
   else if( SNR[w] == 2.5)
   loop = 2e+2;
   else if ( SNR[w] == 3.0)
   loop = 1e+3;
   else  
   loop = 1e+5;

   while (times < loop){
     Uptr = encoder.get_U();
     for ( i = 0; i < k; i++)
       for ( j = 0; j < k; j++)
	 U[i][j] = Uptr[58*i + j];
     
     Y_H_Vptr = encoder.get_Y_H_V();
     for ( i = 0; i < n; i++)
       for ( j = 0; j < n; j++)
	 Y_H_V[i][j] = Y_H_Vptr[65*i + j];
     
     iter = 0;
     while ( iter < iteration)
       {
	 
	 for (i = 0; i < n; i++)
	   for (j = 0; j < n; j++)
	     Lu[i][j] = Le_V[i][j];
	   
 // ##### Below is the computation of Forwardmatric and Backwardmatric for Horizontal ##### 

     node[0].setForwardmatric(1);

     for ( d = 0; d < n; d++)
       {
	 for (i = 0; i < n; i++){
	   for (j = 0; j < n; j++){
	     Lxy[i][j] = Lc*Y_H_V[i][j] + Lu[i][j];
	   }
	 }
	 for (i = 1; i < n + 1 ; i++){  //start to calaulate alpha, i is each time point
	   for ( j = 1; j <= ( int)pow(2, statecomplexfile[i]); j++)
	     {
	       counterbranch = startbranchindex[i-1];
	       nodenum = startnodeindex[i] - 1 + j;
	       temp1 = 0.0; temp2 = 0.0; t = 0.0; counter = 0; 
	       while ( counterbranch < startbranchindex[i])
		 {
		   if (branch[counterbranch].getDestinationnode() == nodenum)
		     {
		       counter++; 
		       if (counter == 1){ 
			 temp1 = log(node[branch[counterbranch].getOriginalnode()].getForwardmatric()) + 0.5*branch[counterbranch].getBit()*Lxy[d][i-1];
			 temp2 = exp(temp1);}
		       else if (counter == 2)
			 { 
			   temp1 = log(node[branch[counterbranch].getOriginalnode()].getForwardmatric()) + 0.5*branch[counterbranch].getBit()*Lxy[d][i-1];
			   temp2 = temp2 + exp(log(node[branch[counterbranch].getOriginalnode()].getForwardmatric()) + 0.5*branch[counterbranch].getBit()*Lxy[d][i-1]);
			 }
		       counterbranch++;
		     }
		   else 
		     counterbranch++;
		 }      
	       node[nodenum].setForwardmatric(temp2);
	     }
	 }
	 
	 //+++++++++++++++++++ Below is the case for the last node: time(i) = n ++++++++++++++++++++++++
	
	 temp1 = log(node[branch[counterbranch].getOriginalnode()].getForwardmatric()) + 0.5*branch[counterbranch].getBit()*Lxy[d][n-1];
	 temp2 = log(node[branch[counterbranch + 1].getOriginalnode()].getForwardmatric()) + 0.5*branch[counterbranch + 1].getBit()*Lxy[d][n-1];
	 temp2 = exp(temp1) + exp(temp2);
	 node[nodenum].setForwardmatric(temp2);
	
//========================== start the computation of Beta ===========================
 
	  node [nodenumber].setBackwardmatric(1);// nodenumber is the number of the last node
	 temp1 = 0.0; temp2 = 0.0; t = 0.0; counter = 0;
	 i = n-1;     //start to calaulate beta for backward, i is each time point
	 for ( j = 0; j < ( int)pow(2, statecomplexfile[i]); j++)
	   {
	     counterbranch = startbranchindex[i];// there are 2 branches going to the last node
	     nodenum = startnodeindex[i] + j;
	     temp1= log(node[branch[counterbranch+j].getDestinationnode()].getBackwardmatric()) + 0.5* branch[counterbranch+j].getBit()*Lxy[d][n-1];
	     temp2 = exp(temp1);
	     node[nodenum].setBackwardmatric(temp2);
	   }
	 for (i = n - 2; i > -1; i--)
	   {
	     for (j = 0; j < (int)pow(2, statecomplexfile[i]); j++)
	       {
		 counterbranch = startbranchindex[i];
		 
		 nodenum = startnodeindex[i] + j;
		 counter = 0;
		 while ( counterbranch < startbranchindex[i+1])
		   {
		     if (branch[counterbranch].getOriginalnode() == nodenum)
		       {
			 counter++;
			 if (counter == 1){
			   temp1 = log(node[branch[counterbranch].getDestinationnode()].getBackwardmatric()) + 0.5*branch[counterbranch].getBit()*Lxy[d][i];
			   temp2 = exp(temp1);
			 }
			 else if (counter == 2)
			   {
			     temp1 = log(node[branch[counterbranch].getDestinationnode()].getBackwardmatric()) + 0.5*branch[counterbranch].getBit()*Lxy[d][i];
			     temp2 = temp2 + exp(log(node[branch[counterbranch].getDestinationnode()].getBackwardmatric()) + 0.5*branch[counterbranch].getBit()*Lxy[d][i]);
			   }
			 counterbranch++;
		       }
		     else
		       counterbranch++;
		   }
		 node[nodenum].setBackwardmatric(temp2);
	       }
	   }
	
	 //++++++ Below is the computation for the Extrinsic value with forward tracing value and backward tracing value ++++++
	 
	 double Le_H_positive[65][65], Le_H_negative[65][65];
	 for (i =0; i < n; i++)
	   for ( j = 0; j < n; j++)
	     Le_H_positive[i][j] = Le_H_negative[i][j] = 0;
	 for (i =0; i < n; i++){
	   if (statecomplexfile[i+1] >= statecomplexfile[i]){
	     counterbranch = startbranchindex[i];
	     while ( counterbranch < startbranchindex[i+1])
	       {
		 if (branch[counterbranch].getBit() == 1){
		   Le_H_positive[d][i] =Le_H_positive[d][i] + ( node[branch[counterbranch].getOriginalnode()].getForwardmatric())* (node[branch[counterbranch].getDestinationnode()].getBackwardmatric());
		   counterbranch++;}
		 else {
		   Le_H_negative[d][i] =Le_H_negative[d][i] + ( node[branch[counterbranch].getOriginalnode()].getForwardmatric())* (node[branch[counterbranch].getDestinationnode()].getBackwardmatric());
		   counterbranch++;}
	       }
	     Le_H[d][i] = log(Le_H_positive[d][i]/Le_H_negative[d][i]);
	   }
	 }
       }
    
     for (i = 0; i < n; i++) //assign the Horizontal softvalues to Lu
       for (j = 0; j < n; j++)
	 Lu[i][j] = Le_H[i][j];
     
 //##### Below is the computation of Forwardmatric and Backwardmatric for Vertical######
     
    for (i = 0; i < n; i++){
       for (j = 0; j < n; j++){
	 Lxy[j][i] = Lc*Y_H_V[i][j] + Lu[i][j];
       }
    }
    for (i = 0; i < n; i++)
      for (d = 0; d < n; d++){
	for (i = 1; i < n + 1 ; i++){  //start to calaulate alpha, i is each time point
	  for ( j = 1; j <= ( int)pow(2, statecomplexfile[i]); j++)
	    {
	     counterbranch = startbranchindex[i-1];
	     nodenum = startnodeindex[i] - 1 + j;
	     temp1 = 0.0; temp2 = 0.0; t = 0.0; counter = 0; 
	     while ( counterbranch < startbranchindex[i])
	       {
		 if (branch[counterbranch].getDestinationnode() == nodenum)
		   {
		     counter++; 
		     if (counter == 1){ 
		       temp1 = log(node[branch[counterbranch].getOriginalnode()].getForwardmatric()) + 0.5*branch[counterbranch].getBit()*Lxy[d][i-1];
		       temp2 = exp(temp1);}
		     else if (counter == 2)
		       { 
			 temp1 = log(node[branch[counterbranch].getOriginalnode()].getForwardmatric()) + 0.5*branch[counterbranch].getBit()*Lxy[d][i-1];
			 temp2 = temp2 + exp(log(node[branch[counterbranch].getOriginalnode()].getForwardmatric()) + 0.5*branch[counterbranch].getBit()*Lxy[d][i-1]);
		       }
		     counterbranch++;
		   }
		 else 
		   counterbranch++;
	       }      
	     node[nodenum].setForwardmatric(temp2);
	   }
	}
	//============ below is the case for the last node: time(i) = n ===============
	
       temp1 = log(node[branch[counterbranch].getOriginalnode()].getForwardmatric()) + 0.5*branch[counterbranch].getBit()*Lxy[d][n-1];
       temp2 = log(node[branch[counterbranch + 1].getOriginalnode()].getForwardmatric()) + 0.5*branch[counterbranch + 1].getBit()*Lxy[d][n-1];
       temp2 = exp(temp1) + exp(temp2);
       node[nodenum].setForwardmatric(temp2);
      
       // ========================== start the computation of Beta ====================	 
       
       node [nodenumber].setBackwardmatric(1);
       temp1 = 0.0; temp2 = 0.0; t = 0.0; counter = 0;
       i = n-1;  //start to calaulate beta for backward, i is each time point
       for ( j = 0; j < ( int)pow(2, statecomplexfile[i]); j++)
	 {
	   counterbranch = startbranchindex[i];// there are 2 branches going to the last node
	   nodenum = startnodeindex[i] + j;
	   temp1= log(node[branch[counterbranch+j].getDestinationnode()].getBackwardmatric()) + 0.5* branch[counterbranch+j].getBit()*Lxy[d][n-1];
	   temp2 = exp(temp1);
	       node[nodenum].setBackwardmatric(temp2);
	 }
       for (i = n - 2; i > -1; i--)
	 {
	   for (j = 0; j < (int)pow(2, statecomplexfile[i]); j++)
	     {
	       counterbranch = startbranchindex[i];
	       nodenum = startnodeindex[i] + j;
	       counter = 0;
	       while ( counterbranch < startbranchindex[i+1])
		 {
		   if (branch[counterbranch].getOriginalnode() == nodenum)
		     {
		       counter++; 
		       if (counter == 1){
			 temp1 = log(node[branch[counterbranch].getDestinationnode()].getBackwardmatric()) + 0.5*branch[counterbranch].getBit()*Lxy[d][i];
			 temp2 = exp(temp1);
		       }
		       else if (counter == 2)
			 {
			   temp1 = log(node[branch[counterbranch].getDestinationnode()].getBackwardmatric()) + 0.5*branch[counterbranch].getBit()*Lxy[d][i];
			   temp2 = temp2 + exp(log(node[branch[counterbranch].getDestinationnode()].getBackwardmatric()) + 0.5*branch[counterbranch].getBit()*Lxy[d][i]);
			 }
		       counterbranch++;
		     }
		   else
		     counterbranch++;
		 }
	       node[nodenum].setBackwardmatric(temp2);
	      
	     }
	 }
   
//++++++ Below is the computation for the Extrinsic value with forward tracing value and backward tracing value +++++++
      
       double Le_V_positive[65][65], Le_V_negative[65][65]; 
       for (i =0; i < n; i++)
	 for ( j = 0; j < n; j++)
	   Le_V_positive[i][j] = Le_V_negative[i][j] = 0;

       for (i =0; i < n; i++){
	 if (statecomplexfile[i+1] >= statecomplexfile[i]){
	   counterbranch = startbranchindex[i];
	   while ( counterbranch < startbranchindex[i+1])
	     {
	       if (branch[counterbranch].getBit() == 1){
		 Le_V_positive[d][i] =Le_V_positive[d][i] + ( node[branch[counterbranch].getOriginalnode()].getForwardmatric())* (node[branch[counterbranch].getDestinationnode()].getBackwardmatric());
		 counterbranch++;}
	       else {
		 Le_V_negative[d][i] =Le_V_negative[d][i] + ( node[branch[counterbranch].getOriginalnode()].getForwardmatric())* (node[branch[counterbranch].getDestinationnode()].getBackwardmatric());
		 counterbranch++;}
	     }
	   Le_V[i][d] = log(Le_V_positive[d][i]/Le_V_negative[d][i]);
	 }
       }
     }

       iter++;}

     for (i = 0; i < n; i++)
       for (j = 0; j < n; j++)
	 LL_U[i][j] = Le_H[i][j] + Le_V[i][j] + Lc*Y_H_V[i][j];
	 
//++++++++++++++++++++ Below is the make a decision for the decoder +++++++++++++++++++++++
     
     int counter_i = 0, counter_j = 0;
     for ( i = 0; i < n; i++){ 
       if  ( statecomplexfile[i+1] >= statecomplexfile[i])
	 {   
	   counter_j = 0;
	   for (j = 0; j < n; j++)
	     {
	       
	       if (statecomplexfile[j+1] >= statecomplexfile[j])
		 {
		   if ( LL_U[i][j] > 0)
		     U_estimate[i - counter_i][j - counter_j] = 1;
		   else
		     U_estimate[i - counter_i][j - counter_j] = 0;
		 }
	       else 
		 counter_j++;
	     }
	 }
       else
	 counter_i++;
     }
     
 // +++++++++++++++++++++++++++  computation for the error bit rate  +++++++++++++++++++++++++++++++
     
     for ( i = 0; i < k; i++){
       for (j = 0; j < k; j++)
	 if ( U[i][j] != U_estimate[i][j])
	   error++;
     }
    
     times++;
   }
      
  
   
 Pe = (double)error/(k*k*times);
 tout << "Number of Iterations = " << iter <<endl;
 tout << "Number of Bits in Error = " << error << endl;
 tout << "Number of Transmitted Bits = " << times*k*k << endl;
 tout << "Bit Error Rate (BER) = " << Pe << endl << endl; 
 }


 tout.close();
 return 0;
}

double findstate(int *index,int length, int Gx[][65],int RM_N, int RM_K)
{
  int i;
  int* Temp = 0;
  Temp = new int[length];
  double result = 0.0;
  if(index[0] == -2)
    {
      for(i=0; i < length; i++)
	Temp[i] = 0;
    }
  else if (index[0] >= 0)
  {
    for(i=0; i < length; i++) 
      Temp[i] = Gx[index[0]][i];
  }
  
  for(i=0; i < RM_K; i++)
    {
      if(index[i+1] == -1)
	break;
      else if(index[i+1] >= 0)
	{
	  newAddXYtoX(Temp,Gx[index[i+1]],length,RM_N);
	}
    }
  
  for(i = length - 1; i >= 0; i--)
      {
	if (Temp[i] == 1){
	  result += pow(2, length - 1 - i);
	}
      }
  delete[] Temp;
  return result;
}

void newAddXYtoX(int *X,int *Y,int leng, int RM_N)
{
  int i; 
  
  int temp[65];
  for(i=0; i < 65; i++)
    temp[i] = 0;
  for(i=0; i < leng; i++){
    temp[i] = X[i];
  }
  for(i = RM_N-leng; i < RM_N; i++)
    {
      X[i-(RM_N-leng)] = temp[i-(RM_N-leng)]^Y[i];
    }
}

void fillid(int *in,int *out,int s,int a, int RM_N)
{
  int i;
  int num = 0;
  for(i=0; i < RM_N; i++)
    {
      if(in[i] == -1)
	break;
      else
	{
	  out[num] = in[i];
	  num++;
	}
    }
  if(a == 1)
    out[num] = s;
}

int findbit(int originaltime,double original,int *index, int Gx[][65], int RM_N, int RM_K)
{
  int i; 
  int arraytemp[65];
  double temp;
  int Temp[65];
  
  temp = original;
  for(i=RM_N-originaltime-1; i >= 0; i--)
    {
      arraytemp[i] = (int)fmod(temp, 2);
      temp = temp / 2;
    }
  if(index[1] == -1)
    {
      return arraytemp[0] * 2 - 1;
    }
  else
    {
      if(index[0] == -2)
	{
	  for(i=0; i < RM_N-originaltime+1; i++)
	    Temp[i] = 0;
	}
      else if(index[0] >= 0)
	{
	  for(i=0; i < RM_N-originaltime+1; i++)
	    Temp[i] = Gx[index[0]][i];
	}
      
      for(i=0; i < RM_K; i++)
	{
	  if(index[i+1] == -1 )
	    break;
	  else if(index[i+1] >= 0)
	    {
	      newAddXYtoX(Temp,Gx[index[i+1]],RM_N-originaltime,RM_N);
	    }
	}
      return Temp[0]*2-1;
    } 
}