/Users/mourrain/Devel/mmx/realroot/include/realroot/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;  
  };
  //====================================================================
  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 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;
    };
  
    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;
    };

    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;
    };
  
    inline 
    void mstore( int d ) 
    {
      sstore(d);
      bckt() = false;
    };

    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 h,a,s,cup,cdw;
      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);
                      //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