#if (!defined(__intutils_h))

#include "capd/capdlib.h"
#include "capd/mpcapdlib.h"
#include "gmp.h"
#include <gmpxx.h>

using namespace std;
using namespace capd;

// for backward compatibility with BIAS/Profil
//inline double Mid(const Interval & x){ return 0.5*(x.rightBound()+x.leftBound());}
inline double Mid(const Interval & x){ return x.mid().leftBound();}
inline double Inf(const Interval & x){ return x.leftBound();}
inline double Sup(const Interval & x){ return x.rightBound();}

inline MpFloat Mid(const MpInterval & x){ return x.mid().leftBound();}
inline MpFloat Inf(const MpInterval & x){ return x.leftBound();}
inline MpFloat Sup(const MpInterval & x){ return x.rightBound();}

DVector Mid(const IVector & x);
DVector Inf(const IVector & x);
DVector Sup(const IVector & x);

MpVector Mid(const MpIVector & x);
MpVector Inf(const MpIVector & x);
MpVector Sup(const MpIVector & x);

DMatrix Mid(const IMatrix & x);
MpMatrix Mid(const MpIMatrix & x);

IMatrix Convert(const DMatrix & x);
MpIMatrix Convert(const MpMatrix & x);

IVector Hull(const IVector & x, const IVector & y);
MpInterval mHull(MpFloat & x);

template <class TVector, class TScalar>
void VectorSetValue(TVector & x, TScalar val){ for (int i=0;i<x.dimension();++i) x[i]=val;}

template <class TMatrix, class TScalar>
void MatrixSetValue(TMatrix & x, TScalar val){
  for (int i=1;i<=x.numberOfRows();++i){
    for (int j=1;j<=x.numberOfColumns();++j) x(i,j)=val;
  }
}

template <class TVector>
typename TVector::ScalarType EuclNorm(const TVector & x){
  typename TVector::ScalarType res=0.0;
  for (int i=0;i<x.dimension();++i) res+=sqr(x[i]);
  return sqrt(res);
}

template <class TVector>
typename TVector::ScalarType EuclDist(const TVector & x, const TVector & y){
  if (x.dimension()!=y.dimension()) return typename TVector::ScalarType(1e100);
  return EuclNorm(x-y);
}

template <class TVector>
typename TVector::ScalarType MaxAbsValue(const TVector & v){
  typename TVector::ScalarType res=0.0;
  for (int d=0;d<v.dimension();++d){
    if (abs(v[d])>res) res=abs(v[d]);
  }
  return res;
}

template <class TVector>
typename TVector::ScalarType VectorAbsMax(const TVector & v){
  typename TVector::ScalarType res=0.0;
  for (int d=0;d<v.dimension();++d){
    if (abs(v[d])>res) res=abs(v[d]);
  }
  return res;
}

template <class TVector>
typename TVector::ScalarType MaxDist(const TVector & x, const TVector & y){
  if (x.dimension()!=y.dimension()) return typename TVector::ScalarType(1e100);
  return MaxAbsValue(x-y);
}

template <class TMatrix>
typename TMatrix::ScalarType MaxAbsValueMatrix(const TMatrix & x){
  typename TMatrix::ScalarType res=0.0;
  for (int i=1;i<=x.numberOfRows();++i){
    for (int j=1;j<=x.numberOfColumns();++j){
      if (abs(x(i,j))>res) res=abs(x(i,j));
    }
  }
  return res;
}

template <class TMatrix>
typename TMatrix::ScalarType Trace(const TMatrix & x){
  typename TMatrix::ScalarType res(0.0);
  for (int i=1;i<=x.numberOfRows();++i) res+=x(i,i);
  return res;
}

template <class TVector>
TVector VectorPart(const TVector & x, int dim, int pos){
  TVector res(dim);
  for (int i=0;i<dim;++i) res[i]=x[pos+i];
  return res;
}

template <class TVector>
TVector VectorRemovePos(const TVector & x, int pos){ // remove selected element from x, returns a vector with the dimension reduced by 1
  int dim=x.dimension();
  TVector res(dim-1,0.0);
  for (int i=0,ip=0;i<dim;++i){
    if (i!=pos) res[ip++]=x[i];
  }
  return res;
}

template <class TVector, class TScalar>
TVector VectorAddPos(const TVector & x, int pos, const TScalar & val){ // adds elemet to x at position pos
  int dim=x.dimension();
  TVector res(dim+1);
  if (pos>=dim){ cerr << "VectorAddPos ERROR: " << endl; return x;}
  for (int i=0, ip=0;i<dim+1;++i){
    if (i==pos) continue;
    res[i]=x[ip];
    ip++;
  }
  res[pos]=val;
  return res;
}

template <class TMatrix>
TMatrix MatrixRemoveRowCol(const TMatrix & M, int r, int c){
  // remove the row r and column c from the matrix M, indexes start with 0
  int rownum=M.numberOfRows();
  int colnum=M.numberOfColumns();
  if (r>=rownum || c>=colnum){ cerr << "MatrixRemoveRowCol ERROR: " << endl; return M;}
  TMatrix Res((rownum-1),(colnum-1),0.0);
  for (int i=0, posi=0;i<rownum;++i){
    if (i==r) continue;
    for (int j=0,posj=0;j<colnum;++j){
      if (j!=c) Res[posi][posj++]=M[i][j];
    }
    posi++;
  }
  return Res;
}

template <class TMatrix>
TMatrix MatrixPart(TMatrix & M, int position){
  int dim=M.numberOfRows();
  if (dim!=M.numberOfColumns()){
    cerr << "MatrixPart: " << "wrong dimensions." << endl;
    return M;
  }
  TMatrix MPart((dim-1),(dim-1),0.0);
  int posi=0;
  for (int i=0;i<dim;++i){
    if (i==position) continue;
    int posj=0;
    for (int j=0;j<dim;++j){
      if (j==position) continue;
      MPart[posi][posj]=M[i][j];
      posj++;
    }
    posi++;
  }
  return MPart;
}

double DiamInterval(Interval & x);
MpFloat DiamInterval(MpInterval & x);
double MaxDiamVector(const IVector & x);
MpFloat MaxDiamVector(const MpIVector & x);
double MaxDiamMatrix(const IMatrix & x);
MpFloat MaxDiamMatrix(const MpIMatrix & x);

#define __intutils_h __intutils_h
#endif