synaps/msolve/BSolver.h

Go to the documentation of this file.

// Copyright (c) 2003  INRIA Sophia-Antipolis (France) and
//                     Department of Informatics and Telecommunications
//                     University of Athens (Greece).
// All rights reserved.
//
// Authors : Elias P. TSIGARIDAS <et@di.uoa.gr>

// Partially supported by INRIA's project "CALAMATA", a
// bilateral collaboration with National Kapodistrian University of
// Athens.
// Partially supported by the IST Programme of the EU as a Shared-cost
// RTD (FET Open) Project under Contract No  IST-2000-26473
// (ECG - Effective Computational Geometry for Curves and Surfaces)
//
//--------------------------------------------------------------------
#ifndef SYNAPS_ARITHM_ALGEBRAIC_BSOLVER_H
#define SYNAPS_ARITHM_ALGEBRAIC_BSOLVER_H
//--------------------------------------------------------------------
#include <synaps/init.h>
#include <synaps/usolve/Algebraic.h>
#include <synaps/usolve/algebraic/sign_fct_2.h>
//#include <synaps/usolve/root_of.h>
#include <synaps/mpol/MPoly_2_2.h>

#include <synaps/usolve/sturm/solve_small.h>
//#include <synaps/solve/sturmsmdeg.h>
#include <synaps/linalg/VectDse.h>

//--------------------------------------------------------------------
__ALGEBRAIC_BEGIN_NAMESPACE
//--------------------------------------------------------------------
//--------------------------------------------------------------------


template<class A>
VectDse<A> make_vector(const A& a, const A& b)
{
  VectDse<A> res(2,AsSize());
  res[0]=a;res[1]=b;
  return res;
}



template < typename SOLVER,
           typename OutputIterator >
OutputIterator
solve_for_circles( const typename SOLVER::MPoly_2_2& f1, 
                   const typename SOLVER::MPoly_2_2& f2, 
                   OutputIterator res, 
                   const SOLVER& mth)
{
    typedef typename SOLVER::Poly     POLY;
    typedef typename SOLVER::RO_t     RO_t;
    typedef typename SOLVER::FT       FT;

    POLY Rx = resultant_y<POLY>(f1, f2);
    POLY Ry = resultant_x<POLY>(f1, f2);
    //  cout << "Rx = " << Rx << endl;
    //  cout << "Ry = " << Ry << endl;

    Seq<RO_t> solX;
    Seq<RO_t> solY;

    solve(Rx, std::back_inserter(solX.rep()), mth);
    solve(Ry, std::back_inserter(solY.rep()), mth);
 
    // The circles have no intersections
    if ((solX.size() == 0) || (solY.size() == 0)) return res;

    // The circles have a tangent intersection
    if ((solX.size() == 1) && (solY.size() == 1)){
        *res++ = std::make_pair( make_vector(solX[0], solY[0]), static_cast<unsigned>(2) );
        return res;
    }

    // Two intersections but with the same X-coordinate
    if ((solX.size() == 1) && (solY.size() == 2)) {
        *res++ = std::make_pair( make_vector(solX[0], solY[0]), static_cast<unsigned>(1) );
        *res++ = std::make_pair( make_vector(solX[0], solY[1]), static_cast<unsigned>(1) );
        return res;
    }

    // Two intersections but with the same Y-coordinate
    if ((solX.size() == 2) && (solY.size() == 1)) {
        *res++ = std::make_pair( make_vector(solX[0], solY[0]), static_cast<unsigned>(1) );
        *res++ = std::make_pair( make_vector(solX[1], solY[0]), static_cast<unsigned>(1) );
        return res;
    }

    // Two distinct intersections
    FT dx = f2.xc() - f1.xc();
    FT dy = f2.yc() - f1.yc();
  
    // From the slope decide how to match the coordinates
    bool slope = (sign(dx) * sign(dy) <= 0);
    if (slope) {
        *res++ = std::make_pair( make_vector(solX[0], solY[0]),static_cast<unsigned>(1) );
        *res++ = std::make_pair( make_vector(solX[1], solY[1]),static_cast<unsigned>(1) );
    } else {
        *res++ = std::make_pair( make_vector(solX[0], solY[1]),static_cast<unsigned>(1) );
        *res++ = std::make_pair( make_vector(solX[1], solY[0]),static_cast<unsigned>(1) );
    }

    return res;
}



template < typename SOLVER,
           typename OutputIterator >
OutputIterator
fast_solve(const typename SOLVER::MPoly_2_2& f1, 
           const typename SOLVER::MPoly_2_2& f2, 
           OutputIterator ii,
           const SOLVER& mth)
{
    typedef typename SOLVER::Poly     POLY;
    typedef typename SOLVER::RO_t     RO_t;
    typedef typename SOLVER::FT       FT;

    POLY Rx = resultant_y<POLY>(f1, f2);
    POLY Ry = resultant_x<POLY>(f1, f2);
    //  cout << "Rx = " << Rx << endl;
    //  cout << "Ry = " << Ry << endl;

    Seq<RO_t> solx = Solve(Rx, mth);
    Seq<RO_t> soly = Solve(Ry, mth);
  
    // Indentify the multiplicities of the roots
    std::vector<int> x_multiplicity(solx.size());
    for (unsigned int i = 0; i < solx.size(); ++i) {
        x_multiplicity[i] = solx[i].multiplicity();
    }

    std::vector<int> y_multiplicity(soly.size());
    for (unsigned int i = 0; i < soly.size(); ++i) {
        y_multiplicity[i] = soly[i].multiplicity();
    }

    //  cout << "--- Testing with StHa -------------" << endl;
    typename Seq<RO_t>::const_iterator tx, ty;
    int i = 0;
    for (tx = solx.begin(); tx != solx.end(); ++tx, ++i) {

        POLY xPoly(tx->poly());
        FT x_lisol = tx->left();
        std::vector<POLY> f1_leftStHa =
            compute_StHa_x_cx2y2(f1, xPoly, numerator(x_lisol), denominator(x_lisol));

        FT x_risol = tx->right();
        std::vector<POLY> f1_rightStHa =
            compute_StHa_x_cx2y2(f1, xPoly, numerator(x_risol), denominator(x_risol));

        std::vector<POLY> f2_leftStHa;
        std::vector<POLY> f2_rightStHa;
        bool second_computed = false;

        int j = 0;
        for (ty = soly.begin(); ty != soly.end(); ++ty, ++j) {
            //    cout << "m =" <<  x_multiplicity[i] << " " << y_multiplicity[j] << endl;
            if ((x_multiplicity[i] == 0) || (y_multiplicity[j] == 0)) continue;
            //    cout << "after " << endl;

            int sgn_a = opposite((tx->degree()-1-tx->which()) % 2 ? NEGATIVE : POSITIVE);
            int sgn_b = opposite(sgn_a);

            int V1;
            if (tx->degree() == 1) {
                V1 = sign_at(f1_leftStHa[1], *ty);
            } else {
                int Va = sign_on_StHa(f1_leftStHa,  sgn_a, *ty);
                int Vb = sign_on_StHa(f1_rightStHa, sgn_b, *ty);
                V1 = Va - Vb;
            }
            //    cout << "V1 = " << V1 << endl;
            if (V1) continue;

            if (!second_computed) {
                //              cout << " In here " << endl;
                f2_leftStHa =
                    compute_StHa_x_cx2y2(f2, xPoly, numerator(x_lisol), denominator(x_lisol));

                f2_rightStHa =
                    compute_StHa_x_cx2y2(f2, xPoly, numerator(x_risol), denominator(x_risol));
                second_computed = true;
            }
            int V2;
            if (tx->degree() == 1) {
                V2 = sign_at(f2_leftStHa[1], *ty);
            } else {
                int Va = sign_on_StHa(f2_leftStHa,  sgn_a, *ty);
                int Vb = sign_on_StHa(f2_rightStHa, sgn_b, *ty);
                V2 = Va - Vb;
            }
            //    cout << "V2 = " << V2 << endl;
            if (V2) continue;

            // We found a root
            *ii++ = make_vector(*tx, *ty);
            --x_multiplicity[i];
            --y_multiplicity[j];
        }
        //      cout << " ---- end x ----- " <<endl;
    }
    return ii;
}



// Compute the solutions of a system of 2 conics, using only bivariate
// @sign_at@ functions.
// Note however that the total number of intersections cannot exceeds the multilicity
// of @X@-coordiante (resp. @Y@). This is taken into account.
// It is very difficult with this functions to compute the multiplicities of the
// intersection points.
template < typename SOLVER,
           typename OutputIterator >
OutputIterator
trivial_solve(const typename SOLVER::MPoly_2_2& f1, 
              const typename SOLVER::MPoly_2_2& f2,
              OutputIterator res,
              const SOLVER& mth)
{
    std::cout << __FUNCTION__ << std::endl;
  
    typedef typename SOLVER::Poly     POLY;
    typedef typename SOLVER::RO_t     RO_t;

    POLY Rx = resultant_y<POLY>(f1, f2);
    POLY Ry = resultant_x<POLY>(f1, f2);

    Seq<RO_t> Xsol = Solve(Rx, mth);
    Seq<RO_t> Ysol = Solve(Ry, mth);

    Seq<int> Xmult(Xsol.size());
    for (int i = 0; i < Xmult.size(); ++i) Xmult[i] = Xsol[i].multiplicity();

    Seq<int> Ymult(Ysol.size());
    for (int i = 0; i < Ymult.size(); ++i) Ymult[i] = Ysol[i].multiplicity();

    for (int i = 0; i < Xsol.size(); ++i) {
        for (int j = 0; j < Ysol.size(); ++j) {
            if (Ymult[j] == 0) continue;
            if (Xmult[i] == 0) break;
            if ( (sign_at(f1, Xsol[i], Ysol[j]) == ZERO) && 
                 (sign_at(f2, Xsol[i], Ysol[j]) == ZERO) ) {
                *res++ = make_vector( Xsol[i], Ysol[j] );
                --Xmult[i];
                --Ymult[j];
            }
        }
    }
    return res;
}




// Check if the 1st subresultant evaluated over x
// is inside the isolating interval of y
template < typename RO_t >
int
check_where(const RO_t& y, const RO_t& x,
            const typename RO_t::POLY& A,
            const typename RO_t::POLY& B)
{
    typedef typename RO_t::RT     RT;
    typedef typename RO_t::FT     FT;
    typedef typename RO_t::POLY   POLY;

    POLY P;
    FT q;


#if 0
    // TODO: Do a more clever approximation!!!

    double xx = to_double(x);
    double  apx = UPOLDAR::eval_horner(-B, xx) /
        UPOLDAR::eval_horner(A, xx);

    if (y.interval().lower() > apx) return SMALLER;
    if (y.interval().upper() < apx) return LARGER;
    return EQUAL;
#endif

    q = y.left();
    P = -B * denominator(q) - A * numerator(q);
    if (sign_at(P, x) * sign_at(A, x) < 0) return SMALLER;

    q = y.right();
    P = -B * denominator(q) - A * numerator(q);
    if (sign_at(P, x) * sign_at(A, x) > 0) return LARGER;

    return EQUAL;

}



template < typename SOLVER,
           typename OutputIterator >
OutputIterator
solve(const typename SOLVER::MPoly_2_2& f1, 
      const typename SOLVER::MPoly_2_2& f2,
      OutputIterator ii, 
      const SOLVER& mth)
{
    typedef typename SOLVER::RT       RT;
    typedef typename SOLVER::FT       FT;
  
    typedef typename SOLVER::Poly     POLY;
    typedef typename SOLVER::RO_t     RO_t;

    POLY Ax, Bx;
    POLY Rx;
    resultant_y(f1, f2, Ax, Bx, Rx);

    POLY Ay, By;
    POLY Ry;
    resultant_x(f1, f2, Ay, By, Ry);


    std::vector<RO_t> solx;
    solve_small_degree(Rx, 1, std::back_inserter(solx), mth);

    std::vector<RO_t> soly;
    solve_small_degree(Ry, 1, std::back_inserter(soly), mth);

  
    bool x_has_multiplicities = false;
    bool y_has_multiplicities = false;

    // Indentify the multiplicities of the roots
    std::vector<int> x_multiplicity(solx.size());
    for (unsigned int i = 0; i < solx.size(); ++i) {
        x_multiplicity[i] = solx[i].multiplicity();
        if (x_multiplicity[i] > 1) x_has_multiplicities = true;
    }

    std::vector<int> y_multiplicity(soly.size());
    for (unsigned int i = 0; i < soly.size(); ++i) {
        y_multiplicity[i] = soly[i].multiplicity();
        if (y_multiplicity[i] > 1) y_has_multiplicities = true;
    }

    // if ( (solx.size() == 2) && (soly.size() == 2))
//     {
//      std::cout << "f1: " << f1 << std::endl;
//      std::cout << "f2: " << f2 << std::endl;

//      for (int i = 0; i < solx.size(); ++i)
//          std::cout << solx[i].multiplicity() << " ";
//      std::cout << std::endl;
//      for (int i = 0; i < soly.size(); ++i)
//          std::cout << soly[i].multiplicity() << " ";
//      std::cout << std::endl;
//     }

    bool revert = false;
    bool solve_by_subres = false;

    if ( (!x_has_multiplicities) && (!y_has_multiplicities) ) {
        solve_by_subres = true;
    }

    POLY A, B, R;
    if (!revert) {
        A = Ax; B = Bx; R = Rx;
    } else {
        std::swap(solx, soly);
        A = Ay; B = By; R = Ry;
    }

    // We are in generic position, use the 1st subresultant
    if (solve_by_subres) {
        for (typename std::vector<RO_t>::iterator xt = solx.begin(); xt != solx.end(); ++xt) {
            int start = 0;
            int end = soly.size() - 1;

            while (start <= end) {
                int mid = static_cast<int>(std::ceil((end-start)/2.0));
                int pos = start + mid;

                int result = check_where(soly[pos], *xt, A, B);
                if (result == EQUAL) {
                    if (!revert) { *ii++ = make_vector(*xt, soly[pos]); }
                    else {      *ii++ = make_vector(soly[pos], *xt); }
                    typename std::vector<RO_t>::iterator jt = soly.begin() + pos;
                    soly.erase(jt);
                    break;
                }
                if (result == SMALLER) {
                    end = pos - 1;
                    continue;
                }
                // result == GREATER
                start = pos + 1;
            }
        }
        return ii;
    }

    // We are not in generic position
    // cout << "--- Testing with StHa -------------" << endl;
    typename Seq<RO_t>::const_iterator tx, ty;
    int i = 0;
    for (tx = solx.begin(); tx != solx.end(); ++tx, ++i) {

        POLY xPoly(tx->poly());
        FT x_lisol = tx->left();
        std::vector<POLY> f1_leftStHa =
            compute_StHa_x_cx2y2(f1, xPoly, numerator(x_lisol), denominator(x_lisol));

        FT x_risol = tx->right();
        std::vector<POLY> f1_rightStHa =
            compute_StHa_x_cx2y2(f1, xPoly, numerator(x_risol), denominator(x_risol));

        std::vector<POLY> f2_leftStHa;
        std::vector<POLY> f2_rightStHa;
        bool second_computed = false;

        int j = 0;
        for (ty = soly.begin(); ty != soly.end(); ++ty, ++j) {
            if ((x_multiplicity[i] == 0) || (y_multiplicity[j] == 0)) continue;

            int sgn_a = opposite((tx->degree()-1-tx->which()) % 2 ? NEGATIVE : POSITIVE);
            int sgn_b = opposite(sgn_a);
            int V1;
            if (tx->degree() == 1) {
                V1 = sign_at(f1_leftStHa[1], *ty);
            } else {
                int Va = sign_on_StHa(f1_leftStHa,  sgn_a, *ty);
                int Vb = sign_on_StHa(f1_rightStHa, sgn_b, *ty);
                V1 = Va - Vb;
            }
            if (V1) continue;

            if (!second_computed) {
                f2_leftStHa =
                    compute_StHa_x_cx2y2(f2, xPoly, numerator(x_lisol), denominator(x_lisol));

                f2_rightStHa =
                    compute_StHa_x_cx2y2(f2, xPoly, numerator(x_risol), denominator(x_risol));
                second_computed = true;
            }
            int V2;
            if (tx->degree() == 1) {
                V2 = sign_at(f2_leftStHa[1], *ty);
            } else {
                int Va = sign_on_StHa(f2_leftStHa,  sgn_a, *ty);
                int Vb = sign_on_StHa(f2_rightStHa, sgn_b, *ty);
                V2 = Va - Vb;
            }
            if (V2) continue;

            // We have found a root
            *ii++ = make_vector(*tx, *ty);
            --x_multiplicity[i];
            --y_multiplicity[j];
        }
        //      cout << " ---- end x ----- " <<endl;
    }


    return ii;
}



template < typename SOLVER,
           typename OutputIterator >
OutputIterator
solve_in_range(const typename SOLVER::MPoly_2_2& f1, 
               const typename SOLVER::MPoly_2_2& f2, 
               OutputIterator ii,
               const typename SOLVER::RO_t& start,
               const typename SOLVER::RO_t& end,
               const SOLVER& mth)
{
    typedef typename SOLVER::RT       RT;
    typedef typename SOLVER::FT       FT;

    typedef typename SOLVER::Poly     POLY;
    typedef typename SOLVER::RO_t     RO_t;

    POLY Ax, Bx;
    POLY Rx;
    resultant_y(f1, f2, Ax, Bx, Rx);

    POLY Ay, By;
    POLY Ry;
    resultant_x(f1, f2, Ay, By, Ry);

    Seq<RO_t> solx = Solve(Rx, mth);
    Seq<RO_t> soly = Solve(Ry, mth);
  
    bool x_has_multiplicities = false;
    bool y_has_multiplicities = false;

    // Indentify the multiplicities of the roots
    std::vector<int> x_multiplicity(solx.size());
    for (unsigned int i = 0; i < solx.size(); ++i) {
        x_multiplicity[i] = solx[i].multiplicity();
        if (x_multiplicity[i] > 1) x_has_multiplicities = true;
    }

    std::vector<int> y_multiplicity(soly.size());
    for (unsigned int i = 0; i < soly.size(); ++i) {
        y_multiplicity[i] = soly[i].multiplicity();
        if (y_multiplicity[i] > 1) y_has_multiplicities = true;
    }

    bool solve_by_subres = false;

    if ( (!x_has_multiplicities) && (!y_has_multiplicities) ) {
        solve_by_subres = true;
    }

    POLY A(Ax);
    POLY B(Bx);
    POLY R(Rx);

    // We are in generic position, use the 1st subresultant
    if (solve_by_subres) {
        for (typename Seq<RO_t>::const_iterator xt = solx.begin(); xt != solx.end(); ++xt) {
            if ((*xt >= start) && (*xt <= end)) {
                int start = 0;
                int end = soly.size() - 1;

                while (start <= end) {
                    int mid = static_cast<int>(std::ceil((end-start)/2.0));
                    int pos = start + mid;

                    int result = check_where(soly[pos], *xt, A, B);
                    if (result == EQUAL) {
                        *ii++ = make_vector(*xt, soly[pos]);
                        break;
                    }
                    if (result == SMALLER) {
                        end = pos - 1;
                        continue;
                    }
                    // result == GREATER
                    start = pos + 1;
                }
            }
        }
        return ii;
    }

    // We are not in generic position
    //  cout << "--- Testing with StHa -------------" << endl;
    typename Seq<RO_t>::const_iterator tx, ty;
    int i = 0;
    for (tx = solx.begin(); tx != solx.end(); ++tx, ++i) {
        if ((*tx >= start) && (*tx <= end)) {

            POLY xPoly(tx->poly());
            FT x_lisol = tx->left();
            std::vector<POLY> f1_leftStHa =
                compute_StHa_x_cx2y2(f1, xPoly, numerator(x_lisol), denominator(x_lisol));

            FT x_risol = tx->right();
            std::vector<POLY> f1_rightStHa =
                compute_StHa_x_cx2y2(f1, xPoly, numerator(x_risol), denominator(x_risol));

            std::vector<POLY> f2_leftStHa;
            std::vector<POLY> f2_rightStHa;
            bool second_computed = false;

            int j = 0;
            for (ty = soly.begin(); ty != soly.end(); ++ty, ++j) {
                if ((x_multiplicity[i] == 0) || (y_multiplicity[j] == 0)) continue;

                int sgn_a = opposite((tx->degree()-1-tx->which()) % 2 ? NEGATIVE : POSITIVE);
                int sgn_b = opposite(sgn_a);
                int V1;
                if (tx->degree() == 1) {
                    V1 = sign_at(f1_leftStHa[1], *ty);
                } else {
                    int Va = sign_on_StHa(f1_leftStHa,  sgn_a, *ty);
                    int Vb = sign_on_StHa(f1_rightStHa, sgn_b, *ty);
                    V1 = Va - Vb;
                }
                if (V1) continue;

                if (!second_computed) {
                    f2_leftStHa =
                        compute_StHa_x_cx2y2(f2, xPoly, numerator(x_lisol), denominator(x_lisol));

                    f2_rightStHa =
                        compute_StHa_x_cx2y2(f2, xPoly, numerator(x_risol), denominator(x_risol));
                    second_computed = true;
                }
                int V2;
                if (tx->degree() == 1) {
                    V2 = sign_at(f2_leftStHa[1], *ty);
                } else {
                    int Va = sign_on_StHa(f2_leftStHa,  sgn_a, *ty);
                    int Vb = sign_on_StHa(f2_rightStHa, sgn_b, *ty);
                    V2 = Va - Vb;
                }
                if (V2) continue;

                // We have found a root
                *ii++ = make_vector(*tx, *ty);
                --x_multiplicity[i];
                --y_multiplicity[j];
            }
            //  cout << " ---- end x ----- " <<endl;
        }
    }
    return ii;
}



template < typename SOLVER, typename O, typename R >
Seq< VectDse< typename SOLVER::RO_t > >
Solve_with_sign( const MPol<typename SOLVER::RT, O, R>& s1, 
                 const MPol<typename SOLVER::RT, O, R>& s2,
                 const SOLVER& mth)
{
    //std::cout << __FUNCTION__ << std::endl;
    assert((NbVar(s1) <= 2) && (NbVar(s2) <= 2));

    typedef typename SOLVER::Poly       Poly;
    typedef typename SOLVER::Poly_2     Poly_2;
    typedef typename SOLVER::RO_t       RO_t;

    Poly_2 f1x = let::convert_to_upoly_2(s1, POLY_wrt_Y);
    Poly_2 f2x = let::convert_to_upoly_2(s2, POLY_wrt_Y);

    SturmSeq<Poly> stha(f1x, f2x, HABICHT());
    Poly Rx = stha[stha.size() - 1][0];
    Seq<RO_t> solx = Solve(Rx, SOLVER());


    Poly_2 f1y = let::convert_to_upoly_2(s1, POLY_wrt_X);
    Poly_2 f2y = let::convert_to_upoly_2(s2, POLY_wrt_X);

    SturmSeq<Poly> sthb(f1y, f2y, HABICHT());
    Poly Ry = sthb[sthb.size() - 1][0];
    Seq< RO_t > soly = Solve(Ry, SOLVER());

#if 0
    std::cout << "Rx: " << Rx << std::endl;
    std::cout << "Ry: " << Ry << std::endl;
    for (int i = 0; i < solx.size(); ++i) {
        std::cout << to_double(solx[i]) << std::endl;
    }
    std::cout << "-------" << std::endl;
    for (int i = 0; i < soly.size(); ++i) {
        std::cout << to_double(soly[i]) << std::endl;
    }
#endif

    typedef typename Seq<RO_t>::const_iterator  IT;
    std::vector< VectDse<RO_t> > sol;
    for (IT xt = solx.begin(); xt != solx.end(); ++xt) {
        for (IT yt = soly.begin(); yt != soly.end(); ++yt) {
            if ((sign_at(f1x, *xt, *yt) == ZERO) && (sign_at(f2x, *xt, *yt) == ZERO)) {
                sol.push_back( make_vector(*xt, *yt) );
            }
        }
    }

    return sol;
}




template < typename SOLVER, typename O, typename R >
Seq< std::pair< typename SOLVER::RO_t, typename SOLVER::RO_t > >
Solve_with_subres( const MPol<typename SOLVER::RT, O, R>& F1, 
                   const MPol<typename SOLVER::RT, O, R>& F2,
                   const SOLVER& mth)
{  
    assert((NbVar(F1) <= 2) && (NbVar(F2) <= 2));
 
    typedef typename SOLVER::RT    RT;

    typedef MPol<RT, O, R>              MPoly;

    typedef typename SOLVER::Poly        Poly;
    typedef typename SOLVER::Poly_2      Poly_2;
    typedef typename SOLVER::RO_t        RO_t;

    Poly_2 f1y = let::convert_to_upoly_2(F1, POLY_wrt_Y);
    Poly_2 f2y = let::convert_to_upoly_2(F2, POLY_wrt_Y);

    Poly_2 f1x = let::convert_to_upoly_2(F1, POLY_wrt_X);
    Poly_2 f2x = let::convert_to_upoly_2(F2, POLY_wrt_X);

    int dg_1 = UPOLDAR::degree( UPOLDAR::lcoeff(f1y) );
    int dg_2 = UPOLDAR::degree( UPOLDAR::lcoeff(f2y) );
    if ((dg_1 >= 1) && (dg_2 >= 1))
    {
        std::cout << "Linear transformation" << std::endl;
        Poly_2 h = construct_X_minus_aY(RT(1));
        Poly_2 tf1y = compose(f1y, h);
        Poly_2 tf2y = compose(f2y, h);

        //Poly_2 tf1y = changeOfVariable(f1y);
        //Poly_2 tf2y = changeOfVariable(f2y);

        //tf1y = changeOfVariable(tf1y);
        //tf2y = changeOfVariable(tf2y);

        Poly_2 tf1x = interchange_variables(tf1y);
        Poly_2 tf2x = interchange_variables(tf2y);

        std::swap(f1y, tf1y);
        std::swap(f2y, tf2y);
        std::swap(f1x, tf1x);
        std::swap(f2x, tf2x);
    }
    if ((dg_1 >= 1) && (dg_2 < 1)) {
        std::cout << "Do a swap" << std::endl;
        std::swap(f1y, f2y);
        std::swap(f1x, f2x);
    }

    //std::cout << "Done init" << std::endl;
    SturmSeq<Poly> StHa_Y(f1y, f2y, HABICHT());
    Poly Rx = StHa_Y[StHa_Y.size() - 1][0];
    //Seq< root_of<RT> > solx = Solve(Rx, SOLVER<RT>());


    SturmSeq<Poly> StHa_X(f1x, f2x, HABICHT());
    Poly Ry = StHa_X[StHa_X.size() - 1][0];
    Seq<RO_t> solY = Solve(Ry, SOLVER());

    Poly Resultant = SquareFree(Rx);

    Seq<Poly> stha_list = CoeffStHa< Seq<Poly> >(StHa_Y);
    Seq<Poly> decomposition_list;
    Poly Fi = Resultant;
    for (int i = 0; i < stha_list.size() -1 ; ++i) {
        std::cout << Fi << ", " << stha_list[stha_list.size()-i-2]
                  << " gcd = " << PGcd(Fi, stha_list[stha_list.size()-i-2]) << std::endl;
        Poly Fin = PGcd(Fi, stha_list[stha_list.size()-i-2]);
        std::cout << " / = " << (Fi / Fin) << std::endl;
        decomposition_list.push_back(Fi / Fin);
        Fi = Fin;
    }

#if 1
    {
        std::cout << "Rx: " << Rx << std::endl;
        Seq< Poly > dec = SquareFreeFactorization(Rx);
        for (int i = 0; i < dec.size(); ++i) {
            std::cout << i << ": " << dec[i] << std::endl;
        }
        std::cout << "-----------------------------" << std::endl;
    }
    {
        std::cout << "Ry: " << Ry << std::endl;
        Seq< Poly > dec = SquareFreeFactorization(Ry);
        for (int i = 0; i < dec.size(); ++i) {
            std::cout << i << ": " << dec[i] << std::endl;
        }
        std::cout << "-----------------------------" << std::endl;
        std::cout << "Done resultants" << std::endl;
    }

    std::cout << "psc: " << std::endl;
    for (int i = 0; i < stha_list.size(); ++i) {
        std::cout << i << ": " << stha_list[i] << std::endl;
    }
    std::cout << std::endl;
    std::cout << "decomp" << std::endl;
    for (int i = 0; i < decomposition_list.size(); ++i) {
        std::cout << i << ": " << decomposition_list[i] << std::endl;
    }
    std::cout << std::endl;
#endif


    Seq< VectDse< RO_t > > solutions;

    //std::cout << "Decomp sz: " << decomposition_list.size() << std::endl;
    for (int j = 0; j < decomposition_list.size(); ++j) {
        if (UPOLDAR::degree(decomposition_list[j]) < 1) continue;

        //cout << "start ------------------------------------------\n";
        //cout << j << ": " << decomposition_list[j] << endl;

        Seq<RO_t> solX = Solve(decomposition_list[j], SOLVER());
        //std::cout << "solX sz: " << solX.size() << std::endl;

        Poly_2 tp = StHa_Y[ StHa_Y.size() - j - 2 ];
        if ((UPOLDAR::degree(tp) >= j) && (j < (StHa_Y.size() - 3))) {
            // RUR
            //cout << StHa_Y.size() << ", ";
            //cout << UPOLDAR::degree(tp) << ", "
            //   << j << ", "
            //   << (UPOLDAR::degree(tp) >= j) << endl;
            //std::cout << j << " tp = " << tp << ", dg: " << UPOLDAR::degree(tp) << endl;
            //cout << "StHa: " << StHa_Y[ StHa_Y.size() - j - 2] << endl;
            Poly B =  tp[j];
            Poly A =  tp[j+1] * (j+1);
            for (int k = 0; k < solX.size(); ++k) {
                int start = 0;
                int end = solY.size() - 1;
                while (start <= end) {
                    int mid = static_cast<int>(std::ceil((end-start)/2.0));
                    int pos = start + mid;

                    int result = check_where(solY[pos], solX[k], A, B);
                    if (result == EQUAL) {
                        solutions.push_back(make_vector(solX[k], solY[pos]));
                        std::cout << "solution\n";
                        std::cout << to_double(solX[k]) << ", "
                                  << to_double(solY[pos])
                                  << " mu ="<< j+1 << std::endl;
                        break;
                    }
                    if (result == SMALLER) {
                        end = pos - 1;
                        continue;
                    }
                    // result == GREATER
                    start = pos + 1;
                }
            }
        } else {
            std::cout << "Degenerate case ?" << std::endl;
        }
        //cout << "end------------------------------------------\n";
    }
    return solutions;
}





// compute the relative position of an algebraic number wrt two rational numbers.
// The algebraic number is in RUR (Rational Univariate Representation)
// i.e  y = -B(r) / A(r)
// This function is needed in a series of computations, and so the type of the
// rational number is a pair, where the first element is the number and the second
// is the order of y wrt to this number, if we know it already.
template < typename SOLVER >
int
pos_wrt_point( std::pair< typename SOLVER::FT, int >& isol_point,
               const typename SOLVER::RO_t& r,
               const typename SOLVER::Poly& A, 
               const typename SOLVER::Poly& B,
               const SOLVER& )
{
    typedef typename SOLVER::RT    RT;
    typedef typename SOLVER::FT    FT;
    typedef typename SOLVER::Poly  Poly;
    typedef typename SOLVER::RO_t  RO_t;

    if (isol_point.second != UNKNOWN) return isol_point.second;    

    int sgnA = sign_at(A, r);
  
    FT q = isol_point.first;
    Poly P  =  B * denominator(q) + A * numerator(q);

    if (sign_at(P, r) * sgnA > 0) isol_point.second = SMALLER;
    else isol_point.second = LARGER;

    return isol_point.second;
}

// IP: are the the intermediate points
template < typename SOLVER >
int
compute_y_in_generic_position( const Seq< typename SOLVER::FT >& IP,
                               const typename SOLVER::RO_t& r,
                               const typename SOLVER::Poly& A, 
                               const typename SOLVER::Poly& B,
                               const SOLVER& )
{
    typedef typename SOLVER::RT    RT;
    typedef typename SOLVER::FT    FT;
    typedef typename SOLVER::Poly  Poly;
    typedef typename SOLVER::RO_t  RO_t;

    for (int i = 0; i < IP.size(); ++i) {
        std::cout << to_double(IP[i]) << "  ";
    }
    std::cout << std::endl;
    double apx = - UPOLDAR::eval_horner(B, to_double(r)) 
        / UPOLDAR::eval_horner(A, to_double(r));
    std::cout << "(" << to_double(r) << ", " << apx << ")" << std::endl;
    
    


    Seq< std::pair<FT, int> > P(IP.size());
    for (int i = 0; i < IP.size(); ++i) {
        P[i] = make_vector( IP[i], UNKNOWN );
    }
    P[0].second = LARGER;
    P[P.size()-1].second = SMALLER;
    
    int start = 0;
    int end = P.size() - 1;
    while (start <= end) {
        if ((end - start) == 1) {
            assert( (P[start].second == LARGER) && (P[end].second == SMALLER));
            return start;
        }
        
        int pos  = static_cast<int>( std::ceil((start+end)/2.0) );
        int result = pos_wrt_point(P[pos], r, A, B, SOLVER());
        
        std::cout << "s, p, e : " << start << ", " << pos << ", " << end << std::endl; 
        std::cout << "In " << to_double(IP[start]) << ", " << to_double(IP[end]) << std::endl; 
        std::cout << "Testing point: " << to_double(P[pos].first) 
                  << " found " << result << std::endl;
        if (result == SMALLER) { 
            end = pos;
            continue;
        }
        if (result == LARGER) {
            start = pos;
            continue;
        }
    }
    std::cout << "We must always find a solution" << std::endl;
    return -1;
}




template < typename SOLVER,
           typename OutputIterator >
OutputIterator
New_solve(const typename SOLVER::MPoly_2_2& f1, 
          const typename SOLVER::MPoly_2_2& f2,
          OutputIterator res, 
          const SOLVER& mth)
{
    std::cout << __FUNCTION__ << std::endl;
    
    typedef typename SOLVER::RT       RT;
    typedef typename SOLVER::FT       FT;
  
    typedef typename SOLVER::Poly     Poly;
    typedef typename SOLVER::RO_t     RO_t;

    Poly Ax, Bx;
    Poly Rx;
    resultant_y(f1, f2, Ax, Bx, Rx);

    Poly Ay, By;
    Poly Ry;
    resultant_x(f1, f2, Ay, By, Ry);


    Seq<RO_t> Xsol = Solve(Rx, SOLVER());
    Seq<RO_t> Ysol = Solve(Ry, SOLVER());

    if ((Xsol.size() == 0) || (Ysol.size() == 0)) return res;
   

    std::vector<int> Xmult( Xsol.size() );
    for (int i = 0; i < Xmult.size(); ++i) Xmult[i] = Xsol[i].multiplicity();
    std::vector<int> possible_x_roots = Xmult;
  
  
    std::vector<int> Ymult( Ysol.size() );
    for (int i = 0; i < Ymult.size(); ++i) Ymult[i] = Ysol[i].multiplicity();
    std::vector<int> possible_y_roots = Ymult;

    unsigned sol_matrix[4][4] = { {0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0} };
            
    Seq<FT> Y_isol_points;
    compute_intermediate_points(Ysol, Y_isol_points);
  

    // std::copy( Xmult.begin(), Xmult.end(),
//             std::ostream_iterator<int>(std::cout, " " ));
//     std::cout << std::endl;
//     std::copy( Ymult.begin(), Ymult.end(),
//             std::ostream_iterator<int>(std::cout, " " ));
//     std::cout << std::endl;  
    Poly A(Ax);
    Poly B(Bx);
    Poly R(Rx);

    for (int i = 0; i < Xsol.size(); ++i) {  
        if (Xsol[i].multiplicity() == 1) {
            std::cout << "We are in generic position" << std::endl;
            int j = compute_y_in_generic_position(Y_isol_points, Xsol[i], A, B, SOLVER());
            sol_matrix[i][j] = 1;
            //std::cout << "found: " << i << ", " << j << std::endl; 
            --possible_x_roots[i];
            --possible_y_roots[j];
            continue;
        }
 
        if (Xsol[i].is_rational()) {
            std::cout << "Check with rational" << std::endl;
            // Xsol[i] is rational so left() == right() == Xsol[i]        
            FT r = Xsol[i].left();
          
            Poly sf1 = specialize_x<Poly>(f1, r);
            Poly sf2 = specialize_x<Poly>(f2, r);

            for (int j = 0; ((j < Ysol.size()) && possible_x_roots[i]); ++j) {
                if (possible_y_roots[j]) {
                    if ((sign_at(sf1, Ysol[j]) == ZERO) && (sign_at(sf2, Ysol[j]) == ZERO)) {
                        --possible_x_roots[i];
                        --possible_y_roots[j];
                        sol_matrix[i][j] = 1;
                    }
                }
            }
            continue;       
        }
 
        std::cout << "Check with sign_at" << std::endl;
        for (int j = 0; ((j < Ysol.size()) && possible_x_roots[i]); ++j) {
            if (possible_y_roots[j]) {
                if ( (sign_at(f1, Xsol[i], Ysol[j]) == ZERO) && 
                     (sign_at(f2, Xsol[i], Ysol[j]) == ZERO) ) {
                    --possible_x_roots[i];
                    --possible_y_roots[j];
                    sol_matrix[i][j] = 1;
                }
            }
        }
        
    }
    for (int i = 0; i < Xsol.size(); ++i) {
        for (int j = 0; j < Ysol.size(); ++j) {
            if (sol_matrix[i][j] == 1)
                *res++ = make_vector( Xsol[i], Ysol[j] );
        }
    }
    return res;
}




//--------------------------------------------------------------------
__ALGEBRAIC_END_NAMESPACE
//--------------------------------------------------------------------
#endif // SYNAPS_ARITHM_ALGEBRAIC_BSOLVER_H