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