namespace GENERIC{ template<class R, class C>
		   void domain(  C & a, C & b, const R & r  ) 
		   {
		     a = C(0);
		     b = C(1);
		   }  
}

template<class R, class C> inline
void domain( C & a, C & b, const V<R>& c )
{
  using namespace GENERIC;
  domain(a,b,rep(c));
};


template<class R, class Point3>  inline
void eval( Point3& p, const V<R> & c, const typename Point3::value_type& t )
{
  eval(p, rep(c),  t );
};

namespace GENERIC { template<class R, class Point3> inline
		    typename Point3::value_type  project( const R & c, const Point3& query )
		    {
		      std::cout << " Point / Parametric Curve Projection TO DO ! \n";
		      return typename Point3::value_type(0);
		    };
};

template<class R, class Point3>
typename Point3::value_type project( const R & c, const Point3& query )
{
  using GENERIC::project; 
  return project(c,query);
};


namespace GENERIC{ template<class R, class Point3> inline
		   void project(  Point3 & result, const R& c, const Point3 & query )
		   {
		     eval(view(c), result, project(view(c),query) );
		   }; 
};

template<class R, class Point3> inline 
void project(  Point3 & result, const V<R> & c, const Point3 & query ) 
{ 
  using GENERIC::project; 
  project(result,rep(c),query); 
};

 

namespace GENERIC{
  template<class R, class Point3>
  void eval( Point3 & p, Point3 & dp, const R & r, const typename Point3::value_type & t )
  {
    typedef typename Point3::value_type C;
    C a,b;
    domain(a,b,view(r));
    C h = (b-a)/C(1000);
    a = t - h;
    b = t + h;
    eval(p,view(r),a);
    eval(dp,view(r),b);
    for ( int i = 0; i < 3; i ++ )
      dp[i] = (dp[i]- p[i])/(C(2)*h);
    eval(p,view(r),t);
  };
};

template<class R, class Point3>  inline
void eval( Point3 & p, Point3 & dp, const V<R> & c, const typename Point3::value_type & t )
{
  eval(p,dp,rep(c),t);
};


namespace GENERIC { template<class R, class Point3> inline
		    void sample( Point3 * smp, const R & c, const typename Point3::value_type * t, int n )
		    {
		      for ( Point3 * o = smp; o != smp + n; eval(*o++,view(c),*t++) );
		    };
};

template<class R, class Point3> 
void sample( Point3 * smp, const V<R> & c, const typename Point3::value_type * tvals, int n ) 
{
  using namespace GENERIC;
  sample(smp,rep(c),tvals,n);
};

namespace GENERIC {template<class R, class C> inline
		   void tsample( C * t, const R & r, int n )
		   {
		     C a,b;
		     domain(a,b,view(r));
		     C h = (b-a)/C(n);
		     for ( C * p = t; p != t + n; *t++ = a, a += h );
		   };
};

template<class R, class C> inline
void tsample( C * t, const V<R> & c, int n )
{
  using namespace GENERIC;
  tsample(t,rep(c),n);
};

namespace GENERIC {template<class R, class Point3> inline
		   void sample( Point3 * smp, typename Point3::value_type * t, const R& r, int n )
		   {
		     tsample(t,view(r),n);
		     sample(smp,view(r),t,n);
		   };
};
  
template<class R, class Point3> 
void sample( Point3 * smp, typename Point3::value_type * t, const V<R> & c, int n )
{
  using namespace GENERIC;
  sample(smp,t,rep(c),n);
};


namespace GENERIC { template<class PointArray, class CoefficientArray, class R>
		    void eps_sample( PointArray& psmp, CoefficientArray& tsmp, const R& crv, const typename CoefficientArray::value_type & eps )
		    {
		      typedef typename CoefficientArray::value_type inputc;
		      typedef typename PointArray::value_type       point;
		      
		      point p,dp,dfp;
		      inputc t, tmax, tprev;
		      
		      domain(t,tmax,crv);
		      
		      while ( t < tmax )
			{
			  eval( p, dfp, crv, t );
			  psmp.push_back(p);
			  tsmp.push_back(t);
			  t = t + eps/sqrt(dotprod(dfp,dfp));
			  do
			    {
			      eval(p,dfp,crv,t);
			      dp = p-psmp.back();
			      f = dotprod(dp,dp)-eps*eps;
			      tprev = t;
			      t = tprev - f/(2*dotprod(dp,dfp));
			    }
			  while( !(std::abs(f) < 1e-12 || std::abs(tprev-t) < 1e-8) );
			};
		      eval(p,crv,tmax);
		      psmp.push_back(p);
		      tsmp.push_back(tmax);
		    };
}

template<class PointArray, class CoefficientArray, class R> inline
void eps_sample( PointArray& p, CoefficientArray& tsmp, const R & crv, const typename CoefficientArray::value_type & eps )
{
  using namespace GENERIC;
  eps_sample(p,t,rep(crv),eps);
};

template<class C, int n = 3>
struct Sample 
{
  typedef fxv<C,n> point_t;
  std::vector<point_t> m_cvals;
  std::vector<C>       m_tvals;
  point_t& operator[]( int i ) { return m_cvals[i]; };
  const point_t& operator[]( int i ) const { return m_cvals[i]; };
  const C& tval( int i ) const { return m_tvals[i]; };
  C & tval( int i ) { return m_tvals[i]; };
  template<class R>
  sample( const V<R>& c, int n  ) : m_cvals(n), m_tvals(n) 
  { 
    tsample(&m_tvals[0],c,n);
    psample(&m_cvals[0],c,&m_tvals[0],n);
  };
  unsigned size() const { return m_cvals.size(); };
};


namespace GENERIC {template<class PointArray, class CoefficientArray, class R>
		   void msample( PointArray & p, CoefficientArray & c, const R & r, int n )
		   {
		     p.resize(n);
		     c.resize(n);
		     sample( &p[0], &c[0], view(r), n );
		   };
};

template<class PointArray, class CoefficientArray, class R> 
void msample( PointArray & p, CoefficientArray & c, const V<R>& c, int n ) 
{
  using namespace GENERIC;
  msample(p,c,rep(c),n);
};

template<class R, class C, int n = 3>
struct MSegmentTree : Sample<C,n> 
{

  typedef unsigned short mcode_t;
  std::vector< mcode_t > m_mrks;
  struct sgnode 
  {
    int a, b ; // sample span
    sgnode * l, * r;
  };
  
  sgnode * m_base;
  
  template<class R> 
  MSegmentTree( const V<R> & c, int nsmp )
  {
    msample( m_cvals, m_tvals, c, n ); // query for monotonic sample
    mcode_t cd = code( smp(0), smp(1) );
    m_mrks.push_back( mrk(tval(0),0,cd) );
    for ( int i = 1; i < m_smp->size()-1; i ++ )
      {
	mcode_t tmp = code( smp(i), smp(i+1) );
	if ( tmp != cd ) { 
	  cd = tmp;
	  m_mrks.push_back( mrk(tval(i), i,cd) );
	};
	};
      m_mrks.push_back( mrk(tval(m_smp->size()-1), m_smp->size()-1, cd) );
  };
};

