std::complex<double> operator*(const bvec &a, const cvec &b)
{
  it_assert_debug(a.size() == b.size(), "operator*(): sizes does not match");
  std::complex<double> temp = 0;
  for (int i = 0;i < a.size();i++) {temp += (double)a(i) * b(i);}
  return temp;
}

void Hamming_Code::decode(const bvec &coded_bits, bvec &decoded_bits)
{
    int length = coded_bits.length();
    int Itterations = floor_i(static_cast<double>(length) / n);
    ivec Hindexes(n);
    bvec temp(n - k);
    bvec coded(n), syndrome(n - k);
    int isynd, errorpos = 0;
    int i, j;

    decoded_bits.set_size(Itterations*k, false);

    for (i = 0; i < n; i++) {
        for (j = 0; j < n - k; j++)
            temp(j) = H(j, i);
        Hindexes(i) = bin2dec(temp);
    }

    //Decode all codewords
    for (i = 0; i < Itterations; i++) {
        coded = coded_bits.mid(i * n, n);
        syndrome = H * coded;
        isynd = bin2dec(syndrome);
        if (isynd != 0) {
            for (j = 0; j < n; j++)
                if (Hindexes(j) == isynd) {
                    errorpos = j;
                };
            coded(errorpos) += 1;
        }
        decoded_bits.replace_mid(k*i, coded.right(k));
    }
}

bvec encode(Convolutional_Code& nsc, int constraint_length, const bvec& encoder_input, int blockSize, bool verbose)
{
  if (verbose) {cout << "input : " << encoder_input << endl;}
  
  int codedLen = 2 * (blockSize + (constraint_length - 1));
  int nBlocks = encoder_input.length() / blockSize;
  ivec window(blockSize);
  for (int j = 0; j < blockSize; j++) {
    window[j] = j;
  }
  
  bvec nsc_coded_bits(codedLen);
  bvec tr_coded_bits;
  for (int j = 0; j < nBlocks; j++) {
    nsc.encode_tail(encoder_input(window), nsc_coded_bits);
    window = window + blockSize;
    tr_coded_bits = concat(tr_coded_bits, nsc_coded_bits);
  }
  
  // Deal with residual sources if remainder exsists
  if (nBlocks*blockSize != encoder_input.length()) {
    bvec residual_bits = encoder_input.get(nBlocks*blockSize, encoder_input.length()-1);
    nsc.encode_tail(residual_bits, nsc_coded_bits);
    tr_coded_bits = concat(tr_coded_bits, nsc_coded_bits);
  }
  
  if (verbose) {cout << "encoder output: " << tr_coded_bits << endl;}
  return tr_coded_bits;
}

void Rec_Syst_Conv_Code::encode_tail(const bvec &input, bvec &tail, bmat &parity_bits)
{
  int i, j, length = input.size(), target_state;
  parity_bits.set_size(length+m, n-1, false);
  tail.set_size(m, false);
  
  encoder_state = 0;
  for (i=0; i<length; i++) {
    for (j=0; j<(n-1); j++) {
      parity_bits(i,j) = output_parity(encoder_state,2*j+int(input(i)));
    }
    encoder_state = state_trans(encoder_state,int(input(i)));
  }
  
  // add tail of m=K-1 zeros
  for (i=0; i<m; i++) {
    target_state = (encoder_state<<1) & ((1<<m)-1);
    if (state_trans(encoder_state,0)==target_state) { tail(i) = bin(0); } else { tail(i) = bin(1); }
    for (j=0; j<(n-1); j++) {
      parity_bits(length+i,j) = output_parity(encoder_state,2*j+int(tail(i)));
    }
    encoder_state = target_state;
  }
  terminated = true;
}

void BERC::count(const bvec &in1, const bvec &in2)
{
  int countlength = std::min(in1.length(), in2.length()) - std::abs(delay)
                    - ignorefirst - ignorelast;

  if (delay >= 0) {
    for (int i = 0; i < countlength; i++) {
      if (in1(i + ignorefirst) == in2(i + ignorefirst + delay)) {
        corrects++;
      }
      else {
        errors++;
      }
    }
  }
  else {
    for (int i = 0; i < countlength; i++) {
      if (in1(i + ignorefirst - delay) == in2(i + ignorefirst)) {
        corrects++;
      }
      else {
        errors++;
      }
    }
  }
}

bvec Extended_Golay::decode(const bvec &coded_bits)
{
    int no_bits = coded_bits.length();
    int no_blocks = (int)floor((double)no_bits/24);
    bvec output(12*no_blocks);
    int i;
    int j;
    bvec S(12),BS(12),r(12),temp(12),e(24),c(24);
    bmat eyetemp = eye_b(12);

    for (i=0; i<no_blocks; i++) {
	r = coded_bits.mid(i*24,24);
	// Step 1. Compute S=G*r.
	S = G*r;
	// Step 2. w(S)<=3. e=(S,0). Goto 8.
	if( weight(S) <= 3 ) {
	    e = concat(S, zeros_b(12)); goto Step8;
	}
	  
	// Step 3. w(S+Ii)<=2. e=(S+Ii,yi). Goto 8.
	for (j=0; j<12; j++) {
			
	    temp = S + B.get_col(j);
	    if ( weight(temp) <=2 ) {
		e = concat(temp, eyetemp.get_row(j)); goto Step8;
	    }
	}

	// STEP 4. Compute B*S
	BS = B*S;

	// Step 5. w(B*S)<=3. e=(0,BS). Goto8.
	if ( weight(BS) <=3 ) {
	    e = concat(zeros_b(12), BS); goto Step8;
	}

	// Step 6. w(BS+Ri)<=2. e=(xi,BS+Ri). Goto 8.
	for (j=0; j<12; j++) {
	    temp = BS + B.get_row(j);
	    if ( weight(temp) <=2 ) {
		e = concat(eyetemp.get_row(j), temp); goto Step8;
	    }
	}

	// Step 7. Uncorrectable erreor pattern. Choose the first 12 bits.
	e = zeros_b(24); goto Step8;
	  
    Step8: // Step 8. c=r+e. STOP
	c = r + e;
	output.replace_mid(i*12, c.left(12));
    }
  
    return output;
}

void Hamming_Code::encode(const bvec &uncoded_bits, bvec &coded_bits)
{
    int length = uncoded_bits.length();
    int Itterations = floor_i(static_cast<double>(length) / k);
    bmat Gt = G.T();
    int i;

    coded_bits.set_size(Itterations * n, false);
    //Code all codewords
    for (i = 0; i < Itterations; i++)
        coded_bits.replace_mid(n*i, Gt * uncoded_bits.mid(i*k, k));
}

vec mix::process(bvec ce, mat x)
{
 vec y;	
 int N;
 bvec one = ("1");

	#if (DEBUG_LEVEL==3)
	cout << "***** mix::process *****" << endl;	
	cout << "ce=" << ce << endl;
	cout << "x=" << x << endl;
	sleep(1000);
	#endif

	N=ce.length();
	if (x.rows()!=N) { 
		throw sci_exception("mix::process - ce.size <> x.rows()", x.rows() );
	}
	if (x.cols()!=2) { 
		throw sci_exception("mix::process - x=[x1,x1] - x.cols()!=2", x.cols() );
	}

	y.set_length(N);
	for (int i=0; i<N; i++) {
		if ( bool(ce[i])) {
			y0 = G.process(one, to_vec(x(i,0)*x(i,1)))(0);
		}
		y[i]=y0;
	}
	#if (DEBUG_LEVEL==3)
	cout << "y=" << y << endl;
	cout << "+++++ mix::process +++++" << endl;	
	sleep(1000);
	#endif
	return (y);
}

bmat binbuff::generate(bvec ce)
{
 bmat y;
 int K,i;

	#if (DEBUG_LEVEL==3)
	cout << "***** binbuff::generate *****" << endl;	
	cout << "ce=" << ce << endl;
	sleep(SLEEP_TIME_MS);
	#endif

	K=ce.length();
	y.set_size(K,symbol_size);
	for ( i=0; i<K; i++) {
		if ( bool(ce(i))) {
			try {
				get(y0,symbol_size);
				put(y0);
			}
			catch (...) {
				throw sci_exception (" binbuff::generate - error for get()/put()");
			}
		}
		y.set_row(i, y0);
	}
	#if (DEBUG_LEVEL==3)
	cout << "y=" << y << endl;
	cout << "+++++ binbuff::generate +++++" << endl;	
	sleep(SLEEP_TIME_MS);
	#endif
	return (y);
}

cvec fir_x::process(bvec ce, cvec x)
{
 cvec y;	
 int N;

	#if (DEBUG_LEVEL==3)
	cout << "***** fir_x::proc *****" << endl;	
	cout << "ce=" << ce << endl;
	cout << "x=" << x << endl;
	sleep(1000);
	#endif

	N=ce.length();
	if (x.length()!=N) { 
		throw sci_exception("fir_x::process - ce.size <> x.size", x.length() );
	}	
	y.set_size(N);
	for (int i=0; i<N; i++) {
		if (bool(ce[i])) y0=update(x[i]);				
		y[i]=y0;
	}
	

	#if (DEBUG_LEVEL==3)
	cout << "y=" << y << endl;
	cout << "+++++ fir_x::proc +++++" << endl;	
	sleep(1000);
	#endif

	return (y);
}

vec amp::process(bvec ce, vec x)
{
 vec y;	
 int N;

	#if (DEBUG_LEVEL==3)
	cout << "***** amp::process *****" << endl;	
	cout << "ce=" << ce << endl;
	cout << "x=" << x << endl;
	sleep(1000);
	#endif

	N=ce.length();
	if (x.length()!=N) { 
		throw sci_exception("amp::process - ce.size <> x.size()", x.length() );
	}
	y.set_length(N);
	for (int i=0; i<N; i++) {
		if ( bool(ce[i])) {
			y0 = Gain*x[i]+Offset;
		}
		y[i]=y0;
	}
	#if (DEBUG_LEVEL==3)
	cout << "y=" << y << endl;
	cout << "+++++ amp::process +++++" << endl;	
	sleep(1000);
	#endif
	return (y);
}

void binbuff::set_output(bvec yout)
{
	if (yout.size() != symbol_size) { 
		throw sci_exception("binbuff::set_output yout.size <> symbol_size ", symbol_size );
	}
	y0 = yout;
}

bool CRC_Code::decode(const bvec &coded_bits, bvec &out) const
{
  out = coded_bits(0, coded_bits.size() - no_parity - 1);
  if (check_parity(coded_bits)) {
    return true;
  }
  else
    return false;
}

void Reed_Solomon::encode(const bvec &uncoded_bits, bvec &coded_bits)
{
  int i, j, iterations = floor_i(static_cast<double>(uncoded_bits.length())
                                 / (k * m));
  GFX mx(q, k), cx(q, n);
  GFX r(n + 1, n - k);
  GFX uncoded_shifted(n + 1, n);
  GF mpow;
  bvec mbit(k * m), cbit(m);

  coded_bits.set_size(iterations * n * m, false);

  if (systematic)
    for (i = 0; i < n - k; i++)
      uncoded_shifted[i] = GF(n + 1, -1);

  for (i = 0; i < iterations; i++) {
    //Fix the message polynom m(x).
    for (j = 0; j < k; j++) {
      mpow.set(q, uncoded_bits.mid((i * m * k) + (j * m), m));
      mx[j] = mpow;
      if (systematic) {
        cx[j] = mx[j];
        uncoded_shifted[j + n - k] = mx[j];
      }
    }
    //Fix the outputbits cbit.
    if (systematic) {
      r = modgfx(uncoded_shifted, g);
      for (j = k; j < n; j++) {
        cx[j] = r[j - k];
      }
    }
    else {
      cx = g * mx;
    }
    for (j = 0; j < n; j++) {
      cbit = cx[j].get_vectorspace();
      coded_bits.replace_mid((i * n * m) + (j * m), cbit);
    }
  }
}

void BERC::estimate_delay(const bvec &in1, const bvec &in2, int mindelay,
                          int maxdelay)
{
  int num, start1, start2;
  int min_input_length = std::min(in1.length(), in2.length());
  int bestdelay = mindelay;
  double correlation;
  double bestcorr = 0;
  for (int i = mindelay; i < maxdelay; i++) {
    num = min_input_length - std::abs(i) - ignorefirst - ignorelast;
    start1 = (i < 0) ? -i : 0;
    start2 = (i > 0) ?  i : 0;
    correlation = fabs(sum(to_vec(elem_mult(in1.mid(start1, num),
                                            in2.mid(start2, num)))));
    if (correlation > bestcorr) {
      bestdelay = i;
      bestcorr  = correlation;
    }
  }
  delay = bestdelay;
}

void Reed_Solomon::decode(const bvec &coded_bits, bvec &decoded_bits)
{
  bvec   cw_isvalid;
  ivec   erasures(0);
  if (!decode(coded_bits, erasures, decoded_bits, cw_isvalid)) {
    for (int i = 0; i < cw_isvalid.length(); i++) {
      if (!cw_isvalid(i)) {
        decoded_bits.replace_mid(i * k * m, zeros_b(k * m));
      }
    }
  }
}

/* Cumulative error counter */
void BERC::count(const bvec &in1, const bvec &in2)
{
	int countlength = std::min(in1.size(), in2.size()) - std::abs(delay)
					- ignorefirst - ignorelast;
	
	if (delay >= 0) {
		for (int i = 0; i < countlength; i++) {
			if (in1[i + ignorefirst] == in2[i + ignorefirst + delay]) {
				corrects++;
			}else {
				errors++;
			}//end if
		}//end for
	}else {
		for (int i = 0; i < countlength; i++) {
			if (in1[i + ignorefirst - delay] == in2[i + ignorefirst]) {
				corrects++;
			}else {
				errors++;
			}//end if
		}//end for
	}//end if
}

void Rec_Syst_Conv_Code::encode(const bvec &input, bmat &parity_bits)
{
  int i, j, length = input.size();
  parity_bits.set_size(length, n-1, false);
  
  encoder_state = 0;
  for (i=0; i<length; i++) {
    for (j=0; j<(n-1); j++) {
      parity_bits(i,j) = output_parity(encoder_state,2*j+int(input(i)));
    }
    encoder_state = state_trans(encoder_state,int(input(i)));
  }
  terminated = false;
}

// Not optimized for speed
bool CRC_Code::check_parity(const bvec &coded_bits) const
{
  int n = coded_bits.size();
  bvec temp;

  if (reverse_parity) {
    temp = concat(coded_bits.left(n - no_parity), reverse(coded_bits.right(no_parity)));
  }
  else {
    temp = coded_bits;
  }

  for (int i = 0; i < temp.size() - polynomial.size() + 1; i++) {
    if (temp(i) == 1) {
      temp.set_subvector(i, temp(i, i + no_parity) + polynomial);
    }
  }

  if (temp(temp.size() - no_parity, temp.size() - 1) == zeros_b(no_parity))
    return true;
  else
    return false;
}

/* Returns the number of errors between in1 and in2. */
double BERC::count_errors(const bvec &in1, const bvec &in2, int indelay,
                          int inignorefirst, int inignorelast)
{
	int countlength = std::min(in1.size(), in2.size()) - std::abs(indelay)
					- inignorefirst - inignorelast;
	int local_errors = 0;

	if (indelay >= 0) {
		for (int i = 0; i < countlength; i++) {
			if (in1[i + inignorefirst] != in2[i + inignorefirst + indelay]) {
				local_errors++;
			}//end if
		}//end for
	}else {
		for (int i = 0; i < countlength; i++) {
			if (in1[i + inignorefirst - indelay] != in2[i + inignorefirst]) {
				local_errors++;
			}//end if
		}//end for
	}//end if

	return local_errors;
}