class common{
public:
  static int N_flavors;
  static int N_ifl;
  static double beta;
  static double U;
  static double mu;
  static int max_size;
  static double minM; // If matrix is singular, it destroys sampling. Very small determinant should never be accepted.
  static double minD; // Smallest determinant should not be accepted
  static int tsample;
  static int warmup;
  static int CleanUpdate;
  static double PChangeOrder;
  static int Ncout;
  static int Naver;
  static int my_rank;
  static int mpi_size;
  static double TwoKinks;
  static int GlobalFlip;
  static double treshold;
  static int SampleGtau;
  static int PreciseP;
  static double minDeltat;
  static bool SampleSusc;
  static void SetParameters(int rank, int size, double mu_, double U_, double beta_, int max_size_, int N_flavors_,
			    int N_ifl_, int tsample_,
			    int warmup_, int CleanUpdate_, double minM_, double minD_, /*double neglectG_, */double PChangeOrder_,
			    int Ncout_, int Naver_, double TwoKinks_, int GlobalFlip_, double treshold_, int SampleGtau_,
			    int PreciseP_, double minDeltat_, bool SampleSusc_)
  {
    my_rank = rank;
    mpi_size = size;
    mu = mu_; U = U_; beta=beta_;
    max_size = max_size_;
    N_flavors = N_flavors_;
    N_ifl = N_ifl_;
    tsample = tsample_;
    warmup = warmup_;
    CleanUpdate = CleanUpdate_;
    minM = minM_;
    minD = minD_;
    PChangeOrder = PChangeOrder_;
    Ncout = Ncout_;
    Naver = Naver_;
    TwoKinks = TwoKinks_;
    GlobalFlip = GlobalFlip_;
    treshold = treshold_;
    SampleGtau = SampleGtau_;
    PreciseP = PreciseP_;
    minDeltat = minDeltat_;
    SampleSusc = SampleSusc_;
  }
};

inline double linear(double x, const pair<double,double>& x1, const pair<double,double>& x2)
{
  return x2.second + (x1.second-x2.second)*(x-x2.first)/(x1.first-x2.first);
}
inline double quadr(double x, const pair<double,double>& x1, const pair<double,double>& x2, const pair<double,double>& x3)
{
  return (x-x2.first)*(x-x3.first)*x1.second/((x1.first-x2.first)*(x1.first-x3.first))+
    (x-x1.first)*(x-x3.first)*x2.second/((x2.first-x1.first)*(x2.first-x3.first))+
    (x-x1.first)*(x-x2.first)*x3.second/((x3.first-x1.first)*(x3.first-x2.first));
}
inline double expin(double x, const pair<double,double>& x1, const pair<double,double>& x2)
{
  if (x2.second*x1.second<0) return x1.second + (x2.second-x1.second)*(x-x1.first)/(x2.first-x1.first);
  double beta = log(x2.second/x1.second)/(x1.first-x2.first);
  double alpha = x1.second*exp(beta*x1.first);
  return alpha*exp(-beta*x);
}