include/shape/bcell2d_algebraic_curve.hpp

Go to the documentation of this file.

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

# include <realroot/Interval.hpp>
# include <realroot/ring_bernstein_tensor.hpp>
# include <realroot/Seq.hpp>
# include <shape/point.hpp>
# include <shape/solver_implicit.hpp>
# include <shape/bcell2d.hpp>
# include <shape/algebraic_curve.hpp>

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

# undef Cell_t
//====================================================================
namespace mmx { namespace shape {
//====================================================================
TMPL struct bcell2d_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 
bcell2d_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);

        tensor::split(left->m_polynomial, right->m_polynomial, v);
        if (left->e_intersections.size() == 0)
        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 ;
            right->n_intersections<< p ;
          } else 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 ;
            right->s_intersections<< p ;
          } else 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);


        tensor::split(left->m_polynomial, right->m_polynomial, v);
        if (left->n_intersections.size() ==0)
          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 ;
            right->w_intersections << p ;
          } else 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 ;
            right->e_intersections << p ;
          } else 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( );

}
//====================================================================
TMPL class mesher2d;
//====================================================================
template<class C, class V=default_env>
class bcell2d_algebraic_curve : public bcell2d<C,REF> {
public:
    typedef C                                          Scalar;
    typedef topology<C,REF>                            Topology;
    typedef typename Topology::Point                   Point;
    typedef Seq<Point*>                                Points;
    //typedef typename bcell2d_algebraic_curve_def<K>::Vertex         Vertex;
    //typedef typename use<bcell2d_algebraic_curve_def,V>::Vector      Vector;
    typedef typename use<point_def,C,V>::Point         Vector;
    typedef typename topology<C,REF>::Edge             Edge;
    typedef bcell<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 mesher2d<C,V>                              Mesher2d;

    bcell2d_algebraic_curve(const Polynomial &, const BoundingBox&, bool inter=true);
    bcell2d_algebraic_curve(const char*, const BoundingBox&);
    bcell2d_algebraic_curve(const typename AlgebraicSet::Polynomial& , const BoundingBox&);
    bcell2d_algebraic_curve(AlgebraicCurve*, const BoundingBox &x);
    bcell2d_algebraic_curve(const SELF& cl)
        : bcell2d<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;
    Seq<Point *>      m_boundary;
    //  AlgebraicCurve*   m_shape;
};

//====================================================================
TMPL inline void 
print(SELF* cv) {
    typedef typename SELF::Point Point;
    std::cout <<"\nBox       :  ["<<cv->xmin()<<","<<cv->xmax()<<"] *"
              <<" ["<<cv->ymin()<<","<<cv->ymax()<<"]";
    std::cout <<"\npoint(s) s: ";
    foreach(Point* p, cv->s_intersections)
        std::cout <<" ["<<p->x()<<" "<<p->y()<<"]";
    std::cout << "\npoint(s) e: ";
    foreach(Point* p, cv->e_intersections)
        std::cout <<" ["<<p->x()<<" "<<p->y()<<"]";
    std::cout << "\npoint(s) n: ";
    foreach(Point* p, cv->n_intersections)
        std::cout <<" ["<<p->x()<<" "<<p->y()<<"]";
    std::cout << "\npoint(s) w: ";
    foreach(Point* p, cv->w_intersections)
        std::cout <<" ["<<p->x()<<" "<<p->y()<<"]";

    std::cout<<std::endl;

    std::cout << "Point(s) sing: ";
    foreach(Point* p, cv->m_singular)
        std::cout <<" ["<<p->x()<<" "<<p->y()<<"]";
    std::cout<<std::endl;
}
//--------------------------------------------------------------------
TMPL
SELF::bcell2d_algebraic_curve(const Polynomial& pol, const BoundingBox& bx, bool inter): bcell2d<C,REF>(bx), m_polynomial(pol)
{
    //std::cout<<"Alegbraic bcell: "<<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::bcell2d_algebraic_curve(const typename AlgebraicSet::Polynomial& p, const BoundingBox & b): bcell2d<C,REF>(b)
{
    //std::cout<<"Alegbraic bcell: "<<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::bcell2d_algebraic_curve(const char* s, const BoundingBox & b): bcell2d<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->s_intersections, m_polynomial, Solver::south_edge, b);
    Solver::edge_point(this->e_intersections, m_polynomial, Solver::east_edge , b);
    Solver::edge_point(this->n_intersections, m_polynomial, Solver::north_edge, b);
    Solver::edge_point(this->w_intersections, m_polynomial, Solver::west_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::bcell2d_algebraic_curve(AlgebraicCurve* cv, const BoundingBox & b)
{
    new (this) SELF(cv->equation(),b );
}

//--------------------------------------------------------------------
TMPL bool 
SELF::is_regular() {
    // Assumption: bcells 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) {
    //Mesher2d* s= dynamic_cast<Mesher2d*>(t);

    int  sgn(1);
    Point *p,*q;

    //std::cout<<"alg. curve, inserting regular "<< *this <<std::endl;

    Scalar eps = 1e-10;//Solver::precision;


    foreach(Point* pt, this->s_intersections) {
      m_boundary<<pt;
    }

    if(this->s_intersections.size()>0) q=this->s_intersections.rep().back(); else q=NULL;
    foreach(Point* pt, this->e_intersections) {
      if(q && pt->dist2(*q) < eps) { 
        q=pt; 
      } else {
        m_boundary<<pt;
      }
    }

    if(this->e_intersections.size()>0) q=this->e_intersections.rep().back(); else q=NULL;
    foreach(Point* pt, this->n_intersections.reversed()) {
      if(q && pt->dist2(*q) < eps) { 
        q=pt; 
      } else {
        m_boundary<<pt;
      }
    }

    if(this->n_intersections.size()>0) q=this->n_intersections.rep().front(); else q=NULL;
    foreach(Point* pt, this->w_intersections.reversed()) {
      if(q && pt->dist2(*q) < eps) {
        q=pt; 
      } else {
        m_boundary<<pt;
      }
    }

    //l= this->intersections();
    //ni=l.size();
    Seq<Point*> l=m_boundary;

//    if( l.size() !=2) {
//        foreach( Point* e, l)
//            std::cout<<"("<< e->x()<<", "<<e->y()<<") ";
//        std::cout<<std::endl;
//    }

    if (l.size()==1 || l.size()==3) {
      //t->insert((BoundingBox*)this,false);
    }
    if(this->m_singular.size()==0)
    {
        foreach(Point* v, l) t->insert(v);

        if (l.size() ==2)
            t->insert_edge(l[0],l[1]);
        else {
            while (l.size()>1)
            {
                //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;
                t->insert_edge(p,q);
                //std::cout<<"edge "<<0<<" "<<q<<std::endl;
                //std::cout<<"l.size()="<<l.size() <<std::endl;
            }
            if (l.size()>0) {
                //t->insert((BoundingBox*)this,false);
                print(this);
                std::cout<< "Warning: point "<<*(l[0])<<" not used"<< 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;
        t->insert(c);
        foreach(Point* p, this->n_intersections) t->insert_edge(p,c);
        foreach(Point* p, this->e_intersections) t->insert_edge(p,c);
        foreach(Point* p, this->s_intersections) t->insert_edge(p,c);
        foreach(Point* p, this->w_intersections) t->insert_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;
    
    //Mesher2d* s = dynamic_cast<Mesher2d*>(t);

    int sgn(1);

    Seq<Point*> l;
    foreach(Point* p, this->n_intersections) {
        t->insert(p); l<<p;
        //qDebug()<<"Point idx:"<<p->index();
    }
    foreach(Point* p, this->e_intersections) {
        t->insert(p); l<<p;
        //qDebug()<<"Point idx:"<<p->index();
    }
    foreach(Point* p, this->s_intersections) {
        t->insert(p); l<<p;
        //qDebug()<<"Point idx:"<<p->index();
    }
    foreach(Point* p, this->w_intersections) {
        t->insert(p); l<<p;
        //qDebug()<<"Point idx:"<<p->index();
    }
    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;
        //qDebug()<<"Point"<<c->index();
        t->insert(c);
        foreach(Point* p, this->n_intersections) t->insert_edge(p,c);
        foreach(Point* p, this->e_intersections) t->insert_edge(p,c);
        foreach(Point* p, this->s_intersections) t->insert_edge(p,c);
        foreach(Point* p, this->w_intersections) t->insert_edge(p,c);
        return true;
    }
    else   // no singular points in the bcell
    {
        if(ni==2) { // 2 points on the boundary
            t->insert_edge(l[0],l[1]);
            return true;
        }
        else // more than one branch in the bcell
        {
            Point *p,*q;
            while (l.size()>0)
            {
                p= l[0];
                q= this->pair(p,sgn);
                l.erase(0);
                l.erase( l.search(q) );
                t->insert_edge(p,q);
            }
        }
    }
    
    //std::cout<<"singular bcell with "<<ni<<" bpts "<<this->m_singular.size()<<" epts"<<std::endl;

    return true;
}

//--------------------------------------------------------------------
TMPL void
SELF::subdivide(Cell*& Left, Cell*& Right, int v, double s) {
    bcell2d_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;

    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 bcell 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;
        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->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 (1)"<< 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 (2)"<< this<<std::endl; 

        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

---