solver_uv_sleeve.hpp

Go to the documentation of this file.

/*******************************************************************
   *   This file is part of the source code of the realroot kernel.
   *   Author(s): J.P. Pavone GALAAD, INRIA
   *              B. Mourrain GALAAD, INRIA
   ********************************************************************/
#ifndef realroot_solver_sleeve_hpp
#define realroot_solver_sleeve_hpp
/********************************************************************/
//====================================================================
# include <realroot/IEEE754.hpp>
# include <realroot/texp_bool.hpp>
# include <realroot/array.hpp>
# include <realroot/solver.hpp>
# include <realroot/solver_binary.hpp>
# include <realroot/polynomial_bernstein.hpp>
# include <realroot/polynomial_tensor.hpp>

# define TMPL    template<class C, class V>
# define SOLVER  solver<C, Sleeve<V> >
//====================================================================
namespace mmx {
template<class V> struct Sleeve{};

template<class C>
struct sleeve_rep {
    C *up, *dw;
    unsigned m_sz;

    sleeve_rep(unsigned n): up(new C[n]), dw(new C[n]), m_sz(n) {};
    unsigned size() const { return m_sz; }
};

struct Monomials;
struct Bernstein;
namespace let {
template<class C, class D> void
assign(sleeve_rep<C>& f, const polynomial<D, with<Bernstein> >& p) {
    assert(f.size()==(unsigned)(degree(p)+1));
    {
        numerics::rounding<C> rnd(numerics::rnd_up());
        for(unsigned i=0;i<f.size();i++) f.up[i]=as<C>(p[i]);
    }
    {
        numerics::rounding<C> rnd(numerics::rnd_dw());
        for(unsigned i=0;i<f.size();i++) f.dw[i]=as<C>(p[i]);
    }
}
template<class C, class D> void
assign(sleeve_rep<C>& f, const polynomial<D, with<MonomialTensor> >& p) {
    assert(f.size()==(unsigned)(degree(p)+1));
    typedef typename kernelof<D>::T::integer  integer;
    typedef typename kernelof<D>::T::rational rational;

    binomials<integer> bn;
    unsigned ip=0, in=0;
    int d=degree(p);
    for(unsigned i=0;i< p.size();i++) {
        if(p[i]>p[ip]) {
            ip = i;
        }
        if(p[i]<p[in]) in=i;
    }
    rational mx = (p[ip]>-p[in]?as<rational>(p[ip]):as<rational>(-p[in]));
    rational c;
    {
        numerics::rounding<C> rnd( numerics::rnd_up() );
        for ( int i = 0; i <= d; i ++ ) {
            c = as<rational>(p[i])/(mx*bn(d,i));
            f.up[i]=as<C>(c);
        }
    }
    {
        numerics::rounding<C> rnd( numerics::rnd_dw() );
        for ( int i = 0; i <= d; i ++ ) {
            c = as<rational>(p[i])/(mx*bn(d,i));
            f.dw[i]=as<C>(c);
        }
    }
    //      polynom< ring<rational,Bernstein> > q;
    //      let::assign(q,p);
    //  assign(f,q);
}

}
//====================================================================
template<class K, class Num>  struct binary_convert {};
template<class K>
struct binary_convert<K,Isolate>: public K {

    typedef Seq<Interval<typename K::ieee> >  Solutions;
    static data_t data;

    template<class output, class real> static inline
    void get(output& sol, const real& first = 0, const real& last =1)
    {
        typedef typename output::value_type root_t;
        //typedef root_t                      as_root;
        //      unsigned s = o.size();
        real a,b;
        //      o.resize(s+2*size());
        real scale = last-first;
        for ( unsigned i = 0; i < data.nb_sol(); i ++ )
        {
            data_t::convert(a,data.m_res[i].a,first, scale);
            data_t::convert(b,data.m_res[i].b,first, scale);
            sol << as<root_t>(Interval<real>(a,b));
        };
    };

    template<class output, class real> static inline
    void get(output & sol, const homography<real>& H)
    {
        typedef typename output::value_type root_t;
        typedef root_t                      as_root;
        typedef typename root_t::value_type bound_t;
        //      unsigned s = o.size();
        bound_t u,v;
        real dt= H.a*H.d-H.b*H.c;
        if(dt>0)
            for (unsigned i = 0; i < data.nb_sol(); i ++ )
            {
                data.convert(u,data.m_res[i].a,H);
                data.convert(v,data.m_res[i].b,H);
                sol << as_root(u,v);
            }
        else
            for ( int i = data.nb_sol()-1; i >=0 ; i --)
            {
                data.convert(u,data.m_res[i].b,H);
                data.convert(v,data.m_res[i].a,H);
                sol << as_root(u,v);
            }
    }
};
template <class K>  data_t binary_convert<K,Isolate>::data(true);
//====================================================================
template<class K>
struct binary_convert<K,Approximate> : public K
{
    typedef unsigned unsigned_t;
    typedef Seq<typename K::ieee>  Solutions;
    typedef typename K::integer integer;
    static data_t data;

    static inline
    integer binomial(const integer& n,const integer& p){return K::binomial(n,p);}

    template<class C> static inline C as_root(const C& a, const C& b) {return (a+b)/2;}

    template<class output, class real> static inline
    void get(output& sol, const real& first = 0, const real& last =1)
    {
        typedef typename output::value_type root;
        real a,b;
        real scale= last-first;
        for (unsigned i = 0; i < data.nb_sol(); i ++ )
        {
            data_t::convert(a,data.m_res[i].a, first, scale);
            data_t::convert(b,data.m_res[i].b, first, scale);
            sol << as<root>(Interval<real>(a,b));
        }
    }

    template<class output, class real> static inline
    void get(output& sol, std::vector<int>& mults, const real& first = 0, const real& last =1)
    {
        get(sol, first, last);
        for (unsigned i = 0; i < data.nb_sol(); i ++ )
            mults.push_back(data.m_res[i].sv);
    }

    template<class output, class real> static inline
    void get(output& sol, const homography<real>& H)
    {
        //  typedef typename output::value_type as_root;
        //root_t u,v;
        unsigned s=data.nb_sol();
        std::vector<real> u(s);
        real v;
        real dt=H.a*H.d-H.b*H.c;

        if(dt>0)
        {
            {
                numerics::rounding<real> rnd(numerics::rnd_dw());
                for ( unsigned i = 0; i < s; i ++ )
                    data.convert(u[i],data.m_res[i].a,H);
            }
            {
                numerics::rounding<real> rnd(numerics::rnd_up());
                for ( unsigned i = 0; i < s; i ++ ) {
                    data.convert(v,data.m_res[i].b,H);
                    sol << as_root(u[i],v);
                }
            }
        }
        else
        {
            {
                numerics::rounding<real> rnd(numerics::rnd_dw());
                for (int i = s-1; i >=0; i -- )
                    data.convert(u[i],data.m_res[i].a,H);
            }
            {
                numerics::rounding<real> rnd(numerics::rnd_up());
                for (int i = s-1; i >=0; i -- ) {
                    data.convert(v,data.m_res[i].b,H);
                    sol << as_root(v,u[i]);
                }
            }
        }
    }

};

template <class K>  data_t binary_convert<K,Approximate>::data(false);

//====================================================================
template < class K>
struct binary_sleeve_subdivision : public K {

    typedef typename K::integer  integer;
    typedef typename K::rational rational;
    //    typedef typename K::floating floating;
    //    typedef typename K::ieee     ieee;

    typedef double                      creal_t;
    typedef unsigned                    sz_t;
    typedef binary_subdivision<K>       Base_t;
    typedef typename Base_t::unsigned_t unsigned_t;
    typedef res_t                       Domain_t;
    
    //    sleeve(): Base_t(1e-6){}
    //    sleeve(creal_t e): Base_t(e){}

    static inline void alloc( sz_t s, sz_t deep  )
    {
        K::data.alloc(s,2*s,deep);
    };

    static void barre( char c, unsigned n )
    {
        char  _bar[ n+1 ];
        std::fill(_bar,_bar+n, c);
        _bar[n] = 0;
        std::cout << _bar << std::endl;
    };

    static void writebounds( creal_t * pup, creal_t * pdw, unsigned s )
    {
        Base_t::print(pup,s); std::cout << std::endl;
        Base_t::print(pdw,s); std::cout << std::endl;
    };

    static inline
    bool glue( sz_t cup, sz_t cdw, int d )
    {
        if ( K::data.bckb() == K::data.m_dmn[d] )
        {
            K::data.bckb() = K::data.m_dmn[d];
            K::data.bckb().m += 1;
            if ( cup ) K::data.m_cup = cup;
            if ( cdw ) K::data.m_cdw = cdw;
            return true;
        };
        return false;
    };

    static inline
    void mstore( sz_t cup, sz_t cdw, int d  )
    {
        if ( K::data.m_res.size() == 0 || !glue(cup,cdw,d)  )
        {
            K::data.mstore(d, std::max(cup,cdw));
            K::data.m_cup = cup;
            K::data.m_cdw = cdw;
            //std::cout << "Push ";
        };
        //      writebck();
    };
    

    static inline
    void dwsplit( creal_t * r, creal_t * l, sz_t s )
    {
        numerics::rounding<creal_t> rnd( numerics::rnd_dw() );
        Base_t::split(r,l,s);
    };

    static inline
    void upsplit( creal_t * r, creal_t * l, sz_t s )
    {
        numerics::rounding<creal_t> rnd( numerics::rnd_up() );
        Base_t::split(r,l,s);
    };

    static void Loop( bool isole = true ) {
        //      std::cout << "LOOP\n";
        numerics::rounding<creal_t> rnd(numerics::rnd_dw());
        creal_t * pup, * pdw;
        int d;
        sz_t a,s,cup,cdw;
        //      unsigned_t v;
        s = K::data.m_s;
        pup = K::data.m_data;
        pdw = K::data.m_data + s;
        K::data.m_dmn[0].m = 0;
        K::data.m_dmn[0].e = 0;
        for( a = 0, d = 0; d >= 0; d--, a -= 2*s   )
        {
            cup = Base_t::sgncnt(pup+a,s);
            cdw = Base_t::sgncnt(pdw+a,s);
            if ( cup || cdw )
            {
                if ( isole && cup == cdw && cup == 1 ){
                    //std::cout << " (b) \n";
                    K::data.sstore(d);
                    continue;
                };
                if ( K::data.m_dmn[d].e == K::data.m_limit-1 ) {
                    // std::cout << " (c) "<<cup<<" "<<cdw<<"\n";
                    if ( cup == cdw && cup == 1 )
                        K::data.sstore(d);
                    else
                        mstore(cup,cdw,d);
                    continue;
                };
                dwsplit(pdw+a,pdw+a+2*s,s);
                upsplit(pup+a,pup+a+2*s,s);
                Base_t::split(K::data.m_dmn[d],K::data.m_dmn[d+1]);
                d += 2;
                a += 4*s;
                continue;
            } else {
                //std::cout<<"  (a) ";
                sz_t k =0, sv=0;
                for (  k = 0; k < K::data.m_s; k ++ ) {
                    if ( ( pup[a+k] > 0 ) != ( pdw[a+k] > 0 ) ) {
                        sv++;
                        if(sv>1) {
                            mstore(cup,cdw,d);
                            break;
                        };
                    }
                }
                //std::cout<<std::endl;
            };
        };
        // std::cout << "END LOOP\n";
    };


    template<class input> static
    void run_loop( const input& in, const creal_t& eps, const texp::true_t&   ) {
        //      binary_subdivision<K> slv(eps);
        sz_t prec = numerics::bitprec(eps);
        alloc(in.size(),prec);
        std::copy(in.begin(),in.end(),K::data.m_data);
        std::copy(in.begin(),in.end(),K::data.m_data+K::data.m_s);
        Loop(K::data.isole);
    };

    const Domain_t& operator[]( int i ) const  { return K::data.m_res[i]; };
    Domain_t&       operator[]( int i ) { return K::data.m_res[i]; };

    
    template<class output, class input> static inline
    void solve_bernstein( output& out, const input& in) {
        run(in,K::data.eps);
        //      solutions<typename output::value_type>::get(this,out);
    }
    template<class VECT,class POL,class Q> static inline
    void init_pol(VECT& ubp, VECT& dbp, const POL& r,  unsigned sz, const Q& u, const Q&v);


    template<class input> static inline
    void run( const input& in, const creal_t& eps)
    {
        typedef typename numerics::is_rounded<creal_t>::result_t round_t;
        run_loop(in,eps,round_t());
    }

    template<class C> static void
    run( const sleeve_rep<C>& p) {
        sz_t prec = numerics::bitprec(K::data.eps);
        alloc(p.size(),prec);
        std::copy(p.up,p.up+p.size(),K::data.m_data);
        std::copy(p.dw,p.dw+p.size(),K::data.m_data+K::data.m_s);
        Loop(K::data.isole);
    }

    template<class input> static
    void run( const input& up, const input& dw) {
        assert(up.size()==dw.size());
        sz_t prec = numerics::bitprec(K::data.eps);
        alloc(up.size(),prec);
        std::copy(up.begin(),up.end(),K::data.m_data);
        std::copy(dw.begin(),dw.end(),K::data.m_data+K::data.m_s);
        Loop(K::data.isole);
    }

    template <typename output, typename POL, typename Q> static
    void solve(output& sol, const POL& r, const Q& u, const Q& v);

    //--------------------------------------------------------------------
    template<class output, class input,class real,class MTH> static inline
    void solve_bernstein( output& sol,
                          const input& up, const input& dw,
                          const real& u, const real& v, const MTH& mth)
    {
        K::data.isole=mth.isole;
        run(up,dw);
        get(sol,u,v);
        //      solutions<typename output::value_type>::get(this,out,u,v);
    }
    
    //--------------------------------------------------------------------
    template <typename output, typename POL ,typename real, typename MTH> static
    void solve_bernstein(output& sol, const POL& r, const real& u,const real& v, const MTH& mth )
    {
        //typedef typename output::value_type    root_t;
        using let::assign;
        
        unsigned sz= r.size();
        std::vector<creal_t> ubp( sz ), dbp( sz );
        {
            numerics::rounding<creal_t> rnd( numerics::rnd_up() );
            for(unsigned i=0;i<sz;i++) ubp[i] =as<creal_t>(r[i]);
        }
        {
            numerics::rounding<creal_t> rnd( numerics::rnd_dw() );
            for(unsigned i=0;i<sz;i++) dbp[i] =as<creal_t>(r[i]);
        }
        K::data.isole=mth.isole;
        run(ubp,dbp);
        get(sol,u,v);
        //      solutions<typename output::value_type>::get(this, sol, u,v);
    }
};
//====================================================================
template<class K>
template<class VECT,class POL,class Q>
void binary_sleeve_subdivision<K>::init_pol(VECT& ubp, VECT& dbp, const POL& r, unsigned sz,const Q& u, const Q&v)
{
    using let::assign;
    
#if 1
    //    PRINT_DEBUG("Conversion using exact arithmetic");
    rational U,V;
    assign(U,u);
    assign(V,v);
    //    std::vector<interval> bp( sz );
    std::vector<rational> bp(sz);
    univariate::convertm2b(bp,r,sz,U,V);
    rational mx = array::max_abs(bp);
    array::div(bp,mx);
    // unsigned i;
    //    for (i = 0; i < sz && bp[i].upper()-bp[i].lower()<as<floating>(1.e-16); i ++ );
    //    if(i<sz) std::cout<<"refine: "<<i<<" " <<(bp[i].upper()-bp[i].lower())<<std::endl;
    //    double l; assign(l, rational(V-U));
    //    std::cout<<"\t init: "<<sign_variation(bp)<<" "<<l<<std::endl;
    {
        numerics::rounding<creal_t> rnd( numerics::rnd_up() );
        for ( unsigned i = 0; i < sz; i ++ ) ubp[i]=as<creal_t>(bp[i]);
    }
    {
        numerics::rounding<creal_t> rnd( numerics::rnd_dw() );
        for ( unsigned i = 0; i < sz; i ++ ) dbp[i]=as<creal_t>(bp[i]);
    };

#else
    //    PRINT_DEBUG("Using rounding arithmetic");
    O ud, vd;
    {
        numerics::rounding<creal_t> rnd( numerics::rnd_dw() );
        ud=as<creal_t>(u);
    };
    {
        numerics::rounding<creal_t> rnd( numerics::rnd_up() );
        vd=as<creal_t>(v);
        BEZIER::monomial_to_bezier(ubp,r,r.size(),ud,vd);
    };
    {
        numerics::rounding<creal_t> rnd( numerics::rnd_dw() );
        BEZIER::monomial_to_bezier(dbp,r,r.size(),ud,vd);
    };

    for ( unsigned i= 0; i < ubp.size(); i ++ )
        if ( dbp[i] > ubp[i] ) {
            std::swap(dbp[i], ubp[i]);
        };
#endif

}
//------------------------------------------------------------------
template <class K>
template <typename output, typename POL ,typename Q>
void binary_sleeve_subdivision<K>::solve(output& sol, const POL& r, const Q& u, const Q& v)
{
    unsigned sz= r.size();
    std::vector<double> ubp(sz), dbp(sz);
    init_pol(ubp,dbp,r,sz,u,v);
    run(ubp,dbp,true);
    //    run(ubp,dbp,this->eps,mth.isolate);
    //    solutions<typename output::value_type>::get(this,sol,u,v);
}

//====================================================================
TMPL
struct solver<C, Sleeve<V> > {
    
    typedef typename kernelof<C>::T          K;
    typedef typename binary_convert<K,V>::Solutions         Solutions;
    // Seq<Interval<typename K::ieee> >                    Solutions;
    typedef binary_sleeve_subdivision<binary_convert<K,V> > base_t;


    template<class POL> static void
    solve(Solutions& sol, const POL& p);

    template<class POL,class T> static void
    solve(Solutions& sol, const POL& p, const T& u, const T& v);

    template<class T, class U, class VRT> static void
    solve(Solutions& sol, const polynomial<U,VRT>& pl, const polynomial<U,VRT>& pu, const T& u, const T& v);

    template<class T, class U, class VRT> static void
    solve(Solutions& sol, std::vector<int>& mult, const polynomial<U,VRT>& pl, const polynomial<U,VRT>& pu, const T& u, const T& v);

    template<class U,class VRT> static void
    solve(Solutions& sol, const polynomial<U,VRT>& pl, const polynomial<U,VRT>& pu);


    static inline void
    set_precision(unsigned p) { binary_convert<K,V>::set_precision(p); }

};

TMPL
template<class POL> void
        SOLVER::solve(Solutions& sol, const POL& p) {
    typedef typename K::ieee real_t;
    sleeve_rep<real_t> f(degree(p)+1);
    let::assign(f,p);
    real_t c;
    for(unsigned i=1;i< f.size();i+=2) {
        c=f.up[i];f.up[i]=-f.dw[i];f.dw[i]=-c;
    }
    base_t::run(f);

    homography<real_t> H(1,0,1,-1);
    base_t::data.root_bound = - bound_root(p,Cauchy<real_t>());
    base_t::get(sol,H);

    for(unsigned i=1;i< f.size();i+=2) {
        c=f.up[i];f.up[i]=-f.dw[i];f.dw[i]=-c;
    }
    base_t::run(f);
    base_t::data.root_bound = - base_t::data.root_bound;
    H.c=-1;H.d=1;
    base_t::get(sol,H);
}
//--------------------------------------------------------------------
TMPL
template<class T, class U, class VRT> void
        SOLVER::solve(Solutions& sol,
                      const polynomial<U,VRT>& pl, const polynomial<U,VRT>& pu,
                      const T& u, const T& v) {
    base_t::run(pl,pu);
    base_t::get(sol,u,v);
}

//--------------------------------------------------------------------
TMPL
template<class T, class U, class VRT> void
        SOLVER::solve(Solutions& sol, std::vector<int>& mult,
                      const polynomial<U,VRT>& pl, const polynomial<U,VRT>& pu,
                      const T& u, const T& v) {
    base_t::run(pl,pu);
    base_t::get(sol,mult,u,v);
}
//--------------------------------------------------------------------
TMPL
template<class U, class VRT> void
        SOLVER::solve(Solutions& sol,
                      const polynomial<U,VRT>& pl, const polynomial<U,VRT>& pu) {
    base_t::run(pl,pu);
    base_t::get(sol,0,1);
}

//--------------------------------------------------------------------
TMPL
template<class POL, class T> void
        SOLVER::solve(Solutions& sol, const POL& p,
                      const T& u, const T& v) {
    //typedef typename POL::Ring                        Ring;
    typedef typename SOLVER::base_t      base_t;
    //typedef typename SOLVER::Solutions Solutions;
    unsigned sz= p.size();
    std::vector<double> ubp(sz), dbp(sz);
    base_t::init_pol(ubp,dbp,p,sz,u,v);
    base_t::run(ubp,dbp);
    base_t::get(sol,u,v);
}


//--------------------------------------------------------------------
template<class POL, class M, class B>
typename solver<B, Sleeve<M> >::Solutions
solve (const POL& p, const Sleeve<M>& slv, const B& u, const B& v) {
    typedef solver<B, Sleeve<M> >                       Solver;
    typedef typename Solver::Solutions                  Solutions;
    Solutions sol;
    Solver::solve(sol,p, u,v);
    return sol;
}
//--x------------------------------------------------------------------
template<class POL, class M>
typename solver<typename POL::Scalar, Sleeve<M> >::Solutions
solve (const POL& p, const Sleeve<M>& slv) {
    typedef typename POL::Scalar                             Scalar;
    typedef solver<Scalar, Sleeve<M> >                       Solver;
    typedef typename solver<Scalar, Sleeve<M> >::Solutions   Solutions;
    Solutions sol;
    Solver::solve(sol,p);
    return sol;
}

} //namespace mmx
//====================================================================
# undef TMPL
# undef SOLVER
//====================================================================
#endif  //realroot_solver_sleeve_hpp