#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().
00023 { 00024 //std::cout <<"Bezier"<<std::endl; 00025 typedef typename Poly::value_type RT; 00026 unsigned deg = degree(f); 00027 tensor::bernstein<RT> bz(deg+1); 00028 let::assign(bz, f); 00029 // The Bernstein form of f w.r.t. the interval [0,\infinity] 00030 00031 int sn = sign(bz[0]); 00032 if(sn == 0) return 0;// That is, zero is a root, this should not be the case 00033 00034 //std::cout << "bz[0] = "<< bz[0] << std::endl; 00035 00036 RT lowerBoundN =-1, lowerBoundD=1;// Numerator and denominator of lower bound 00037 //Now we find the least positive intersection of the convex hull. 00038 for(int i=1; i <= deg; i++){ 00039 if(sn * sign(bz[i]) < 0){ 00040 /* 00041 std::cout << "bz["<<i<<"] = "<< bz[i] << std::endl; 00042 std::cout << "bz[0]*i "<< bz[0] * i << std::endl; 00043 std::cout << "(deg*(bz[0]-bz[i]))" << (deg*(bz[0]-bz[i])) << std::endl; 00044 */ 00045 if(lowerBoundN * deg *(bz[0]-bz[i]) > bz[0]*i*lowerBoundD 00046 || lowerBoundN == -1){ 00047 lowerBoundN = bz[0]*i; 00048 lowerBoundD = deg*(bz[0]-bz[i]); 00049 } 00050 } 00051 } 00052 lowerBoundD -= lowerBoundN;// Have to scale back the result to the orignial setting of [0, B] 00053 long K = bit_size(lowerBoundN)- bit_size(lowerBoundD) -1; 00054 if(K < 0) 00055 return RT(0); 00056 else if(lowerBoundN == lowerBoundD) 00057 return 1; 00058 else{// K > 0 00059 return pow( RT(2), K); 00060 } 00061 }