|
realroot_doc 0.1.1
|
|
realroot_doc 0.1.1
|
#include <contfrac_lowerbound.hpp>
A lower bound using the Bernstein representation of the polynomial. The idea is to compute the Bernstein form of the polynomial over the positive axis and find the least positive intersection point of the convex hull with the axis. From the convex hull property it follows that this quantity is a lower bound on the positive roots of the polynomial. The value returned is the floor of the intersection point.
Definition at line 19 of file contfrac_lowerbound.hpp.
| Poly::value_type lower_bound | ( | const Poly & | f | ) | [inline] |
Definition at line 22 of file contfrac_lowerbound.hpp.
References mmx::assign(), mmx::bit_size(), mmx::degree(), mmx::pow(), and mmx::sign().
{
//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);
}
}
