cell_mv_bernstein.hpp

Go to the documentation of this file.

/*****************************************************************************
 * M a t h e m a g i x
 *****************************************************************************
 * cell for bernstein systems
 * 2010-02-16
 *****************************************************************************
 *               Copyright (C) 2010 INRIA Sophia-Antipolis
 *****************************************************************************
 * Comments :
 ****************************************************************************/
# ifndef realroot_cell_mv_bernstein_hpp
# define realroot_cell_mv_bernstein_hpp

# include <realroot/Seq.hpp>
# include <realroot/sign_variation.hpp>
# include <realroot/Interval.hpp>
# include <realroot/polynomial_bernstein.hpp>
# include <realroot/polynomial_sparse.hpp>

# undef C
# define TMPL template<typename C>
# define SELF cell_mv_bernstein<C>

//====================================================================
namespace mmx {

 
TMPL 
class cell_mv_bernstein {
  
public:

  typedef C Scalar;
  typedef polynomial<C, with<Bernstein> > Polynomial;

  cell_mv_bernstein(void) ;
  cell_mv_bernstein(const Polynomial& p) {
    m_doms<<Interval<double>(0,1);
    m_equs<<p;
  }
  cell_mv_bernstein(const Seq<Polynomial>& E) {
    for(unsigned i=0;i<E.size();i++) m_doms<<Interval<double>(0,1);
    for(unsigned i=0;i<E.size();i++) m_equs<<E[i];
  }

  cell_mv_bernstein(const Seq<Polynomial>& E, const Seq<Interval<double> >& D) {
    for(unsigned i=0;i<D.size();i++) m_doms<<D[i];
    for(unsigned i=0;i<E.size();i++) m_equs<<E[i];
  }

  template<class POL>
  cell_mv_bernstein(const Seq<POL>& E, const Seq<Interval<double> >& D) {
    Seq<double> dmn;
    for(unsigned i=0;i<D.size();i++) {
      m_doms<<D[i];
      dmn <<D[i].lower() <<D[i].upper();
    }

    for(unsigned i=0;i<E.size();i++) {
      Polynomial p;
      let::assign(p,E[i],dmn);
      m_equs<<p;
    }
  }

  ~cell_mv_bernstein(void) ;

  unsigned nbeq() const {return m_equs.size();}
  unsigned nbeq()       {return m_equs.size();}

  Polynomial  equation(unsigned i) const {return m_equs[i];}
  Polynomial& equation(unsigned i)       {return m_equs[i];}

  unsigned nbvar() const {return m_doms.size();}
  unsigned nbvar()       {return m_doms.size();}

  Seq<Interval<double> >  domain() const {return m_doms;}
  Seq<Interval<double> >& domain()       {return m_doms;}

  Interval<double>  domain(unsigned i) const {return m_doms[i];}
  Interval<double>& domain(unsigned i)       {return m_doms[i];}

  C size(void) const {
    C s=0;
    for(unsigned i=0;i<m_doms.size();i++)
      s=std::max(s,m_doms[i].upper()-m_doms[i].lower());
    return s;
  }


private:
  Seq<Polynomial>         m_equs;
  Seq<Interval<double> >  m_doms;
};

//--------------------------------------------------------------------
TMPL
SELF::cell_mv_bernstein(void)  {}

TMPL
SELF::~cell_mv_bernstein(void) {}

//--------------------------------------------------------------------
struct binary_approx {
  
  static double m_eps;

  template<class Cell, class Stack>
  static void subdivide(Cell* cl, Stack* stack);

  template<class Cell>
  static bool reduce(Cell* cl);

  template<class Cell>
  static bool regular(Cell* cl);

} ;

double binary_approx::m_eps=1e-6;

template<class Cell, class Stack> void 
binary_approx::subdivide(Cell* cl, Stack* st) {
  //std::cout<<"Subdivide "<<cl->equation(0)<< " "<<cl->domain(0)<<std::endl;

  Cell* left  = new Cell(*cl);
  Cell* right = new Cell(*cl);

  unsigned v=0; 
  typename Cell::Scalar s=cl->domain(0).upper()-cl->domain(0).lower(),s0;
  for (unsigned i=0;i<cl->nbvar();i++)
    if((s0=cl->domain(i).upper()-cl->domain(i).lower())>s) {
    s=s0;v=i;
    }
  typename Cell::Scalar m=(cl->domain(v).upper()+cl->domain(v).lower())/2;

  for (unsigned i=0;i<cl->nbeq();i++)
    tensor::split(left->equation(i), right->equation(i), v);
  left->domain(v).upper()=m;
  right->domain(v).lower()=m;

  //  std::cout<<"==> "<<left->equation(0)<< " "<<left->domain(0)<<std::endl;
  //  std::cout<<"==> "<<right->equation(0)<< " "<<right->domain(0)<<std::endl;
  //  std::cout<<std::endl;

  st->push(left); 
  st->push(right);
};

template<class Cell> bool
binary_approx::reduce(Cell* cl) {  
  if(cl->size() < m_eps) return false;
  for(unsigned i=0;i<cl->nbeq();i++)
    if(!has_sign_variation(cl->equation(i))) return false;
  return true; 
}

template<class Cell> bool
binary_approx::regular(Cell* cl) { 
  for(unsigned i=0;i<cl->nbeq();i++)
    if(!has_sign_variation(cl->equation(i))) return false;
  return true;
}

//--------------------------------------------------------------------
struct binary_isolate: binary_approx {
  
  template<class Cell>
  static bool reduce(Cell* cl);

} ;

template<class Cell> bool
binary_isolate::reduce(Cell* cl) {  
  if(cl->size() < m_eps) return false;
  std::cout<<"isolate"<<std::endl;
  for(unsigned i=0;i<cl->nbeq();i++)
    if(!has_sign_variation(cl->equation(i))) return false;
  
  
  return true; 
}

//====================================================================
} // namespace mmx
//====================================================================
# undef TMPL
# undef SELF
# endif