include/shape/algebraic_curve_fcts.hpp

Go to the documentation of this file.

/*****************************************************************************
 * M a t h e m a g i x
 *****************************************************************************
 * AlgebraicCurve
 * Bernard Mourrain
 *****************************************************************************
 *               Copyright (C) 2009 INRIA Sophia-Antipolis
 *****************************************************************************
 * Comments :
 ****************************************************************************/
# ifndef algebraic_curve_fcts_hpp
# define algebraic_curve_fcts_hpp

# include <realroot/Seq.hpp>
# include <realroot/ring_monomial_tensor.hpp>
# include <realroot/ring_bernstein_tensor.hpp>
# include <realroot/solver_uv_continued_fraction.hpp>
# include <shape/point.hpp>
# include <shape/vertex.hpp>
# include <shape/solver_implicit.hpp>
# include <shape/algebraic_curve.hpp>

# define TMPL template<class C,class V>
# define SELF algebraic_curve<C,V>
# define PointSet point_set<C,3,V>
# define BoundingBox bounding_box<C,V>
//====================================================================
namespace mmx {
namespace shape {

TMPL Seq<typename topology<C,V>::Point *>
extremal(const SELF& c, const BoundingBox& b)  {
  typedef typename topology<C,V>::Point Point;
  Seq<Point*> res;
  solver_implicit<C,V>::extremal(res,c.equation(),b);
  return res;
}

TMPL PointSet
plot(const SELF& c, const BoundingBox& bx, int N= 500)  {
  
  typedef double                      Scalar;
  typedef typename PointSet::Point    Point;
  typedef typename SELF::Polynomial   Polynomial;
  typedef typename Polynomial::Scalar Rational;
  typedef solver<Rational, ContFrac<Approximate> > Solver; 

  Point pt(0,0,0);
  PointSet pts;

  Rational a0(bx.xmin()), a1(bx.xmax()), b0(bx.ymin()), b1(bx.ymax());   
  Rational dx= (a1-a0)/N;
  Seq<Polynomial> sy = coefficients(c.equation(),1);
  int d = sy.size()-1;
  polynomial<Rational, with<MonomialTensor> > P(0,d);

  for (int i=0; i<N;i++,a0+=dx) {
    pt[0] = as<Scalar>(a0);
    for (int u=0;u<=d;u++) P[u] = sy[u](a0);
    typename Solver::Solutions sol; Solver::solve(sol,P,b0,b1);
    for(unsigned u=0;u<sol.size();u++) {
      pt[1] = as<Scalar>((lower(sol[u])+upper(sol[u]))/2);
      pts<<pt;    
    }
  }
  Rational dy= (b1-b0)/N;
  a0= Rational(bx.xmin());
  Seq<Polynomial> sx = coefficients(c.equation(),0);
  d = sx.size()-1;
  polynomial<Rational, with<MonomialTensor> > Q(0,d);
  for (int i=0; i<N;i++,b0+=dy) {
    pt[1] = as<Scalar>(b0);
    for (int u=0;u<=d;u++) Q[u] = sx[u]((Rational)0,b0);
    typename Solver::Solutions sol; Solver::solve(sol,Q,a0,a1);
    for(unsigned u=0;u<sol.size();u++) {
      pt[0] = as<Scalar>((lower(sol[u])+upper(sol[u]))/2);
      pts<<pt;    
    }
  }

  return pts;
}
//====================================================================
  } ; // namespace mmx
} ; // namespace mmx
//====================================================================
# undef TMPL
# undef SELF
# undef PointSet
# undef BoundingBox
# endif // shape_algebraic_curve_fcts_hpp

---