synaps/usolve/sturmsmdeg.h

// 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_ASOLVER_H
#define SYNAPS_ARITHM_ALGEBRAIC_ASOLVER_H

#include <synaps/usolve/Algebraic_config.h>
#include <synaps/usolve.h>
#include <synaps/usolve/Method_name.h>
#include <synaps/usolve/Construct_rootof.h>
#include <synaps/upol/SturmSeq.h>
#include <synaps/upol/SquareFree.h>
#include <synaps/usolve/Util.h>
#include <synaps/upol/Bound.h>
#include <synaps/usolve/method_name.h>
#include <stack>


__ALGEBRAIC_BEGIN_NAMESPACE

#if 0

//--------------------------------------------------------------------------------

template < typename RT, typename POLY >
void
compute_intermediate_points(const Seq< root_of<RT, POLY> >& V,
                            Seq< typename NumberTraits<RT>::FT >& Points)
{
    typedef root_of<RT, POLY>    RO_t;
    typedef NumberTraits<RT>    NTR;
    typedef typename NTR::FT     FT;
    typedef typename NTR::FIT    FIT;

    typedef typename Seq<RO_t>::const_iterator  IT1;

    IT1 it = V.begin();
    Points.push_back( FT(it->left() - 1) );

    ++it;
    for (; it != V.end(); ++it) {
        while ( ((it-1)->right() >= it->left())  ) {
            (it-1)->refine_interval();
            it->refine_interval();
        }
        Points.push_back( FT(((it - 1)->right() + it->left())/2) );
    }
    Points.push_back(FT((it-1)->right() + 1));
    return;
}

//--------------------------------------------------------------------------------


//--------------------------------------------------------------------------------
//--------------------------------------------------------------------------------
//--------------------------------------------------------------------------------


template < typename RT,
           typename OutputIterator >
OutputIterator
solve( const typename Sturm<RT>::Poly& F,
       OutputIterator  sol,
       unsigned mult,
       const Sturm<RT>& mth)
{
    typedef typename Sturm<RT>::FT    FT;
    typedef typename Sturm<RT>::FIT   FIT;
    typedef typename Sturm<RT>::RO_t  RO_t;
    typedef typename Sturm<RT>::Poly  POLY;

    typedef Quadtuple<FT>   Quad;

    std::stack<Quad>   Work;
    std::vector<FIT>   result;

    SturmSeq<RT> sth(F, UPOLDAR::diff(F), SUBRES());

    FT b = max_bound(sth[0], Cauchy());
    FT a = -b;
    POLY h = sth[0];

    int Va = Var(sth, a);
    int Vb = Var(sth, b);
    int V = Va - Vb;

    if (V == 0) return sol;
    if (V == 1) {
        *sol++ = RO_t(sth[0], a, b, mult);
        return sol;
    }
    Work.push( Quad(a, b, Va, Vb) );

    POLY x(2, AsSize()); x[1] = RT(1);

    while (!Work.empty()) {
        Quad q = Work.top();
        Work.pop();
    
        FT c = (q.a() + q.b())/2;
        if (sign(UPOLDAR::eval_horner(sth[0], c)) == ZERO) {
            result.push_back( FIT(c, c) );

            POLY t(2, AsSize());
            t[1] = denominator(c);
            t[0] = -numerator(c);
            h = EUCLIDIAN::pquo(h, t);

            POLY f = translate(sth[0], c)/x;
            FT M = min_bound(f, Cauchy());
            //int Vc_plus = Count(sth, FT(c + M));
            int Vc_plus = Var(sth, FT(c + M));
            //int Vc_minus = Var(sth, FT(c - M));
            int Vc_minus = Vc_plus + 1;

            if (q.Va() == Vc_minus + 1) {
                result.push_back( FIT(q.a(), FT(c - M)) );
                //continue;
            }
            if (q.Va() > Vc_minus + 1) {
                Work.push( Quad(q.a(), FT(c - M), q.Va(), Vc_minus) );
                //continue;
            }
            if (q.Vb() == Vc_plus - 1) {
                FIT I(FT(c + M), q.b());
                result.push_back( FIT(FT(c + M), q.b()) );
                //continue;
            }
            if (q.Vb() < Vc_plus - 1) {
                Work.push( Quad(c + M, q.b(), Vc_plus, q.Vb()) );
                //continue;
            }
        } else {
            int Vc = Var(sth, c);
            if (q.Va() == Vc + 1) {
                result.push_back( FIT(q.a(), c) );
                //continue;
            }
            if (q.Vb() == Vc - 1) {
                result.push_back( FIT(c, q.b()) );
                //continue;
            }

            if (q.Vb() < Vc - 1) {
                Work.push( Quad(c, q.b(), Vc, q.Vb()) );
                //continue;
            }
            if (q.Va() > Vc + 1) {
                Work.push( Quad(q.a(), c, q.Va(), Vc) );
                //continue;
            }
        }

    }
    std::sort(result.begin(), result.end());
    typename std::vector<FIT>::const_iterator it;
    for (it = result.begin(); it != result.end(); ++it) {
        if (singleton(*it) ) {
            RO_t r = make_root_of_1(RT(-numerator(it->lower())), denominator(it->lower()), mth);
            r.set_multiplicity(mult);
            *sol++ = r;
        } else {
            *sol++ = RO_t(h, it->lower(), it->upper(), mult);
        }
    }

    return sol;
}




//----------------------------------------------------------------------------


template < typename RT,
           typename OutputIterator >
OutputIterator
sturm_in_interval( const SturmSeq<RT>                  StHa,
                   const typename Sturm<RT>::FIT       I,
                   OutputIterator                      sol,
                   unsigned                            mult,
                   const Sturm<RT>&                    MTH)
{
    typedef typename Sturm<RT>::FT    FT;
    typedef typename Sturm<RT>::FIT   FIT;
    typedef typename Sturm<RT>::Poly  Poly;
    typedef typename Sturm<RT>::RO_t  RO_t;

    typedef Quadtuple<FT>   Quad;

    std::stack<Quad>   Work;
    std::vector<FIT>   result;

    FT a = lower(I);
    FT b = upper(I);

    Poly h = StHa[0];

    int Va = Var(StHa, a);
    int Vb = Var(StHa, b);
    int V = Va - Vb;

    if (V == 0) return sol;
    if (V == 1) {
        *sol++ = RO_t(StHa[0], a, b, mult);
        return sol;
    }
    Work.push( Quad(a, b, Va, Vb) );

    Poly x(2, AsSize()); x[1] = RT(1);

    while (!Work.empty()) {
        Quad q = Work.top();
        Work.pop();
   
        FT c = (q.a() + q.b())/2;
        if (sign(UPOLDAR::eval_horner(StHa[0], c)) == ZERO) {
            result.push_back( FIT(c, c) );

            Poly t(2, AsSize());
            t[1] = denominator(c);
            t[0] = -numerator(c);
            h = EUCLIDIAN::pquo(h, t);

            Poly f = translate(StHa[0], c)/x;
            FT M = min_bound(f, Cauchy());
            int Vc_plus = Var(StHa, FT(c + M));
            //int Vc_minus = Var(StHa, FT(c - M));
            int Vc_minus = Vc_plus + 1;

            if (q.Va() == Vc_minus + 1) { result.push_back( FIT(q.a(), FT(c - M)) ); }
            if (q.Va() > Vc_minus + 1) { Work.push( Quad(q.a(), FT(c - M), q.Va(), Vc_minus) ); }

            if (q.Vb() == Vc_plus - 1) { result.push_back( FIT(FT(c + M), q.b()) ); }
            if (q.Vb() < Vc_plus - 1) { Work.push( Quad(c + M, q.b(), Vc_plus, q.Vb()) ); }
        } else {
            int Vc = Var(StHa, c);
            if (q.Va() == Vc + 1) { result.push_back( FIT(q.a(), c) ); }
            if (q.Vb() == Vc - 1) { result.push_back( FIT(c, q.b()) ); }

            if (q.Vb() < Vc - 1) { Work.push( Quad(c, q.b(), Vc, q.Vb()) ); }
            if (q.Va() > Vc + 1) { Work.push( Quad(q.a(), c, q.Va(), Vc) ); }
        }
    }
    std::sort(result.begin(), result.end());
    typename std::vector<FIT>::const_iterator it;
    for (it = result.begin(); it != result.end(); ++it) {
        if (singleton(*it) ) {
            RO_t r = make_root_of_1(it->lower(), MTH);
            r.set_multiplicity(mult);
            *sol++ = r;
        } else {
            *sol++ = RO_t(h, it->lower(), it->upper(), mult);
        }
    }
    return sol;
}





template < typename NT,
           typename OutputIterator >
OutputIterator
solve( const typename Sturm<NT>::Poly& f,
       OutputIterator  sol,
       const Sturm<NT>& mth)
{
    //std::cout << "solve<Sturm>" << std::endl;
    typedef typename Sturm<NT>::RT    RT;
    typedef typename Sturm<NT>::FT    FT;
    typedef typename Sturm<NT>::FIT   FIT;
    typedef typename Sturm<NT>::RO_t  RO_t;
    typedef typename Sturm<NT>::Poly  Poly;

    if (UPOLDAR::degree(f) < 5) {
        return solve(f, sol, Small_degree<RT>());
    }
  
    Seq<Poly> PL = SquareFreeFactorization(f);

    for (int i = 0; i < PL.size(); ++i) {
        if (UPOLDAR::degree(PL[i]) < 5) {
            //std::cout << "small\n";
            solve(PL[i], sol, i+1, Small_degree<NT>());
        } else {
            //std::cout << "big\n";
            solve(PL[i], sol, i+1, Sturm<NT>());
        }
    }   
    return sol; 
}

#endif


__ALGEBRAIC_END_NAMESPACE


#endif // SYNAPS_ARITHM_ALGEBRAIC_ASOLVER_H