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

Go to the documentation of this file.

/*******************************************************************
 *   This file is part of the source code of the realroot kernel.
 *   Author(s): Elias P. TSIGARIDAS <et@di.uoa.gr>
 *              J.P. Pavone, GALAAD, INRIA
 *              V. Sharma, GALAAD, INRIA
 ********************************************************************/
#ifndef realroot_UPOL_BOUND_H
#define realroot_UPOL_BOUND_H

#include <realroot/rounding_mode.hpp>
#include <realroot/texp_rationalof.hpp>
#include <realroot/univariate.hpp>
#include <vector>
#include <stack>
#include <queue>


namespace mmx {

using univariate::lcoeff;

#ifndef MMX_ABS_FCT
#define MMX_ABS_FCT
  template<typename T> inline T abs(const T& x) { return (x < 0 ? -x : x); }
#endif


  template<class T>
  T pow2(int i)
  {
    assert(i>=0);
    if(i == 1) return 2;
    if(i > 0)
    {
      T y(2);
      T z(1);
      unsigned int n = i;
      unsigned int m;
      while (n > 0) {
        m = n / 2;
        if (n %2 )
          {
            z *= y;
          }
        y = y * y;
        n = m;
      }
      return z;
    }
    return T(1);
  }


  template<class C> struct NISP
  {
    //--------------------------------------------------------------------
    template <typename POLY> static
    C upper_bound(const POLY& p)
    {
      typedef typename texp::rationalof<C>::T rational;
      //  typedef typename NTR::T                FT;
      //typedef double FT;
      C max(0), sum(0);
    
      unsigned i;
      unsigned n = p.size();
    
      // check the degree 
      i = n; while ( p[--i] == 0 ) ; n = i+1;
      if ( p[i] > 0 )
        {
          i = 0;
          while ( p[i] > 0 && i < n ) i++;
          if ( i == n ) return 0;
        
          for ( unsigned k = i+1; k < n; k ++ )
            if ( p[k] > 0 ) sum += p[k];
        
          for ( unsigned k = i; k < n; k ++ )
            {
              if ( p[k] < 0 ) 
                { rational tmp = -p[k]/sum; if ( tmp > max ) max = tmp; }
              else 
                sum -= p[k];
            };
        }
      else
        {
          i = 0;
          while ( p[i] < 0 && i < n ) i++;
          if ( i == n ) return 0;
          for ( unsigned k = i+1; k < n; k ++ )
            if ( p[k] < 0 ) sum += p[k];
          for ( unsigned k = i; k < n; k ++ )
            if ( p[k] > 0 ) { rational tmp = -p[k]/sum; if ( tmp > max ) max = tmp; }
            else sum -= p[k];
        };
      max += 1;
      return max;
    }


    template <typename POLY> static
    C lower_bound(const POLY& p)
    {
      typedef typename texp::rationalof<C>::T rational;
      using univariate::degree;

      POLY r(p);
      reverse(r,degree(r));
      rational m = NISP<rational>::upper_bound(r);
      return denominator(m)/numerator(m);
    }
  };

  template<class FT> struct NISN
  {
    template <typename POLY> static
    FT upper_bound(const POLY& p)
    {
      POLY tmp = p;
      for ( unsigned i = 1; i < p.size(); i += 2 ) 
        tmp[i] = -tmp[i];
      return  -NISP<FT>::upper_bound(tmp);
    }
  };

  //====================================================================
  template < typename T >
  struct abs_max : public std::unary_function<T, void>
  {
    abs_max() : max(T(-1)) {}
    void operator() (const T& x)
    {
        // using std::abs;
      T temp = abs(x);
      if (temp > max) max = temp;
    }
    T max;
  };

  template<class C> struct Cauchy 
  {
    //--------------------------------------------------------------------
    template<class Poly> static
    C upper_bound(const Poly& p)
    {
      //  DEBUG::print("Cauchy");
      typedef typename Poly::value_type     RT;
      typedef typename texp::rationalof<RT>::T rational;

      //using std::abs;
      abs_max<RT> M = std::for_each(p.begin(), p.end()-1, abs_max<RT>());
      RT an =  abs(lcoeff(p));
      //      rational val= rational(M.max)/an + rational(1); 
      C res= C(M.max/an) + C(1); 
      //      C res; let::assign(res, val); 
      return res;
    }

    template<class Poly>
    C lower_bound(const Poly& f) const
    {
      using univariate::degree;

      typedef C  RT;
      long deg = degree( f);
      long l = 0;
    
      for (int i = 1; i <= deg; ++i) {
        if ( sign( f[i]) * sign( f[0]) < 0 ) {
          ++l;
        }
      }
    
      long ba0  = bit_size(f[0]) - 1; /* puiss de 2 < an */
      bool bound_set = false;
      long K = 0;
    
      for (int i = 1; i <= deg; ++i) {
        if ( sign( f[i]) * sign( f[0]) < 0 ) {
          long bai  = bit_size( RT( l * f[i])) - 1;
          long p = bai - ba0 - 1 ;
          long k = i;
          long q = p / k;
          long r = p % k;
        
          long rq = 0;
          if (r <= k-2) {
            rq = q + 1;
          } else {
            rq = q  + 2;
          }
        
          if ( (bound_set == false) || (K < rq) ) {
            K = rq;
            bound_set = true;
          }
        }
      }
      if ( K < 0 ) {return pow2<RT>( abs( K) );}
    
      return 0;
    }

    // Return only the power of 2.
    template<class Poly>
    long lower_power_2(const Poly& f) const
    {
      using univariate::degree;
 
      typedef typename Poly::value_type RT;
      
      long deg = degree( f);
      long l = 0;
    
      for (int i = 1; i <= deg; ++i) {
        if ( sign( f[i]) * sign( f[0]) < 0 ) {
          ++l;
        }
      }
    
      long ba0  = bit_size(f[0]) - 1; /* puiss de 2 < an */
      bool bound_set = false;
      long K = 0;
    
      for (int i = 1; i <= deg; ++i) {
        if ( sign( f[i]) * sign( f[0]) < 0 ) {
          long bai  = bit_size( RT( l * f[i])) - 1;
          long p = bai - ba0 - 1 ;
          long k = i;
          long q = p / k;
          long r = p % k;
        
          long rq = 0;
          if (r <= k-2) {
            rq = q + 1;
          } else {
            rq = q  + 2;
          }
        
          if ( (bound_set == false) || (K < rq) ) {
            K = rq;
            bound_set = true;
          }
        }
      }
      if ( K < 0 ) { abs( K);}

      return -1;
    }

  };

  template<>
  template<class Poly>
  double Cauchy<double>::upper_bound(const Poly& p)
  {
    typedef typename Poly::value_type coeff_t;
    numerics::rounding<double> rnd( numerics::rnd_up() );
    //using std::abs;
    coeff_t M=0;
    for(unsigned i=0;i< p.size(); i++) M = std::max(M, coeff_t(abs(p[i])));
    //    abs_max<coeff_t> M = std::for_each(p.begin(), p.end()-1, abs_max<coeff_t>());
    
    coeff_t an = coeff_t(abs(univariate::lcoeff(p)));
    //std::cout <<M<<" "<<an<<" "<<std::endl;

    double res=as<double>(coeff_t(M/an)) + 1; 
    //    std::cout <<res<<std::endl;
    return res;
  }
  //====================================================================

  // /** 
  //  * @brief Prototype function, Computes an upper bound for the real roots.
  //  * 
  //  * @param p Univariate Polynomial
  //  * @param mth Method that we use in order to compute the root bound
  //  * @param AsType<T> The return type
  //  */
  // template < typename POLY,
  //       typename M,
  //       typename T >
  // T bound(const POLY& p, const M& mth, texp::As<T>)
  // {
  //   return T(0);
  // }



  template <typename FT, typename POLY>
  FT bound_root(const POLY& p, const Cauchy<FT>& m)
  {
    return m.upper_bound(p);
  }


  template < typename FT, typename POLY>
  FT max_bound(const POLY& p, const Cauchy<FT>& m)
  {
    return bound_root(p,m);
  }

  template < typename FT, typename POLY>
  FT
  min_bound(const POLY& p, Cauchy<FT>)
  {
    using univariate::degree;

    POLY f(p);
    univariate::reverse(f, degree(p));
    FT tmp =  bound_root(f, Cauchy<FT>());
    //    let::convert(denominator(tmp), numerator(tmp), texp::As<FT>());
    return FT(1)/tmp;
  }




  //====================================================================
  template<class RT> struct HongBound
  {
    template < typename Poly > static
    RT lower_bound( const Poly& f)
    {
      using univariate::degree;

      unsigned deg = degree(f);
    
      long lB=0, gB=0;
      bool localBoundSet = false, globalBoundSet = false;
      long temp, q;
      int s = sign(f[0]);
      //std::cout << "f[0]= " << f[0] << std::endl;
      for(int i= deg; i > 0; i--){
        //std::cout << "f["<< i << "]" << " = " << f[i] << std::endl;
        if(sign(f[i]) * s < 0){
          for(int k=i-1; k >= 0; k--){
            if(sign(f[k]) * s > 0){
              temp = bit_size( RT(f[i]) ) - bit_size( RT(f[k]) ) - 1;
              q = temp /(i-k);
              if(!localBoundSet || lB > q + 2){// Choose the minimum among localBounds
                localBoundSet = true;
                lB = q+2;
              }
            }
          }
          localBoundSet = false;
          if(!globalBoundSet || gB < lB){// Choose the maximum among globalBounds
            globalBoundSet = true;
            gB = lB;
          }
        }
        //std::cout << "Cout gb after " << i << " loop = "<< gB << std::endl;
      }
      if ( gB+1 < 0 ) return pow2<RT>( abs( gB+1) );
      return 0;
    }
  };

  //====================================================================
  template<class RT> 
  struct AkritasBound
  {

    template < typename Poly > static
    RT lower_bound( const Poly& f)
    {
      using univariate::degree;

      //    std::cout <<"Poly is " << f << std::endl;
      typedef std::pair<RT, int>        Coeff_deg;
      Coeff_deg T1, T2;

      unsigned deg=degree(f);
      int s=sign(f[0]);//, len;
      long B=0;// The logarithm of the bound to be returned
      bool boundSet=false;// The bound is not set initially
      long temp,q;
      int posCounter=0, negCounter=0, nPosCounter=0;
      // The coefficients from posCounter to negCounter-1 have the same sign
      // as f[0]; the coefficients from negCounter to nPosCounter-1 have the
      // opposite sign of f[0].
      std::queue< Coeff_deg > Neg;// Queue that stores the negative (between
      // negCounter and nPosCounter) coefficients.
      std::stack< Coeff_deg > Pos; // Stack that stores the positive (between
      // posCounter and negCounter) coefficients.

      int l1, l2;// size of the above two queues

      while(posCounter != int(deg+1)){
        // Find the first change in sign larger than posCounter and assign negCounter
        // to this value. If no such sign variation occurs then
        // set negCounter to -1.
        negCounter = -1;// By default there is no such coefficient
        //Pos.push(Coeff_deg(RT(f[posCounter]), posCounter));// Push the first "positive" coefficient
        Pos.push(Coeff_deg(RT(f[posCounter]), posCounter));// Push the first "positive" coefficient
        for(unsigned i=posCounter + 1;i <= deg; i++){
          if(sign(f[i]) * s < 0){
            negCounter=i; break;
          }
          if(sign(f[i]) *s > 0){// Store the list of "positive" coefficients
            //Pos.push(Coeff_deg(RT(f[i]), i));// Note the actual degree is deg-i, but
            // later we subtract degrees so we store only i
            Pos.push(Coeff_deg(RT(f[i]), i));
          }
        }
        if(negCounter == -1){
          //    std::cout <<"Returning the lower bound " << B << std::endl;
          if ( B < 0 ) 
            return pow2<RT>( abs(B) );
          else
            return 0;
        }
        // Now find the next coefficient that has the same sign as f[0].
        nPosCounter = deg+1;// By default this is deg+1
        Neg.push(Coeff_deg(RT(f[negCounter]), negCounter));
        for(unsigned i=negCounter + 1;i <= deg; i++){
          if(sign(f[i]) * s > 0){
            nPosCounter=i; break;
          }
          if(sign(f[i]) *s < 0)// Store the list of "negative" coefficients
            Neg.push(Coeff_deg(RT(f[i]), i));
        }

        //      std::cout <<" Assigning to the queues done" << std::endl;
        // Start computing the bound
        // If the number of "positive" terms from posCounter to negCounter -1 are
        // greater or equal to  the number of "negative" terms from negCounter to
        // nPosCounter -1 then we have the straightforward matching of "negative"
        // coefficients with the "positive" ones.
        l1 = Neg.size(); l2 = Pos.size();
        if(l2 >= l1){
          //    std::cout <<"Pos length " << l2 << " greater than Neg length " 
          //              << l1 << std::endl;
          // or equally negCounter - posCounter >= nPosCounter - negCounter
          for(int i=0; i < l1; i++){
            T1= Neg.front(); Neg.pop();
            //T2= Pos.front(); Pos.pop();
            T2=Pos.top();Pos.pop();
            //    std::cout <<"Neg coeff " << T1.first << " Pos coeff " << T2.first
            //              << " Deg difference " << T1.second - T2.second 
            //              << std::endl;
            temp = bit_size( T1.first ) - bit_size( T2.first ) - 1;
            q = temp/(T1.second - T2.second);
            //    std::cout << "Difference between logs " << temp 
            //              << " temporary bound value " << B << std::endl;
            if(!boundSet || B < q + 2){// Choose the maximum
              boundSet = true;
              B = q+2;
            }
          }
        }else{//Else split the coefficient f[posCounter] to compensate for
          // the paucity of "positive" terms and then do the matching.
          //    std::cout <<" Pos length "<< l2 << " smaller than Neg length " 
          //              << l1 << std::endl;
        
          //T2=Pos.front(); Pos.pop();
          T2= Pos.top();Pos.pop();
          for(int i=0; i <= l1-l2; i++){
            T1= Neg.front(); Neg.pop();
            temp = bit_size( T1.first ) - bit_size( T2.first ) + bit_size(RT(l1-l2+1)) - 1;
            q = temp/(T1.second-T2.second);
            if(!boundSet || B < q + 2){// Choose the maximum
              boundSet = true;
              B = q+2;
            }
          }
          //    std::cout <<" Splitting the leading coefficient done" << std::endl;
          l1 = Neg.size();// Now the size of Neg and Pos should be the same
          for(int i=0; i < l1; i++){
            T1= Neg.front(); Neg.pop();
            //T2= Pos.front(); Pos.pop();
            Pos.pop();
            temp = bit_size( T1.first ) - bit_size( T2.first ) - 1;
            q = temp/(T1.second - T2.second);
            if(!boundSet || B < q + 2){// Choose the maximum
              boundSet = true;
              B = q+2;
            }
          }
        }// end else
        //      std::cout <<" Matching the coefficients done" << std::endl;
        posCounter = nPosCounter; Neg.empty(); Pos.empty();
      }//end while
      //    std::cout <<"Returning the lower bound " << B << std::endl;
      if ( B < 0 ) return pow2<RT>( abs(B) );
      return 0;
    }


      template < typename Poly > static
      long lower_power_2( const Poly& f)
      {
        using univariate::degree;

          //    std::cout <<"Poly is " << f << std::endl;
          typedef std::pair<RT, int>    Coeff_deg;
          Coeff_deg T1, T2;

          unsigned deg=degree(f);
          int s=sign(f[0]);//, len;
          long B=0;// The logarithm of the bound to be returned
          bool boundSet=false;// The bound is not set initially
          long temp,q;
          int posCounter=0, negCounter=0, nPosCounter=0;
          // The coefficients from posCounter to negCounter-1 have the same sign
          // as f[0]; the coefficients from negCounter to nPosCounter-1 have the
          // opposite sign of f[0].
          std::queue< Coeff_deg > Neg;// Queue that stores the negative (between
          // negCounter and nPosCounter) coefficients.
          std::stack< Coeff_deg > Pos; // Stack that stores the positive (between
          // posCounter and negCounter) coefficients.

          int l1, l2;// size of the above two queues

          while(posCounter != int(deg+1)){
              // Find the first change in sign larger than posCounter and assign negCounter
              // to this value. If no such sign variation occurs then
              // set negCounter to -1.
              negCounter = -1;// By default there is no such coefficient
              //Pos.push(Coeff_deg(RT(f[posCounter]), posCounter));// Push the first "positive" coefficient
              Pos.push(Coeff_deg(RT(f[posCounter]), posCounter));// Push the first "positive" coefficient
              for(unsigned i=posCounter + 1;i <= deg; i++){
                  if(sign(f[i]) * s < 0){
                      negCounter=i; break;
                  }
                  if(sign(f[i]) *s > 0){// Store the list of "positive" coefficients
                      //Pos.push(Coeff_deg(RT(f[i]), i));// Note the actual degree is deg-i, but
                      // later we subtract degrees so we store only i
                      Pos.push(Coeff_deg(RT(f[i]), i));
                  }
              }
              if(negCounter == -1){
                  //    std::cout <<"Returning the lower bound " << B << std::endl;
                  if ( B < 0 ) return abs( B);
                  else return -1;
              }
              // Now find the next coefficient that has the same sign as f[0].
              nPosCounter = deg+1;// By default this is deg+1
              Neg.push(Coeff_deg(RT(f[negCounter]), negCounter));
              for(unsigned i=negCounter + 1;i <= deg; i++){
                  if(sign(f[i]) * s > 0){
                      nPosCounter=i; break;
                  }
                  if(sign(f[i]) *s < 0)// Store the list of "negative" coefficients
                      Neg.push(Coeff_deg(RT(f[i]), i));
              }

              //      std::cout <<" Assigning to the queues done" << std::endl;
              // Start computing the bound
              // If the number of "positive" terms from posCounter to negCounter -1 are
              // greater or equal to  the number of "negative" terms from negCounter to
              // nPosCounter -1 then we have the straightforward matching of "negative"
              // coefficients with the "positive" ones.
              l1 = Neg.size(); l2 = Pos.size();
              if(l2 >= l1){
                  //    std::cout <<"Pos length " << l2 << " greater than Neg length " 
                  //              << l1 << std::endl;
                  // or equally negCounter - posCounter >= nPosCounter - negCounter
                  for(int i=0; i < l1; i++){
                      T1= Neg.front(); Neg.pop();
                      //T2= Pos.front(); Pos.pop();
                      T2=Pos.top();Pos.pop();
                      //          std::cout <<"Neg coeff " << T1.first << " Pos coeff " << T2.first
                      //                    << " Deg difference " << T1.second - T2.second 
                      //                    << std::endl;
                      temp = bit_size( T1.first ) - bit_size( T2.first ) - 1;
                      q = temp/(T1.second - T2.second);
                      //          std::cout << "Difference between logs " << temp 
                      //                    << " temporary bound value " << B << std::endl;
                      if(!boundSet || B < q + 2){// Choose the maximum
                          boundSet = true;
                          B = q+2;
                      }
                  }
              }else{//Else split the coefficient f[posCounter] to compensate for
                  // the paucity of "positive" terms and then do the matching.
                  //    std::cout <<" Pos length "<< l2 << " smaller than Neg length " 
                  //              << l1 << std::endl;
        
                  //T2=Pos.front(); Pos.pop();
                  T2= Pos.top();Pos.pop();
                  for(int i=0; i <= l1-l2; i++){
                      T1= Neg.front(); Neg.pop();
                      temp = bit_size( T1.first ) - bit_size( T2.first ) + bit_size(RT(l1-l2+1)) - 1;
                      q = temp/(T1.second-T2.second);
                      if(!boundSet || B < q + 2){// Choose the maximum
                          boundSet = true;
                          B = q+2;
                      }
                  }
                  //    std::cout <<" Splitting the leading coefficient done" << std::endl;
                  l1 = Neg.size();// Now the size of Neg and Pos should be the same
                  for(int i=0; i < l1; i++){
                      T1= Neg.front(); Neg.pop();
                      //T2= Pos.front(); Pos.pop();
                      Pos.pop();
                      temp = bit_size( T1.first ) - bit_size( T2.first ) - 1;
                      q = temp/(T1.second - T2.second);
                      if(!boundSet || B < q + 2){// Choose the maximum
                          boundSet = true;
                          B = q+2;
                      }
                  }
              }// end else
              //      std::cout <<" Matching the coefficients done" << std::endl;
              posCounter = nPosCounter; Neg.empty(); Pos.empty();
          }//end while
          //    std::cout <<"Returning the lower bound " << B << std::endl;
          if ( B < 0 ) return abs( B);
          else return -1;

      }
  };
    //====================================================================

} //namespace mmx
//--------------------------------------------------------------------
#endif // realroot_UPOL_BOUND_H