inline void NState::Evolve(const NState& C0, const ClusterData& cluster, const function2D<double>& exp_)
{
  M.resize(C0.M.size_N(),C0.M.size_Nd());
  istate = C0.istate;
  if (exp_[istate].size()!=M.size_N()) cerr<<"Exponents are not of right size!"<<endl;
  
  // finds biggest entry in the matrix of the state and extract new exponents from the biggest entry
  double max_exp=-100000;
  for (int i=0; i<M.size_N(); i++){
    double maxi = fabs(C0.M[i][0]); // maxi is the biggest entry in the row
    for (int j=1; j<M.size_Nd(); j++)	if (fabs(C0.M(i,j))>maxi) maxi = fabs(C0.M(i,j));
    double expn = log(maxi) + exp_(istate,i);// each row needs to be multiplied with different time evolution
    if (expn>max_exp) max_exp = expn;        // because atomic energies are different within the state
  }
  exponent = max_exp + C0.exponent; // this is the maximal exponent which is assigned to the time evolved state
  for (int i=0; i<M.size_N(); i++){ // The entries in the matrix need to be updated accordingly
    double expo = exp(exp_(istate,i)-max_exp);
    for (int j=0; j<M.size_Nd(); j++) M(i,j) = C0.M(i,j)*expo;
  }
}

inline void NState::apply(const function<function2D<double> >& FM, const function<int>& Fi, const NState& C)
{
  int orig_state = C.istate;
  int new_state = Fi[orig_state];
  istate = new_state;
  exponent = C.exponent;
  //    if (istate!=0) M.SMProduct(FM[orig_state],C.M);
  if (istate!=0) Multiply(M,FM[orig_state],C.M);
}

inline bool NState::empty() const
{
  if (istate==0) return true;
  for (int i=0; i<M.size_N(); i++)
    for (int j=0; j<M.size_Nd(); j++)
      if (fabs(M(i,j))>common::minM) return false;
  return true;
}

inline void NState::Print(ostream& out) const
{
  for (int i=0; i<M.size_N(); i++){
    for (int j=0; j<M.size_Nd(); j++){
      out<<setw(10)<<M(i,j)<<" ";
    }
    out<<endl;
  }
}

inline Number NState::TraceProject(const NState& s){
  if (s.istate!=istate) return 0.0;
  if (M.size_N()!=M.size_Nd()){cerr<<"Non quadratic state in TraceProject!"<<endl;}
  double sum=0;
  for (int i=0; i<M.size_N(); i++) sum += M(i,i);
  return Number(sum,exponent+s.exponent);
}

inline Number NState::Project_to_Pra(const NState& s, function<Number>& Proj){
  if (s.istate!=istate) {
    for (int i=0; i<M.size_N(); i++) Proj[i]=0;
    return 0.0;
  }
  if (M.size_N()!=M.size_Nd()){cerr<<"Non quadratic state in TraceProject!"<<endl;}
  double sum=0;
  for (int i=0; i<M.size_N(); i++){
    sum += M(i,i);
    Proj[i] = Number(M(i,i),exponent+s.exponent);
  }
  return Number(sum,exponent+s.exponent);
}

inline Number NState::ScalarProduct(const NState& s)
{
  if (s.istate!=istate) return 0.0;
  double sum=0;
  for (int i=0; i<M.size_N(); i++)
    for (int j=0; j<M.size_Nd(); j++)
      sum += M(i,j)*s.M(i,j);
  return Number(sum, exponent+s.exponent);
}

inline void NState::SetPraState(int i, const ClusterData& cluster)
{
  istate = cluster.praState(i);
  int msize = cluster.msize(istate);
  M.resize(msize,msize);
  M = 0;
  for (int l=0; l<msize; l++) M(l,l)=1;
  exponent = 0;
}