/Users/mourrain/Devel/mmx/shape/include/shape/cell2d_algebraic_curve.hpp

Go to the documentation of this file.

/*****************************************************************************
 * M a t h e m a g i x
 *****************************************************************************
 * cell2d_algebraic_curve
 * 2008-03-20
 * Julien Wintz & Bernard Mourrain
 *****************************************************************************
 *               Copyright (C) 2008 INRIA Sophia-Antipolis
 *****************************************************************************
 * Comments :
 ****************************************************************************/

# ifndef shape_cell2d_algebraic_curve_hpp
# define shape_cell2d_algebraic_curve_hpp

# include <realroot/Interval.hpp>
# include <realroot/polynomial_bernstein.hpp>
# include <realroot/Seq.hpp>
# include <shape/point.hpp>
# include <shape/solver_implicit.hpp>
# include <shape/cell2d.hpp>

# include <stack>
# define  STACK std::stack<Point*>
# define TMPL template<class C, class V>
# define TMPL1 template<class V>
# define SELF cell2d_algebraic_curve<C,V>
# define REF REF_OF(V)

# undef Cell_t
//====================================================================
namespace mmx {
namespace shape {
//====================================================================
// TMPL struct cell2d_algebraic_curve;

// TMPL struct with_cell2d_algebraic_curve { 
//   typedef cell2d_algebraic_curve<K> Cell2dAlgebraicCurve; 
// };
//--------------------------------------------------------------------
TMPL struct cell2d_subdivisor
{
  //  typedef polynomial< Interval<double>, with<Bernstein> > Polynomial;
  //  typedef polynomial< double, with<Bernstein> > Polynomial;
  //  typedef polynomial< mmx::GMP::rational, with<Bernstein> > Polynomial;
  //  typedef solver_implicit<double,with_idx<double> > Solver;
  
  template<class CELL> static void 
  subdivide(CELL & cl, CELL*&, CELL*&, int v, double s);
};

//--------------------------------------------------------------------
TMPL
template<class CELL> void 
cell2d_subdivisor<C,V>::subdivide(CELL& cl, CELL*& left, CELL*& right, int v, double s) {

  typedef typename topology<C,V>::Point Point;
  typedef solver_implicit<C,V> Solver;


  double eps =Approximate().eps;
  if(v==0) {//direction v=0
       left = new CELL(cl);   left ->set_xmax(s);
       right= new CELL(cl);   right->set_xmin(s);

       //      left = new CELL(cl.m_polynomial, BoundingBox(cl.xmin(),s, cl.ymin(), cl.ymax(), cl.zmin(), cl.zmax()));
       //      right= new CELL(cl.m_polynomial, BoundingBox(s,cl.xmax(), cl.ymin(), cl.ymax(), cl.zmin(), cl.zmax()));

      tensor::split(left->m_polynomial, right->m_polynomial, v);
      Solver::edge_point(left->e_intersections, 
                         left->m_polynomial, Solver::east_edge, left->boundingBox());   
      right->w_intersections=left->e_intersections;
      left ->w_intersections=cl.w_intersections;
      right->e_intersections=cl.e_intersections;
      foreach(Point* p,cl.n_intersections) {
          if(p->x() < s) 
              left ->n_intersections << p ;
          else
              right->n_intersections << p ;
      }
      foreach(Point* p,cl.s_intersections) {
          if(p->x() < s) 
              left ->s_intersections<< p ;
          else
              right->s_intersections<< p ;
      }

      foreach(Point* p,cl.m_singular) {

        if(p->x() < s+eps) 
          left ->m_singular << p ;
        //      else 
        if(p->x() > s-eps) 
          right->m_singular << p ;
      }
      
      /*  Update neighbors  */
      cl.connect0(left, right);
      left->join0(right);


    } else {//direction v==1
        
      left = new CELL(cl);  left->set_ymax(s);
      right= new CELL(cl);  right->set_ymin(s);
      
      //left = new CELL(cl.m_polynomial, BoundingBox(cl.xmin(),cl.xmax(), cl.ymin(), s, cl.zmin(), cl.zmax()));
      //right= new CELL(cl.m_polynomial, BoundingBox(cl.xmin(),cl.xmax(), s, cl.ymax(), cl.zmin(), cl.zmax()));
      
      tensor::split(left->m_polynomial, right->m_polynomial, v);
      Solver::edge_point(left->n_intersections, 
                         left->m_polynomial, Solver::north_edge, left->boundingBox());   
      right->s_intersections=left->n_intersections;
      left ->s_intersections=cl.s_intersections;
      right->n_intersections=cl.n_intersections;
      foreach(Point* p,cl.w_intersections) {
          if(p->y() < s) 
              left ->w_intersections << p ;
          else
            right->w_intersections << p ;
      }
      
      foreach(Point* p,cl.e_intersections) {
          if(p->y() < s) 
              left ->e_intersections << p ;
          else
              right->e_intersections << p ;
      }
      
      foreach(Point* p,cl.m_singular) {
        if(p->y() < s+eps)
          left->m_singular << p ;
        //      else 
        if(p->y() > s-eps)
          right->m_singular << p ;
      } 
        
      foreach(Point* p,cl.m_xcritical) {
        if(p->y() < s+eps)
          left->m_xcritical << p ;
        //      else 
        if(p->y() > s-eps)
          right->m_xcritical << p ;
      } 

      foreach(Point* p,cl.m_ycritical) {
        if(p->y() < s+eps)
          left->m_ycritical << p ;
        //      else 
        if(p->y() > s-eps)
          right->m_ycritical << p ;
      } 
      
      /*  Update neighbors  */
      cl.connect1(left, right);
      left->join1(right);

    }   
    /* disconnect parent */
    cl.disconnect( );

}
//====================================================================
template<class C, class V=default_env>
struct cell2d_algebraic_curve : public cell2d<C,REF> {
public:
  typedef topology<C,REF>    Topology;
  typedef typename Topology::Point          Point;
  //  typedef typename cell2d_algebraic_curve_def<K>::Vertex         Vertex;
  //  typedef typename use<cell2d_algebraic_curve_def,V>::Vector      Vector;
  typedef point<C,REF> Vector;
  typedef typename topology<C,REF>::Edge                          Edge;
  typedef cell<C,REF>           Cell;
  typedef algebraic_set<C,REF>                    AlgebraicSet;
  typedef algebraic_curve<C,REF>                  AlgebraicCurve;
  typedef polynomial< Interval<C>, with<Bernstein> > Polynomial;
  typedef bounding_box<C,REF>                   BoundingBox;

  typedef solver_implicit<C,V>                  Solver;
  typedef topology2d<C,V>                       Topology2d;

  cell2d_algebraic_curve(const Polynomial &, const BoundingBox&, bool inter=true);
  cell2d_algebraic_curve(const char*, const BoundingBox&);
  cell2d_algebraic_curve(const typename AlgebraicSet::Polynomial& , const BoundingBox&);
  cell2d_algebraic_curve(AlgebraicCurve*, const BoundingBox &x);
  cell2d_algebraic_curve(const SELF& cl)
    : cell2d<C,REF>(cl.boundingBox()), m_polynomial(cl.m_polynomial) {}

  Polynomial equation(void) { return  m_polynomial; }
  
  bool  is_active (void) ;
  bool  is_regular(void) ;
  bool  is_intersected(void) {return (this->nb_intersect()>0);};

  virtual Point * pair(Point * p, int & sgn) ;
  virtual Point * starting_point( int sgn) ;

  virtual bool insert_regular (Topology*);// = 0;
  virtual bool insert_singular(Topology*);// = 0;
  virtual void subdivide(Cell*& left, Cell*& right, int v, double s);
  Vector gradient(const Point& p);

public:
  Polynomial        m_polynomial;
  Seq<Point *>      m_xcritical;
  Seq<Point *>      m_ycritical;
//  AlgebraicCurve*   m_shape;
};

//====================================================================
TMPL inline void 
print(SELF* cv) {
  typedef typename SELF::Point Point;
  std::cout << "point(s) b: "; 
  foreach(Point* p, cv->n_intersections)  
    std::cout <<" ["<<p->x()<<" "<<p->y()<<"]";
  foreach(Point* p, cv->s_intersections)  
    std::cout <<" ["<<p->x()<<" "<<p->y()<<"]";
  foreach(Point* p, cv->w_intersections)  
    std::cout <<" ["<<p->x()<<" "<<p->y()<<"]";
  foreach(Point* p, cv->e_intersections)  
    std::cout <<" ["<<p->x()<<" "<<p->y()<<"]";
  std::cout<<std::endl;
  
  std::cout << "Point(s) s: "; 
  foreach(Point* p, cv->m_singular)  
    std::cout <<" ["<<p->x()<<" "<<p->y()<<"]";
  std::cout<<std::endl;
}
//--------------------------------------------------------------------
TMPL
SELF::cell2d_algebraic_curve(const Polynomial& pol, const BoundingBox& bx, bool inter): cell2d<C,REF>(bx), m_polynomial(pol)
{
  //std::cout<<"Alegbraic cell: "<<bx<<std::endl;
  
  if (inter)
  {
    //std::cout<<"compute inter "<<bx<<std::endl;
  Solver::edge_point(this->n_intersections, m_polynomial, Solver::north_edge, bx);
  Solver::edge_point(this->s_intersections, m_polynomial, Solver::south_edge, bx);
  Solver::edge_point(this->w_intersections, m_polynomial, Solver::west_edge , bx);
  Solver::edge_point(this->e_intersections, m_polynomial, Solver::east_edge , bx);
  }
    //Solver::extremal(this->m_singular, m_polynomial, bx); 

}
//--------------------------------------------------------------------  
TMPL
SELF::cell2d_algebraic_curve(const typename AlgebraicSet::Polynomial& p, const BoundingBox & b): cell2d<C,REF>(b)
{
  //std::cout<<"Alegbraic cell: "<<p<<std::endl;
    Seq<mmx::GMP::rational> bx;
    bx<<as<mmx::GMP::rational>(b.xmin());
    bx<<as<mmx::GMP::rational>(b.xmax());
    bx<<as<mmx::GMP::rational>(b.ymin());
    bx<<as<mmx::GMP::rational>(b.ymax());
    
      tensor::bernstein<mmx::GMP::rational> polq(p.rep(),bx);

      let::assign(m_polynomial.rep(),polq);

    Solver::edge_point(this->n_intersections, m_polynomial, Solver::north_edge, b);
    Solver::edge_point(this->s_intersections, m_polynomial, Solver::south_edge, b);
    Solver::edge_point(this->w_intersections, m_polynomial, Solver::west_edge , b);
    Solver::edge_point(this->e_intersections, m_polynomial, Solver::east_edge , b);

    Solver::extremal(this->m_singular, p, b); 

    // foreach(Point* p,this-> s_intersections) 
    //   std::cout<<"s_ boundary point: "<<p->x()<<", "<<p->y() <<std::endl;
    // foreach(Point* p,this-> e_intersections) 
    //   std::cout<<"e_ boundary point: "<<p->x()<<", "<<p->y() <<std::endl;
    // foreach(Point* p,this-> n_intersections) 
    //   std::cout<<"n_ boundary point: "<<p->x()<<", "<<p->y() <<std::endl;
    // foreach(Point* p,this-> w_intersections) 
    //   std::cout<<"w_ boundary point: "<<p->x()<<", "<<p->y() <<std::endl;

    // foreach(Point* p,this-> m_singular) 
    //   std::cout<<"extremal point: "<<p->x()<<", "<<p->y() <<std::endl;

    //std::cout<<"Found "<<this->nb_intersect()<<" boundary and "<<this->m_singular.size()<<" extremal points" <<std::endl;
  }
//--------------------------------------------------------------------
TMPL
SELF::cell2d_algebraic_curve(const char* s, const BoundingBox & b): cell2d<C,REF>(b)
{
    Seq<mmx::GMP::rational> bx;
    bx<<as<mmx::GMP::rational>(b.xmin());
    bx<<as<mmx::GMP::rational>(b.xmax());
    bx<<as<mmx::GMP::rational>(b.ymin());
    bx<<as<mmx::GMP::rational>(b.ymax());
    
    typename AlgebraicSet::Polynomial pol(s);
       tensor::bernstein<mmx::GMP::rational> polq(s,bx);

       let::assign(m_polynomial.rep(),polq);

    Solver::edge_point(this->n_intersections, m_polynomial, Solver::north_edge, b);
    Solver::edge_point(this->s_intersections, m_polynomial, Solver::south_edge, b);
    Solver::edge_point(this->w_intersections, m_polynomial, Solver::west_edge , b);
    Solver::edge_point(this->e_intersections, m_polynomial, Solver::east_edge , b);

    Solver::extremal(this->m_singular, pol, b); 
    foreach(Point* p,this-> m_singular) 
      std::cout<<"**extremal point: "<<p->x()<<" "<<p->y() <<std::endl;
}
//--------------------------------------------------------------------  
TMPL
SELF::cell2d_algebraic_curve(AlgebraicCurve* cv, const BoundingBox & b)
{
  new (this) SELF(cv->equation(),b );
}
//--------------------------------------------------------------------
TMPL bool 
SELF::is_regular() {
// Assumption: cells with extremal points (dx=dy=0) that do not belong on the curve do not intersect the curve (assumes deep enough subdivision)

    if(this->m_singular.size()>1) return false;

//    if (0)
    if (this->m_singular.size()==0) {
      if(!has_sign_variation(m_polynomial)) return true;

      Polynomial dx= diff(m_polynomial,0);
      if(!has_sign_variation(dx)) return true;

      Polynomial dy= diff(m_polynomial,1);
      if(!has_sign_variation(dy)) return true;

      int n= this->nb_intersect();
      return   ( n==2 );
      //return   ( n%2==0);
      //return   (n == 4 || n == 2 || n == 0 );
    }

    //return true;

    //--------------------------------------------------------------------
    if(this->m_singular.size()==1) {
//    Polynomial 
//      dxf = diff(m_polynomial,0),
//      dyf = diff(m_polynomial,1);
//       std::cout<< "__________________"<<std::endl;
//       std::cout<< m_polynomial;std::cout<< std::endl;
//       std::cout<< "__________________"<<std::endl;
//       std::cout<< dxf;std::cout<< std::endl;
//       std::cout<< "__________________"<<std::endl;
//       std::cout<< dyf;std::cout<< std::endl;
//       std::cout<< "__________________"<<std::endl;

//      int d;
//      d = 1-2*topological_degree(dxf,dyf);
//      d = 0;
      return false; //(d==ni);
    }    
    
    return true;
}
//--------------------------------------------------------------------
TMPL bool 
SELF::is_active() {

    return (has_sign_variation(this->m_polynomial));

//   if(this->m_singular.size()>0) 
//       return true;
//     return (this->nb_intersect()>0);
}
//--------------------------------------------------------------------
TMPL bool 
SELF::insert_regular(Topology* t) {
    Topology2d* s= dynamic_cast<Topology2d*>(t);
    Seq<Point*> l;
    int ni, sgn(1);
    Point *p,*q;
    l= this->intersections();
    ni=l.size();

    std::cout<<"alg. curve, inserting regular "<< *this <<std::endl;
    foreach( Point* e, l)
      std::cout<<"("<< e->x()<<", "<<e->y()<<")"<<std::endl;

    //assert(ni%2==0);

    if(this->m_singular.size()==0) 
    {
      foreach(Point* v, l)
        s->insert(v);


      while (l.size()>0)
      {
        //std::cout<<"L = "<< l <<std::endl;      
        p= l[0];
        q= this->pair(p,sgn);
        l.erase(0);
        //std::cout<<"erased "<< p <<" "<<std::endl;int k=l.search(q);
        l.erase( l.search(q) );
        //std::cout<<"erased "<< q <<" ~"<<k<<std::endl;
        s->insert(p,q);
        //std::cout<<"edge "<<0<<" "<<q<<std::endl;
        //std::cout<<"l.size()="<<l.size() <<std::endl;
      }
    }
    else
    {
      //std::cout<<"Singular:"<<std::endl;    
      //std::cout<<this<<std::endl;    
      t->insert((BoundingBox*)this,false);
      
      Point *c=this->m_singular[0];
      s->m_specials<<c;
      s->insert(c);
      foreach(Point* p, this->n_intersections) s->insert(new Edge(p,c));
      foreach(Point* p, this->e_intersections) s->insert(new Edge(p,c));
      foreach(Point* p, this->s_intersections) s->insert(new Edge(p,c));
      foreach(Point* p, this->w_intersections) s->insert(new Edge(p,c));
    }

    //std::cout<<"ok "<<std::endl;

    return true;

}
//--------------------------------------------------------------------
TMPL bool 
SELF::insert_singular(Topology* t) {

//  std::cout<<"Singular!!"<<std::endl;    

    //Mark singular box with a cross
    //t->insert((BoundingBox*)this,true);return true;
    
    Topology2d* s = dynamic_cast<Topology2d*>(t);
    int sgn(1);

    Seq<Point*> l;
    foreach(Point* p, this->n_intersections) {
        s->insert(p); l<<p;  
    }
    foreach(Point* p, this->e_intersections) {
      s->insert(p); l<<p;
    }
    foreach(Point* p, this->s_intersections) {
      s->insert(p); l<<p;
    }
    foreach(Point* p, this->w_intersections) {
      s->insert(p); l<<p;
    }
    int ni=l.size();
    
    if(this->m_singular.size()==1) {

        //std::cout<<"Singular:"<<std::endl;    
        //std::cout<<this<<std::endl;    
        //t->insert((BoundingBox*)this,false);
        

      Point *c=this->m_singular[0];
      s->m_specials<<c;
      s->insert(c);
      foreach(Point* p, this->n_intersections) s->insert(new Edge(p,c));
      foreach(Point* p, this->e_intersections) s->insert(new Edge(p,c));
      foreach(Point* p, this->s_intersections) s->insert(new Edge(p,c));
      foreach(Point* p, this->w_intersections) s->insert(new Edge(p,c));
      return true;
    }
    else   // no singular points in the cell
    {
        if(ni==2) { // 2 points on the boundary
            s->insert(l[0],l[1]);
            return true;
        }
        else // more than one branch in the cell
        {
            Point *p,*q;
            while (l.size()>0)
            {
                p= l[0];
                q= this->pair(p,sgn);
                l.erase(0);
                l.erase( l.search(q) );
                s->insert(p,q);
            }
        }
    }
    
    //std::cout<<"singular cell with "<<ni<<" bpts "<<this->m_singular.size()<<" epts"<<std::endl;    

    return true;
  }

//--------------------------------------------------------------------
TMPL void
SELF::subdivide(Cell*& Left, Cell*& Right, int v, double s) {
  cell2d_subdivisor<C,V>::subdivide( *this,(SELF*&)Left,(SELF*&)Right,v,s);
}
//--------------------------------------------------------------------
TMPL  typename SELF::Vector
SELF::gradient(const Point& p) {
  Polynomial 
    dx= diff(m_polynomial,0),
    dy= diff(m_polynomial,1);
  double 
    u0= (p[0]-this->xmin())/(this->xmax()-this->xmin()),
    u1= (p[1]-this->ymin())/(this->ymax()-this->ymin());
  Vector v;
  v[0] = dx(u0,u1);
  v[1] = dy(u0,u1);
  return v;
}

//--------------------------------------------------------------------
TMPL typename SELF::Point*
SELF::starting_point(int sgn) {

//  std::cout<<"startingPoint.AlgCurve"<<std::endl;

  Seq<Point*> l, all;
      unsigned a;
      all = this->intersections();

      // std::cout<<this<<" , "<<a<<std::endl;
      // foreach(Point*p,all)
      //   std::cout<<p->x()<<" "<<p->y()<<" "<<p->z()<<std::endl;


       // if (l.size()==4) // two duplicate corners
       // {
       //   all<< l[1];
       //   all<< l[0];
       // }
       // else
       //   all=l;


      a   = all.size();
      if (a==1) return all[0];

      int ev(0);
      int u ;//position of p
      unsigned
          s=this->s_intersections.size() ,// ~0
          e=this->e_intersections.size() ,// ~1
          //n=this->n_intersections.size() ,// ~2
          w=this->w_intersections.size() ;// ~3

      u= ( 0<s   ? 0 :
         ( 0<s+e ? 1 :
         ( 0<a-w ? 2 :
                   3 )));

      int * sz = this->m_polynomial.rep().szs();
      int * st = this->m_polynomial.rep().str();

          switch (u){
          case 0:
              ev= (this->m_polynomial[0] >0 ? 1:-1);
         break;
          case 1:
              ev= (this->m_polynomial[(sz[0]-1)*st[0]] >0 ? 1:-1);      
         break;
          case 2:
              ev= (this->m_polynomial[sz[0]*sz[1]-1]>0 ? 1:-1);
          break;
          case 3:
              ev= (this->m_polynomial[(sz[1]-1)*st[1]] >0 ? 1:-1);
          break;
          }

          if    (ev*sgn>0)  {  
              return all[0]; 
          }
          else              {  
              return all[1];
          }
}

//--------------------------------------------------------------------
TMPL typename SELF::Point*
SELF::pair( typename SELF::Point* p, int & sgn) {
     // Pair: returns a neighboring point on the component (p,sgn)
     // The candidates are the 2 neighbor intersections of p in *this

      Seq<Point*> all;
      int a;
      all= this->intersections();
      a   = all.size();
      if (a==2)// regular cell with a single branch
         { return (all[0]==p ? all[1]: all[0]);  }
      if (a==1)
      { std::cout<<"Only 1 intersection point in "<< *this <<" (i.e."<<*all[0] <<")"<<std::endl;

//  if  ( abs( l[0]->x()-this->xmin()<EPSILON) && 
//        abs( l[0]->y()-this->ymin()<EPSILON) )


          return (p); }


      int
          s=this->s_intersections.size() ,// ~0
          e=this->e_intersections.size() ,// ~1
          //n=this->n_intersections.size() ,// ~2
          w=this->w_intersections.size() ;// ~3
      int j,k,i= all.search(p);        
    
      //if (this->is_regular() )    { 
      if (this->m_singular.size()==0  )    { 


     Vector grq, grp = this->gradient(*p);       
        foreach(Point* q, all)
            if ( q != p) 
            {
                grq = this->gradient(*q);
                if(grp[0]*grq[0]>0  && grp[1]*grq[1]>0 )
                    return (q);
            }

        //box has 2 or more "identical" branches

        //p is part of critical branch

       std::cout<<"...maybe pair Trouble"<< this<<std::endl; 

        for (int v=0;v<a;v++)
            for (int w=v+1;w<a;v++)
            {
                grp= this->gradient(*all[v]);
                grq= this->gradient(*all[w]);
                if( grp[0]*grq[0]>0  && grp[1]*grq[1]>0 )
                {
                    for (int u=0; u<a;u++)
                        if ( u!=v && u!=w && all[u]!=p )
                            return (all[u]);
                }                
            }
        
        std::cout<<"empty BOX:("<<this->m_singular.size()<<",a="<<a<<")"<<this <<std::endl;         
        foreach(Point*v,all) {
            Vector gr= this->gradient(*v);
            //std::cout<<"("<<v->x()<<","<<v->y()<<"): ";
            std::cout<<(gr[0]>0?1:-1)<<","<<(gr[1]>0?1:-1) <<std::endl;
        }
        
        //look non-tangent direction
        foreach(Point* q, all)
            if ( q != p) 
            {
                grq = this->gradient(*q);
                if( abs(grp[0]) < 0.017)
                { if (grp[1]*grq[1]>0)
                        return (q);
                }
                else if ( abs(grp[1])< 0.01)
                { if (grp[0]*grq[0]>0)
                        return (q);
                }
            }

        std::cout<<"...pair Trouble"<< this<<std::endl; 



     // foreach(Point* q, all)
     //     if ( q != p) 
     //     {
                  return (all[0]);
     //     }
    }
       else {//singular box

           //std::cout<< "SBOX: " <<this <<std::endl;
           //std::cout<< "      " <<all <<std::endl;
           //for (unsigned y=0;y<4;y++)
           //std::cout<<y <<": ("<<all[y]->x()<<","<<all[y]->y()<<")"<<std::endl;

      // Cell with 1 self-intersection
      int ev(0);
      int u, v;//side of p, ln resp.
      j= ( i!=0   ? i-1 : a-1 );
      k= ( i!=a-1 ? i+1 : 0   );

      //std::cout<<"j,i,k (left,p,right)="<<j <<" "<<i << " "<<k<< std::endl; 

      Point *ln= all[j], 
            *rn= all[k] ;

      //std::cout<<"i= "<< i<<", j=" <<j<< std::endl; 
      // std::cout<<"("<<a<<")"<<s<<", "<<e<<", "<<n<<", "<<w<< std::endl; 

      u= ( i<s   ? 0 :
         ( i<s+e ? 1 :
         ( i<a-w ? 2 :
                   3 )));
      v= ( j<s   ? 0 :
         ( j<s+e ? 1 :
         ( j<a-w ? 2 :
                   3 )));

      int * sz = this->m_polynomial.rep().szs();
      int * st = this->m_polynomial.rep().str();

          switch (u){
          case 0:
              ev= (this->m_polynomial[0] >0 ? 1:-1);
              if (v==0 && j%2==0) // p, ln on the same face
                  ev*=-1;
         break;
          case 1:
              ev= (this->m_polynomial[(sz[0]-1)*st[0]] >0 ? 1:-1);      
              if (v==1 && (j-s)%2==0) // p, ln on the same face
                  ev*=-1;
          break;
          case 2:
              ev= (this->m_polynomial[sz[0]*sz[1]-1]>0 ? 1:-1);
              if (v==2 && (j-s-e)%2==0) // p, ln on the same face
                  ev*=-1;
          break;
          case 3:
              ev= (this->m_polynomial[(sz[1]-1)*st[1]] >0 ? 1:-1);
              if (v==3 && (j-a+w)%2==0) // p, ln on the same face
                  ev*=-1;
          break;
          }

      // std::cout<<"Cell"<<this <<"("<<(sgn==1 ? "+": "-")<<")" <<std::endl;
      // std::cout<<"u= "<< u<<", v=" <<v<<", ev=" <<ev << std::endl; 
      // std::cout<<"ln= "<< ln->x()<<","<<ln->y() << std::endl; 
      // std::cout<<"p = "<< p->x()<<","<<p->y() << std::endl; 
      // std::cout<<"rn= "<< rn->x()<<","<<rn->y() << std::endl; 
      // if    (ev*sgn>0)  std::cout<<"result is ln"<< std::endl;
      // else              std::cout<<"result is rn"<< std::endl;



      if    (ev*sgn>0)  return ln;
      else              return rn;
      }
}

//====================================================================
} ; // namespace shape
} ; // namespace mmx
# undef Solver
# undef SELF
# undef REF
# undef TMPL
# undef STACK
# endif // shape_cell2d_algebraic_curve_hpp