/Users/mourrain/Devel/mmx/realroot/include/realroot/solver_uv_bspline.hpp

Go to the documentation of this file.

/*******************************************************************
 *   This file is part of the source code of the realroot kernel.
 *   Author(s): J. Seland SINTEF
 *              B. Mourrain GALAAD, INRIA
 ********************************************************************/
#ifndef realroot_solver_bspline_hpp
#define realroot_solver_bspline_hpp
/********************************************************************/
#include <realroot/univariate.hpp>
#include <realroot/polynomial_bernstein.hpp>
#include <realroot/solver.hpp>
//--------------------------------------------------------------------
#define NO_ROOT -1.0
#define TMPL template<class Real>
namespace mmx {
  
struct Bspline{};

//-----------------------------------------------------------------------------
TMPL
Real knotSum(Real* t, int k, int d)
{
  Real ret = t[k+1];
  for ( int i=2; i<=d; ++i) ret+=t[k+i];
  return ret;
}

//-----------------------------------------------------------------------------
TMPL
struct solver_bspline {
  Real *m_t, *m_c;
  unsigned m_n, maxn, m_k ;
  int m_d;
  Real eps;

  solver_bspline(unsigned n, unsigned d, unsigned N=100): m_n(n), maxn(N+n), m_d(d), eps(1e-7) { 
    m_t= new Real[2*n+N]; 
    for(unsigned i=0; i< 2*n+N; i++) m_t[i]=Real(1);
    m_c = new Real[n+N];
    for(unsigned i=0; i< n+N; i++) m_c[i]=Real(0);
  }
  
  int  insert_knot ( const Real& x, int mu) ;
  Real first_root  ();

};

//-----------------------------------------------------------------------------
// insert knot x in interval mu by Boehms algorithm
// Note: t,c must have capacity at least n+1, n+d+2 respectively
// require that x>=t[mu]
TMPL
int solver_bspline<Real>::insert_knot( const Real& x, int mu) 
{
  mu = std::max(mu,m_d);
  while ( x>=m_t[mu+1]) mu++;
  
  for ( int i=m_n; i>mu; i--) {
    m_c[i] = m_c[i-1];
  }
  
  for ( int i=mu; i>=mu-m_d+1; i--) {
    const Real alpha = (x-m_t[i])/(m_t[i+m_d]-m_t[i]);
    m_c[i] = (1.0-alpha) * m_c[i-1] + alpha * m_c[i];
  }

  m_t[m_n+m_d] = 1.0;
  
  for ( int i=m_n; i>mu+1; --i)
    m_t[i] = m_t[i-1];
  m_t[mu+1] = x;
  
  //    n++;
  return mu;
}
  
// --------------------------------------------------------------------------
TMPL
Real solver_bspline<Real>::first_root(void)
{

  unsigned k = 1;
  Real x = NO_ROOT;
    
  for ( ;  m_n < maxn; ++m_n ) {
    
    while( k<m_n && (m_c[k-1]*m_c[k] > 0.0) )
      ++k;
    
    
    // Off end?
    if ( k >= m_n ) {
      //if (x!= NO_ROOT)
      //x = -x; // NO_ROOT;
      x = NO_ROOT;
      break;
    }
    
    // Interval converged?
    const Real diff = m_t[k+m_d]-m_t[k];

    if (diff<eps) {
      x = m_t[k];
      break;
    }
    
    const Real av = knotSum(m_t,k-1, m_d);
    const Real lambda = m_c[k-1]/(m_c[k-1]-m_c[k]);
    x = (av + lambda*diff)/Real(m_d);

    // Stopping criterion
    const Real e = std::max(x, m_t[k+m_d-1]) - std::min(x, m_t[k+1] );
    // const Real e = max(fabs(x-t[k+1]),fabs(x-t[k+m_d-1]));
    if (e < eps)
      break;
    
    // Refine spline
    insert_knot(x,k);
  }
  
  //  dx = d*(m_c[k]-m_c[k-1])/(t[k+m_d]-t[k]);
  m_k  = k;
  return x;
}
//--------------------------------------------------------------------
template<class Ring>
struct solver<Ring, Bspline> {
  typedef double                                    solution_t;
  typedef solver_bspline<typename Ring::Scalar>     base_t;

  template<class POL,class T> static solution_t
  first_root(const POL& p, const T& u, const T& v);

  template<class POL> static solution_t
  first_root(const POL& p);

};
//--------------------------------------------------------------------
template<class Ring> 
template<class POL, class T> 
typename solver<Ring,Bspline>::solution_t
solver<Ring,Bspline>::first_root(const POL& p,
                                 const T& u, const T& v) {

  typedef typename Ring::scalar_t                     Real;
  typedef typename POL::Ring                          ring_t;
  typedef typename solver<Ring, Bspline>::base_t      base_t;
  typedef GMP::rational                               rational;
  
  rational U,V;
  let::assign(U,u);
  let::assign(V,v);
  int sz= degree(p)+1;

  std::vector<rational> bp(sz);
  univariate::convertm2b(bp,p,sz,U,V);

  solver_bspline<Real> slv(sz,sz-1);

  for ( int i = 0; i < sz; i ++ ) {
    slv.m_c[i]=as<Real>(bp[i]);
    slv.m_t[i]=Real(0);
  }

  return u +(v-u)*slv.first_root();

}

template<class Ring> 
template<class POL> 
typename solver<Ring,Bspline>::solution_t
solver<Ring,Bspline>::first_root(const POL& p) {

  typedef typename Ring::Scalar                       Real;
  typedef typename POL::Ring                          ring_t;
  typedef typename solver<Ring, Bspline>::base_t      base_t;
  typedef GMP::integer                                integer;
  typedef GMP::rational                               rational;


  int d= degree(p), sz=d+1;
  solver_bspline<Real> slv(sz,sz-1);

  std::vector<rational> bp(sz);
  binomials<integer> bn;
  rational c;
  if(p[0]==0) 
    return 0;
  else if(p[0]<0) {
    numerics::rounding<Real> rnd( numerics::rnd_up() );
    for(int i=0;i<sz;i++) {
      c= as<rational>(p[i])/bn(d,i);
      slv.m_c[i]=as<Real>(c);
    }
  } else {
    numerics::rounding<Real> rnd( numerics::rnd_dw() );
    for(int i=0;i<sz;i++) {
      c= as<rational>(p[i])/bn(d,i);
      slv.m_c[i]=as<Real>(c);
    }
  }
  for ( int i = 0; i < sz; i ++ ) {
    slv.m_t[i]=Real(0);
  }

  Real x= slv.first_root();
 
  return x/(1-x);

}
//====================================================================
} //namespace mmx
#undef TMPL
#undef NO_ROOT
#endif  //realroot_solver_bspline_hpp