/Users/mourrain/Devel/mmx/realroot/include/realroot/solver_ucf.hpp

Go to the documentation of this file.

/********************************************************************
 *   Author(s): A. Mantzaflaris, GALAAD, INRIA
 ********************************************************************/
#ifndef _realroot_SOLVER_FIRST_ROOT_CONFRAC_H
#define _realroot_SOLVER_FIRST_ROOT_CONFRAC_H
/********************************************************************/

# include <realroot/Interval.hpp>

# include <realroot/homography.hpp>
# include <realroot/sign_variation.hpp>
# include <realroot/univariate_bounds.hpp>
# include <stack>
# include <realroot/solver.hpp>

# define NO_ROOT -1.0
# define TMPL  template<class POL>
# define TMPLR template<class Real, class POL>
# define N_ITER 10000

namespace mmx {

    struct CFfirstIsolate{};
    struct CFfirstApproximate{};
    struct CFfirstFloor{};
    struct CFallIsolate{};
    struct CFseparate{};
    struct CFdecide{};

    //bit_size support for QQ

    // inline long bit_size(QQ q) 
    // {
    //          long n,d;
    //          n= bit_size(  numerator(q) );
    //          d= bit_size(denominator(q) );
    //          return( n>d ? n : d ) ;
    // }


/* ******************************************************** */
/* Represents the polynomial on an interval [u,v].
   Keeps homography information to map [0,inf] back to [u,v].
   Main operations are shift, contraction and reciprocal.
 */
TMPL
class interval_rep
{
    typedef typename POL::scalar_t C;

    POL   f;
    homography<C> hg;

public:

    interval_rep(){ }

    interval_rep(const POL p) 
        {
            f  = p;
            hg = homography<C>() ;
        }

    interval_rep(const POL p, homography<C> h) 
        {
            f  = p;
            hg = h;
        }


    interval_rep(const POL p, C u,  C v ) 
        {
            f  = p;
            hg = homography<C>() ;

            shift_box           ( u );
            contract_box      ( v-u );
            reverse_and_shift_box (1);
            //print();
        }

    /* Access the equation */
    POL inline poly() {return f;};
    /* Access the homography */
    homography<C> hom() { return hg; }

    /* Access the domain */
    template<class Real>
    Interval<Real> domain() 
        { 
            Real l, r;
            Interval<Real>  s ;
            //lim to 0
            if ( hg.b!=0 && hg.d!=0 )
                l= as<Real>(hg.b)/as<Real>(hg.d);
            else if ( hg.d==0 )
                l= 10000000;
            else if ( hg.b==0 )
                l= 0 ;
            else 
                l= as<Real>(hg.a)/as<Real>(hg.c) ;
            //lim to inf
            if ( hg.a!=0 && hg.c!=0 )
                r= as<Real>(hg.a)/as<Real>(hg.c);
            else if ( hg.c==0 )
                r=10000000;
            else if ( hg.a==0 )
                r= 0 ;
            else 
                r= as<Real>(hg.b)/as<Real>(hg.d);
            
            if ( l<=r ) s= Interval<Real>(l,r);
            else        s= Interval<Real>(r,l);
    
            return s; 
        }

    /* Access the size of the domain */
    template<class Real>
    Real inline width() 
        { 
            return this->template domain<Real>().width(); 
        }

    /* Access the point where t\in\RR_+ corresponds */
    template<class Real>
    Real inline point(C & t)
        {
            return ( as<Real>(hg.a)*as<Real>(t) + as<Real>(hg.b) ) /
                   ( as<Real>(hg.c)*as<Real>(t) + as<Real>(hg.d) ) ;
        }
    
    C middle()
        {
            C t(hg.d / hg.c);
            
            if (t>0)
                return ( t ) ;
            else
                return C(1);
        }

    unsigned sign_var()
        {
            return sign_variation(f.begin(),f.end()) ;
        }


    // Shift the polynomial by t
    void shift_box(const C& t)
        {
            int s, k,j, sz= f.size();
            
            for ( s = sz, j = 0; j <= s-2; j++ )
                for( k = s-2; k >= j; k-- ) 
                    f[k] += t*f[k+1];
            
            //update homography
            hg.shift_hom(t);
        };
    
    // x = t*x
    void contract_box(const C & t )
        {
            int s, sz= f.size();
            
            for ( s = 0 ; s < sz; ++s )
                f[s] *= pow(t,s) ;
            
            //update homography
            hg.contract_hom(t);
        };
    
    // x = 1/x   -- parameter t=1 unused for now
    void reverse_box(const C & t=1 )
        {
            unsigned s, l= f.size()-1;
            
            for (s = 0 ; s <= l/2; s++ )
                std::swap(f[s], f[l-s]);
            
            //update homography
            hg.reciprocal_hom(t);
        };
    
    // x = 1/(x+1)   -- parameter t=1 unused for now
    void inline reverse_and_shift_box(const C & t=1)
        {
            reverse_box(t);
            shift_box  (C(1));
        };
    
    void print() 
        {
            std::cout << "----------Interval---------------"    << "\n" ;
            std::cout << poly() << "\n" ;
            std::cout<< "("<<hg.a <<"x + " << hg.b<<")/("<<hg.c<<"x+ "<<hg.d << ")"<<", sv:"<<sign_var() <<std::endl;
            std::cout << domain() << "\n" ;
            std::cout << "-------------------------------"    << "\n" ;
        }
    
}; /* Interval_rep */
    

/* ******************************************************** */
    /* 
       Solver interface
    */
    TMPLR
    struct solver_cffirst {
        
        interval_rep<POL> ir;
        Real eps;
        
        solver_cffirst(interval_rep<POL> p){
            ir = p;
            eps= as<Real>(1e-7);
        } 
        
        interval_rep<POL>           first_root();
        Interval<Real>      first_root_isolate();
        Real            first_root_approximate();
        Real                  first_root_floor();

        Seq<interval_rep<POL> >      all_roots();
        Seq<Interval<Real> > all_roots_isolate();
        Seq<Real>           all_roots_separate();
    };
   
    /* 
       Main routine:
       Returns an "interval_rep" containing the minimum >0 root.
    */
    TMPLR
    interval_rep<POL> solver_cffirst<Real,POL>::first_root()
    {
        typedef interval_rep<POL> BOX;
        std::stack<BOX * > uboxes;
        typedef typename POL::scalar_t C;
        
        POL zero(0);
        Interval<Real> mid(0);
        BOX * b, * tmp;
        POL p;
        unsigned s, iters=0;
        
        Interval<Real> I;
        AkritasBound<C> lb;
        //HongBound<C> lb; 
        //Kioustelidis<C> lb; 
        //NISP<C> lb; 
        //Cauchy<C> lb; 
        
        C lower;
        b= new BOX(ir);
        uboxes.push (b);
        
        while ( !uboxes.empty() && iters++ < N_ITER )
            
        {
            b = uboxes.top() ;
            p = b->poly() ;
            uboxes.pop();
            
            if ( p==zero && (!uboxes.empty()) ) {
                return *b;}
            
            s = b->sign_var() ;
            I = b->template domain<Real>();
            
            if ( s==1 && (!uboxes.empty()) ) {
                return *b;}
            
            if ( s > 0 ) 
            {
                if ( b->template width<Real>() < eps )
                {std::cout<< "first_root: multiple root?"<<b->template domain<Real>() <<std::endl;}
                else
                {
                    lower= lb.lower_bound(p) ;
                    if ( lower!=0 )
                    {
                        if ( p( lower ) == 0 )
                            return *b;
                        
                        tmp= new BOX( *b ) ;
                        free(b);
                        tmp->shift_box( lower );
                        uboxes.push (tmp);
                    }
                    else
                    {
                        // right box
                        tmp = new BOX( *b ) ;
                        tmp->shift_box( 1 );
                        uboxes.push ( tmp );                    
                        //middle if needed
                        if ( p( Real(1) ) == Real(0) )
                        {   tmp = new BOX(zero, b->hom() ) ;
                            uboxes.push ( tmp ); }
                        //left
                        tmp = new BOX( *b ) ;
                        tmp->reverse_and_shift_box( 1 );
                        tmp->reverse_box( 1 );
                        uboxes.push ( tmp );
                        
                        free(b);
                    }
                }
            }
        }
        /* -1 = No positive root */
        return BOX(-1);
    }

    /* 
       Returns a list of "interval_rep"s isolating the positive roots,
       in increasing order.
    */
  TMPLR
  Seq<interval_rep<POL> > solver_cffirst<Real,POL>::all_roots()
  {
    typedef interval_rep<POL> BOX;
    std::stack<BOX * > uboxes;
    typedef typename POL::scalar_t C;

    Real ev(0);
    Seq<interval_rep<POL> > sols;
    POL zero(0);
    Interval<Real> mid(0);
    BOX * b, * tmp;
    POL p;
    unsigned s, iters=0;
        
    Interval<Real> I;
    AkritasBound<C> lb;
    //HongBound<C> lb; 
    //Kioustelidis<C> lb; 
    //NISP<C> lb; 
    //Cauchy<C> lb; 
    
    C lower;
    b= new BOX(ir);
    uboxes.push (b);
    
    while ( !uboxes.empty() && iters++ < N_ITER )
    {
      b = uboxes.top() ;
      p = b->poly() ;
      uboxes.pop();
      
      if ( p==zero && (!uboxes.empty()) ) {
        sols<< *b;}
      
      I = b->template domain<Real>();
      // Interval<Real> ev= p( I ) ;
      // if ( ev.m * ev.M > 0  ) continue;
      s = b->sign_var() ;
      if ( s==1 && (!uboxes.empty()) ) {
        sols << *b; continue; }
      
      if ( s > 0 ) 
      {
        if ( (!uboxes.empty())  && b->template width<Real>() < eps*.1 )
        {std::cout<<
            "all_roots: multiple root?"<<b->template domain<Real>() <<std::endl; 
          std::cout<< b->poly()<<" , "<<b->template width<Real>()<<std::endl;
          std::cout<<"source: "<<ir.poly()<<std::endl;

        }
        else
        {
          lower= lb.lower_bound(p) ;
          if ( lower!=0 )
          {
            
            let::assign(ev,lower);
            if ( p(ev) == 0 )
              sols<< *b;
            
            tmp= new BOX( *b ) ;
            free(b);
            tmp->shift_box( lower );
            uboxes.push (tmp);
          }                     
          else
          {
            // right box
            tmp = new BOX( *b ) ;
            tmp->shift_box( 1 );
            uboxes.push ( tmp );                        
            //middle if needed
            ev=1;
            if ( p( ev ) == Real(0) )
            {   tmp = new BOX(zero, b->hom() ) ;
              uboxes.push ( tmp ); }
            //left
            tmp = new BOX( *b ) ;
            tmp->reverse_and_shift_box( 1 );
            tmp->reverse_box( 1 );
            uboxes.push ( tmp );
            
            free(b);
          }
        }
      }
    }
    return sols;
  }

    /*
       Return isolation box of the first positive root, [-1,-1] if none exists.
     */
    TMPLR
    Interval<Real> solver_cffirst<Real,POL>::first_root_isolate()
    {
        interval_rep<POL> a( first_root() );
        typename POL::scalar_t t(1);
        
        if ( a.poly()== POL(-1) ) 
            return (0);
        else 
            if (a.poly()==POL(0) ) 
                return( a.template point<Real>(t) );
        else
            return ( a.template domain<Real>() );
    }

    TMPLR
    Real solver_cffirst<Real,POL>::first_root_approximate()
    {
        typedef interval_rep<POL> BOX;
        typedef typename POL::scalar_t C;
        BOX * l, * r= new BOX( first_root() ) ;
        C t(1);
        
        if ( r->poly()== POL(-1) ) 
            return (0);
        else 
            if (r->poly()==POL(0) ) 
                return( r->template point<Real>(t) );
        else
        {
            /* Approximate */
            while ( r->template width<Real>() > eps )
            {
                t= r->middle();
                
                if ( r->poly()( t ) == 0 ) return( r->template point<Real>(t) );
                l= new BOX( *r) ;
                l->shift_box( t );
                
                if ( l->sign_var() ) 
                {   
                    free(r);
                    r=l; 
                    continue;
                }
                else 
                {   
                    free(l);
                    r->contract_box(t);
                    r->reverse_and_shift_box(1);
                    r->reverse_box(1);
                }  
            }
            
            return ( r->template domain<Real>() );
        }
    }
    
    TMPLR
    Real solver_cffirst<Real,POL>::first_root_floor()
    {
        typedef interval_rep<POL> BOX;
        typedef typename POL::scalar_t C;

        BOX * l, * r= new BOX( first_root() ) ;
        Interval<Real> I;
        C t(1);
        if ( r->poly()== POL(-1) ) 
            return (0);
        else 
            if (r->poly()==POL(0) ) 
                return( r->template point<Real>(t) );
            else
            {
                
                I= r->template domain<Real>();
                if ( r->poly() == POL(0) ) return as<Real> (floor(as<double>(r->template point<Real>(t))) );
                while ( rfloor(I.m)!=rfloor(I.M) )
                {
                    t= r->middle();
                    if ( r->poly()( t ) == 0 ) return as<Real>( floor(as<double>(r->template point<Real>(t))) ); 
                    
                    l= new BOX( *r) ;
                    l->shift_box( t );
                    
                    if ( l->sign_var() ) 
                    {   
                        free(r);
                        r=l; 
                        I= r->template domain<Real>(); 
                        continue;
                    }
                    else 
                    {   
                        free(l);
                        r->reverse_and_shift_box();
                        r->reverse_box();
                        I= r->template domain<Real>();
                    }  
                }
            }
        return as<Real>( floor(as<double>(I.m)) );
    }/*first_root_floor*/

    TMPLR
    Seq<Interval<Real> > solver_cffirst<Real,POL>::all_roots_isolate()
    {
        Seq<interval_rep<POL> > a( all_roots() );
        Seq<Interval<Real> > s;
        typename POL::scalar_t t(1);

        for ( unsigned i=0; i<a.size(); ++i)
            if (a[i].poly()==POL(0) ) 
                s<< Interval<Real>(a[i].template point<Real>(t) );
            else
                s<< a[i].template domain<Real>();

        return s;
    }

    TMPLR
    Seq<Real> solver_cffirst<Real,POL>::all_roots_separate()
    {
        Seq< Interval<Real> > ints;
        Seq<Real> sep;
        
        ints= all_roots_isolate();
        for (unsigned i=1; i<ints.size(); ++i)
            sep << (ints[i-1].M + ints[i].m)*Real(0.5) ;

        return sep;
    }

//--------------------------------------------------------------------
template<class Ring>
struct solver<Ring, CFfirstIsolate > {

    typedef typename Ring::scalar_t      base_t;

    template<class Real, class POL> static Interval<Real>
    solve(const POL& p,const base_t& u,const base_t& v);

    template<class Real, class POL> static Interval<Real>
    solve(const POL& p);
};

template<class Ring> 
template<class Real, class POL> 
Interval<Real>
solver<Ring,CFfirstIsolate >::solve(const POL& p,
const base_t & u,
const base_t & v) {

  interval_rep<POL>  i(p,u,v);
  solver_cffirst<Real,POL> slv(i);

  return slv.first_root_isolate();
}

template<class Ring> 
template<class Real, class POL> 
Interval<Real>
solver<Ring,CFfirstIsolate >::solve(const POL& p) {

  interval_rep<POL> i(p);
  solver_cffirst<Real,POL> slv(i);

  return slv.first_root_isolate();
}


//--------------------------------------------------------------------
    template<class Ring>
    struct solver<Ring, CFfirstApproximate > {

        typedef typename Ring::scalar_t   base_t;

        template<class Real, class POL> static Real
        solve(const POL& p, const base_t& u, const base_t& v);
        
        template<class Real, class POL> static Real
        solve(const POL& p);
};

    template<class Ring> 
    template<class Real, class POL>
    Real

    solver<Ring,CFfirstApproximate >::solve(const POL& p,
                                 const base_t& u, const base_t& v) {

  
        interval_rep<POL>  i(p,u,v);
        solver_cffirst<Real,POL> slv(i);

        return slv.first_root_approximate();
}

    template<class Ring> 
    template<class Real, class POL> 
    Real
    solver<Ring,CFfirstApproximate >::solve(const POL& p) {


        interval_rep<POL> i(p);
        solver_cffirst<Real,POL> slv(i);

  return slv.first_root_approximate();
}


//--------------------------------------------------------------------
    template<class Ring>
    struct solver<Ring, CFfirstFloor> {
        typedef typename Ring::Scalar  base_t;

        template<class Real, class POL> static Real
        solve(const POL& p, const base_t& u, const base_t& v);

        template<class Real,class POL> static Real
        solve(const POL& p);
};

    template<class Ring> 
    template<class Real,class POL> 
    Real
    solver<Ring,CFfirstFloor >::solve(const POL& p,
                                 const base_t& u, const base_t& v) {

        interval_rep<POL>  i(p,u,v);
        solver_cffirst<Real,POL> slv(i);
        
        return slv.first_root_floor();
}

    template<class Ring> 
    template<class Real, class POL> 
    Real
    solver<Ring,CFfirstFloor >::solve(const POL& p) {
        
        interval_rep<POL> i(p);
        solver_cffirst<Real,POL> slv(i);
        
        return slv.first_root_floor();
}


//--------------------------------------------------------------------
    template<class Ring>
    struct solver<Ring, CFallIsolate > {

        typedef typename Ring::Scalar   base_t;

        template<class Real, class POL> static Seq<Interval<Real> >
        solve(const POL& p,const base_t& u,const base_t& v);

        template<class Real,class POL> static Seq<Interval<Real> >
        solve(const POL& p);
};

    template<class Ring> 
    template<class Real, class POL> 
    Seq<Interval<Real> >
    solver<Ring,CFallIsolate >::solve(const POL& p,
                                 const base_t& u, const base_t& v) {

        interval_rep<POL>  i(p,u,v);
        solver_cffirst<Real,POL> slv(i);

  return slv.all_roots_isolate();
}

    template<class Ring> 
    template<class Real,class POL> 
    Seq<Interval<Real> >
    solver<Ring,CFallIsolate >::solve(const POL& p) {

        interval_rep<POL> i(p);
        solver_cffirst<Real,POL> slv(i);

        return slv.all_roots_isolate();
}


//--------------------------------------------------------------------
    template<class Ring>
    struct solver<Ring, CFseparate > {

        typedef typename Ring::Scalar  base_t;

        template<class Real,class POL> static Seq<Real>
        solve(const POL& p,const base_t& u,const base_t& v);

        template<class Real, class POL> static Seq<Real>
        solve(const POL& p);
};

    template<class Ring> 
    template<class Real, class POL> 
    Seq<Real>
    solver<Ring,CFseparate >::solve(const POL& p,
                                 const base_t& u, const base_t& v) {

        interval_rep<POL>  i(p,u,v);
        solver_cffirst<Real,POL> slv(i);
        
        return slv.all_roots_separate();
}

    template<class Ring> 
    template<class Real,class POL> 
    Seq<Real>
    solver<Ring,CFseparate >::solve(const POL& p) {

        interval_rep<POL> i(p);
        solver_cffirst<Real,POL> slv(i);
        
        return slv.all_roots_separate();
}


//--------------------------------------------------------------------
    template<class Ring> 
    struct solver<Ring,CFdecide> {

        template<class POL> static bool
        solve(const POL& p);
};

    template<class Ring> 
    template<class POL> 
    bool
    solver<Ring,CFdecide>::solve(const POL& p) {

        interval_rep<POL> i(p);
        solver_cffirst<QQ,POL> slv(i);
        
        return ( slv.all_roots_isolate().size() );
}


} //namespace realroot

#undef TMPL
#undef TMPLR
#undef NO_ROOT
#undef N_ITER
#endif // _realroot_SOLVER_FIRST_ROOT_CONFRAC_H