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\c 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/GMP.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)
    {
          //std::cout<< "Cauchy" <<std::endl;
          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));
          //std::cout<< "Max "<< M.max <<" " <<an<< std::endl;
          //      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;
  }

  //====================================================================
  template <typename FT, typename POLY>
  FT bound_root(const POLY& p, const Cauchy<FT>& m)
  {
    return Cauchy<FT>::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