template<class Rand>
  class ContinuousRandWithLagrangeInterpolation {
  public:
    int num;
    double x_max;
    Rand rand;
    std::vector<double> x, f, den;
    ContinuousRandWithLagrangeInterpolation(double _x_max, int _num):
      rand(), x_max(_x_max), num(_num),
      x(_num+1,0), f(_num+1,0), den(_num+1,0)
    {
      init();
    }
    void init(){
      int i;
      for(i=0;i<num+1;i++){
	x[i] = x_max/num*i;
	f[i] = rand();
      }
      for(i=0;i<num+1;i++){
	den[i] = mult(x[i], i);
      }
    }
    double mult(const double& y, const int& j){
      int i;
      double m=1;
      for(i=0;i<num+1;i++){
	if(i != j){
	  m *= ( y - x[i] );
	}
      }
      return m;
    }
    double operator()(const double& y) {
      double p=0;
      int i;
      for(i=0;i<num+1;i++){
	p += f[i] * mult(y, i) / den[i];
      }
      return p;
    }
  };