solver_binary.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_binary_hpp
#define realroot_solver_binary_hpp
/********************************************************************/
#include <realroot/Seq.hpp>
#include <realroot/binomials.hpp>
#include <realroot/univariate_bounds.hpp>
#include <realroot/rounding_mode.hpp>
#include <realroot/numerics_hdwi.hpp>
#include <realroot/sign_variation.hpp>
#include <realroot/solver.hpp>
#include <realroot/homography.hpp>
#include <realroot/scalar_bigunsigned.hpp>
#include <assert.h>
#include <vector>
#include <algorithm>
//====================================================================
namespace mmx {

#ifndef WITH_AS
#define WITH_AS
template<typename T,typename F>
struct cast_helper {
    static inline T cv (const F& x) { return x; }
};

template<typename T,typename F> inline
T as (const F& x) {  return cast_helper<T,F>::cv (x); }
#endif

//  template<class C> inline void operator<<(std::vector<C>& v, const C& x) { v.push_back(x) ; }

template<class C, int N> struct Interval;
//  template<class C> struct with_rounding { typedef meta::false_t T;};
//  template<class C, int N> struct with_rounding<Interval<C,N> > { typedef typename with_rounding<C>::T T;};
//  template<>               struct with_rounding<double> { typedef meta::true_t T;};

//  template<class real>  struct Homography
//   {
//     real a,b,c,d;
//     Homography(): a(1),b(0),c(0),d(1) {}
//     Homography(const Homography& H): a(H.a),b(H.b),c(H.c),d(H.d) {}
//     Homography(const real& A, const real& B, const real& C, const real&D): a(A),b(B),c(C),d(D) {}
//   };

//====================================================================
struct bound
{
    unsigned  e;
    unsigned  m;

    //     bound() :e(0), m(0){}
    //     bound(const bound& b) :e(b.e), m(b_m.hpp){}
    //     bound & operator=(const bound& b) {e=b.e; m=b_m.hpp;return *this;}

    bool operator==( const bound& b )
    {
        if ( e < b.e )
        {
            return (m << (b.e-e)) == b.m;
        };
        return (b.m << (e-b.e)) == m;
    };
};
//--------------------------------------------------------------------
//    struct dmn_t
//    {
//      unsigned   h;
//      unsigned_t v;
//    };
//====================================================================
struct res_t {
    bound a;
    bound b;
    bool  t;
    int   sv;
};
//====================================================================
struct data_t {
    typedef double creal_t;
    typedef unsigned sz_t;
    typedef unsigned unsigned_t;

    creal_t  *m_data;
    bound    *m_dmn ;
    sz_t      m_limit;
    sz_t      m_s;
    sz_t      m_cup;
    sz_t      m_cdw;
    creal_t   eps, root_bound;
    bool      isole;
    std::vector<res_t> m_res;

    data_t(): eps(1e-6) { m_data = 0; };
    data_t(bool isl): eps(1e-6), isole(isl) { m_data = 0; };
    data_t(const creal_t& e): eps(e) { m_data = 0; };
    ~data_t() { Free(); };
    
    inline
    void Free()
    {
        if ( m_data )
        {
            delete[] m_data;
            delete[] m_dmn;
            m_res.resize(0);
        };
    }

    inline
    void alloc( sz_t s, sz_t cs, sz_t deep  )
    {
        Free();
        m_s = s;
        m_limit = std::min((sz_t)numerics::hdwi<sz_t>::nbit,deep);
        m_data  = new creal_t [ cs*m_limit ];
        m_dmn   = new bound   [ m_limit ];
        m_res.resize(0);
        //      std::cout<<"Alloc "<<m_limit<<std::endl;
    }

    inline const bound& bcka() const { return m_res.back().a; };
    inline       bound& bcka()       { return m_res.back().a; };
    inline const bound& bckb() const { return m_res.back().b; };
    inline       bound& bckb()       { return m_res.back().b; };
    
    inline        bool& bckt()       { return m_res.back().t; };
    inline        bool  bckt() const { return m_res.back().t; };
    
    inline        int & bcks()       { return m_res.back().sv; };
    inline        int   bcks() const { return m_res.back().sv; };

    inline void set_back_a(const bound& a) { m_res.back().a = a;};
    inline void set_back_b(const bound& b) { m_res.back().b = b;};
    
    inline void setbck( const bound& a, const bound& b, bool type )
    {
        bcka() = a;
        bckb() = b;
        bckt() = type;
    };

    inline
    void pshbck()
    {
        m_res.push_back(res_t());
    };
    
    inline
    void sstore( int d, int t = 1 )
    {
        // m_data + d*m_s;
        pshbck();
        bcka() = m_dmn[d];
        bckb() = bcka();
        bckb().m += t;
        bckt() = true;
        bcks() = 1;
    };

    inline
    void bstore( int d, int t = 1 )
    {
        // m_data + d*m_s;
        pshbck();
        bcka() = m_dmn[d];
        bcka().m += t;
        bckb() = bcka();
        bckt() = true;
        bcks() = 1;
    };

    inline
    void estore( int d, int t = 1 )
    {
        // m_data + d*m_s;
        pshbck();
        bckb() = m_dmn[d];
        bckb().m += t;
        bcka() = bckb();
        bckt() = true;
        bcks() = 1;
    };

    inline
    void mstore( int d , int sv)
    {
        sstore(d);
        bckt() = false;
        bcks() = sv;
    };

    sz_t size() const  { return m_res.size(); };
    sz_t nb_sol()      { return m_res.size();}


    void writedomain( int d )
    {
        std::cout << "Domaine-" << d << std::endl;
        std::cout << "\thauteur = " << m_dmn[d].e << std::endl;
        std::cout << "\tvaleur  = " << m_dmn[d].m << std::endl;
    };

    void writebck()
    {
        std::cout << "Back " << std::endl;
        std::cout << "\tA = " << bcka().e << " , " << bcka().m  << std::endl;
        std::cout << "\tB = " << bckb().e << " , " << bckb().m  << std::endl;
    };
    
    inline creal_t get_root_bound() { return root_bound;}

    template<class real> static
    void convert( real& a, const bound& b,
                  const real & first = 0, const real & scale = 1  )
    {
        unsigned m(b.m);
        numerics::hdwi<unsigned_t>::reverse(b.e+1,m);
        a = first;
        real x = scale;
        for ( unsigned i = 0; i < b.e + 1; x /= 2, i ++ )
        {
            if ( m & 1 ) a += x;
            m >>= 1;
        };
    }

    template<class real, class coeff>
    void convert( real& a, const bound& b, const homography<coeff>& H)
    {
        unsigned_t m(b.m);
        numerics::hdwi<unsigned_t>::reverse(b.e+1,m);
        a = H.b;
        real d = -H.d;
        real x = H.a, y = -H.c;
        for ( unsigned i = 0; i < b.e + 1; x /= 2, y/=2, i ++ )
        {
            if ( m & 1 ) {a += x; d+=y;}
            m >>= 1;
        };
        //      std::cout <<"\n\t d= "<<d<<" a= "<<a<< std::endl;
        if(d==0 )
        {
            //  std::cout <<"\n\t rb= "<<this->root_bound<< std::endl;
            a = real(root_bound);
        }
        else
        {
            a /= real(-d);
            //std::cout <<"\t d!=0 a= "<<a<< std::endl;
        }
    };


};

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

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

    typedef unsigned sz_t;
    //  typedef unsigned unsigned_t;
    static const sz_t bitquo = numerics::bit_resolution<C>::nbit / numerics::hdwi<sz_t>::nbit;
    static const sz_t bitrem = numerics::bit_resolution<C>::nbit % numerics::hdwi<sz_t>::nbit;

    typedef bigunsigned<bitquo+(sz_t)(bitrem!=0)> unsigned_t;
    typedef C creal_t;

    static data_t data;

    inline
    static void split( creal_t * r, creal_t * l, sz_t sz )
    {
        creal_t * er, * p;
        er = r + (sz-1);
        for ( er = r + (sz-1); er != r; er--, l++ )
        {
            *l = *r;
            for ( p = r; p != er; p ++ )
                *p = (*p+*(p+1))/2;
        };
        *l = *r;
    };

    static sz_t sgncnt ( creal_t const *  b , sz_t sz )
    {
        creal_t const * last = b + sz;
        sz_t c;
        int prv = (((*b>0)<<1)|(*b<0));
        int crr;
        for ( c = 0; b != last && c<2; crr = (((*b>0)<<1)|(*b<0)), c += prv != crr, prv = crr, b ++ ) ;
        //      std::cout<<"sgncnt "<< c<<std::endl;
        return c;
    };

    static inline void split( bound& r, bound& l )
    {
        r.m <<= 1;
        l.m = r.m;
        r.m |=  1;
        r.e ++;
        l.e = r.e;
    };

    static inline
    void print( creal_t * p, sz_t  n )
    {
        std::cout << "[";
        for ( sz_t i = 0; i < n-1; i ++ )
            std::cout << p[i] << ",";
        std::cout << p[n-1] << " ]";
    };
    
    void Loop( bool isole = true )
    {
        //      std::cout <<"Starting simple solver loop"<<std::endl;
        creal_t * pup, * pdw;
        int d;
        sz_t a,s,cup;
        unsigned_t v;
        s = data.m_s;
        pup = data.m_data;
        data.m_dmn[0].m = 0;
        data.m_dmn[0].e = 0;
        for( a = 0, d = 0; d >= 0; d--, a -= s   )
        {
            cup = sgncnt(pup+a,s);
            if ( cup )
            {
                if ( isole && cup == 1 ) { data.sstore(d); continue; };
                if ( data.m_dmn[d].e == data.m_limit-1 )
                {
                    if ( cup == 1 )
                    {
                        if ( *(pup+a) != 0 )
                        {
                            if ( *(pup+a+s-1) == 0 )
                                data.sstore(d,0);
                            else
                                data.sstore(d);
                        };
                    }
                    else
                    {
                        data.mstore(d,cup);
                        //std::cout<<"s "<<cup<<" "<<d<<" "<<data.m_limit<<std::endl;
                    }
                    continue;
                };
                split(pup+a,pup+a+s,s);
                split(data.m_dmn[d],data.m_dmn[d+1]);
                d += 2;
                a += 2*s;
            };
        };
        //std::cout <<"End of simple solver loop"<<std::endl;
    };

    //--------------------------------------------------------------------
    template<class output>
    void get_flags( output& o )
    {
        unsigned s = o.size();
        o.resize( s + o.size() );
        for ( sz_t i = 0; i < data.m_res.size(); o[s+i] = data.m_res[i++].t )
            ;
    };
    //--------------------------------------------------------------------
    template<class input>
    void run_loop( const input& in, const creal_t& eps, bool isole, texp::false_t)
    {
        sz_t prec = numerics::bitprec(eps);
        data.alloc(in.size(),prec);
        std::copy(in.begin(),in.end(),data.m_data);
        Loop(isole);
    };

};

template <class K>  data_t binary_subdivision<K>::data;
//====================================================================
} //namespace mmx
#endif //realroot_solver_binary_hpp