contfrac_lowerbound.hpp

Go to the documentation of this file.

#ifndef _synaps_usolve_contfrac_lower_bound_h_
#define _synaps_usolve_contfrac_lower_bound_h_

#include <realroot/univariate_bounds.hpp>
#include <realroot/tensor_bernstein.hpp>

namespace mmx {
  // #define bit_size(a) mpz_sizeinbase(a.get_mpz_t(),2)

//--------------------------------------------------------------------
struct BezierBound
{
  template < typename Poly >
  typename Poly::value_type lower_bound( const Poly& f)
  {
    //std::cout <<"Bezier"<<std::endl;
    typedef typename Poly::value_type    RT;
    unsigned deg = degree(f);
    tensor::bernstein<RT> bz(deg+1); 
    let::assign(bz, f);  
    // The Bernstein form of f w.r.t. the interval [0,\infinity]
    
    int sn = sign(bz[0]);
    if(sn == 0) return 0;// That is, zero is a root, this should not be the case
    
    //std::cout << "bz[0] = "<< bz[0] << std::endl;
    
     RT lowerBoundN =-1, lowerBoundD=1;// Numerator and denominator of lower bound
    //Now we find the least positive intersection of the convex hull.
    for(int i=1; i <= deg; i++){
      if(sn * sign(bz[i]) < 0){
    /*
      std::cout << "bz["<<i<<"] = "<< bz[i] << std::endl;
      std::cout << "bz[0]*i "<< bz[0] * i << std::endl;
      std::cout << "(deg*(bz[0]-bz[i]))" << (deg*(bz[0]-bz[i])) << std::endl;
    */
    if(lowerBoundN * deg *(bz[0]-bz[i]) > bz[0]*i*lowerBoundD 
       || lowerBoundN == -1){
      lowerBoundN = bz[0]*i;
      lowerBoundD = deg*(bz[0]-bz[i]);
    }
      }
    }
    lowerBoundD -= lowerBoundN;// Have to scale back the result to the orignial setting of [0, B]
    long K = bit_size(lowerBoundN)- bit_size(lowerBoundD) -1;
    if(K < 0)
      return RT(0);
    else if(lowerBoundN == lowerBoundD)
      return 1;
    else{// K > 0
      return pow( RT(2), K);    
    }
  }
};


  
} //namespace mmx
#endif //_synaps_usolve_contfrac_lower_bound_h_