|
realroot_doc 0.1.1
|
|
realroot_doc 0.1.1
|
00001 /******************************************************************* 00002 * This file is part of the source code of the realroot kernel. 00003 * Description: The Interval Newton method to approximate a root of a 00004 * polynomial given an initial interval that contains it. 00005 * Author(s): 00006 * G. Tzoumas <geotz@di.uoa.gr> 00007 * Department of Informatics and Telecommunications 00008 * University of Athens (Greece). 00009 ********************************************************************/ 00010 #ifndef _realroot_solve_inewton_hpp_ 00011 #define _realroot_solve_inewton_hpp_ 00012 00020 //#include <numerix/interval.hpp> 00021 //#include <realroot/Interval.hpp> 00022 #include <realroot/polynomial_tensor.hpp> 00023 #include <realroot/Seq.hpp> 00024 00025 namespace mmx { 00026 00027 /* // perhaps it might come handy in the future -- I leave it in comments 00028 namespace inewton { 00029 00030 template<class C> struct IntervalTraits { }; 00031 00032 template<> template<class B> struct IntervalTraits<Interval<B> > { 00033 typedef Interval<B> interval_type; 00034 typedef B boundary_type; 00035 00036 static inline boundary_type upper(const interval_type& x) { return x.upper(); } 00037 static inline boundary_type lower(const interval_type& x) { return x.lower(); } 00038 static inline boundary_type width(const interval_type& x) { return x.width(); } 00039 static inline boundary_type median(const interval_type& x) { return x; } 00040 static inline int sign(const interval_type& x) { return ::mmx::sign(x); } 00041 static inline bool intersect(interval_type &r, const interval_type &a, const interval_type &b) { 00042 return ::mmx::intersect(r, a, b); 00043 } 00044 }; 00045 00046 template<> template<class B> struct IntervalTraits<interval<B> > { 00047 typedef interval<B> interval_type; 00048 typedef B boundary_type; 00049 00050 static inline boundary_type upper(const interval_type& x) { return ::mmx::upper(x); } 00051 static inline boundary_type lower(const interval_type& x) { return ::mmx::lower(x); } 00052 static inline boundary_type width(const interval_type& x) { return upper(x)-lower(x); } 00053 static inline boundary_type median(const interval_type& x) { return ::mmx::center(x); } 00054 static inline int sign(const interval_type& x) { return ::mmx::sign(x); } 00055 static inline bool intersect(interval_type &r, const interval_type &a, const interval_type &b) { 00056 if ((upper(a) < lower(b)) || (lower(a) > upper(b))) return false; 00057 r = interval_type(std::max(lower(a), lower(b)), std::min(upper(a), upper(b))); 00058 return true; 00059 } 00060 }; 00061 00062 } */ 00063 00064 template<class INT, class CT> struct interval_newton 00065 { 00066 typedef typename ring<CT, MonomialTensor>::Polynomial polynomial; 00067 typedef typename INT::boundary_type BT; 00068 00069 static inline BT width(const INT& x) { return mmx::upper(x) - mmx::lower(x); } 00070 static inline BT median(const INT& x) { return mmx::lower(x) + width(x)/BT(2); } 00071 static bool singleton(const INT& x) { return mmx::lower(x) == mmx::upper(x); } 00072 00073 static INT bisection_iteration(const polynomial& p, const INT approx, int& status_) 00074 { 00075 BT x = median(approx); 00076 INT lval, rval, mval; 00077 lval = p(INT(mmx::lower(approx))); 00078 rval = p(INT(mmx::upper(approx))); 00079 mval = p(INT(x)); 00080 // std::cerr << lval << ' ' << rval << ' ' << mval << ' ' << x << std::endl; 00081 // assert( mmx::sign(lval)*mmx::sign(rval) <= 0 ); 00082 00083 status_ = 0; 00084 if (mval == 0) return INT(x); 00085 if (lval == 0) return INT(mmx::lower(approx)); 00086 if (rval == 0) return INT(mmx::upper(approx)); 00087 00088 if (mmx::sign(lval) == 0 || mmx::sign(rval) == 0) { 00089 status_ = -3; return approx; 00090 } 00091 00092 if (mmx::sign(lval) == mmx::sign(rval)) { 00093 status_ = -1; return approx; 00094 } 00095 00096 if (mmx::sign(mval) == 0) return approx; 00097 00098 if (mmx::sign(lval) == mmx::sign(mval)) return INT(x,mmx::upper(approx)); 00099 return INT(mmx::lower(approx), x); 00100 } 00101 00102 /* 00103 status: 00104 -4 = no convergence (should not happen!) 00105 -3 = precision not enough 00106 -2 = no root exists 00107 -1 = no root found (cannot guarantee existence) 00108 0 = root found (cannot guarantee uniqueness) 00109 1 = unique root found 00110 */ 00111 // caution with coeffs including zero! 00112 static INT newton_iteration(const polynomial& p, const polynomial& dp, 00113 INT approx, BT delta_, int& status_) 00114 { 00115 INT slope, testapprox, newapprox, ver; 00116 BT oldw = width(approx); 00117 BT neww; 00118 int z; 00119 00120 #ifdef INEWTONDEBUG 00121 int i = 0; 00122 #endif 00123 00124 INT lval, rval; 00125 00126 lval = p(INT(mmx::lower(approx))); 00127 rval = p(INT(mmx::upper(approx))); 00128 #ifdef INEWTONDEBUG 00129 std::cerr << "approx = " << approx << " lval = " << lval << " rval = " << rval << std::endl; 00130 #endif 00131 if (lval == 0) { status_ = 0; return INT(mmx::lower(approx)); } 00132 if (rval == 0) { status_ = 0; return INT(mmx::upper(approx)); } 00133 if (mmx::sign(lval) == 0 || mmx::sign(rval) == 0) { 00134 status_ = -3; return approx; 00135 } 00136 00137 status_ = 0; 00138 slope = dp(approx); 00139 00140 #ifdef INEWTONDEBUG 00141 std::cerr << "slope = " << slope << width(slope) << std::endl; 00142 #endif 00143 00144 z = mmx::sign(slope); 00145 if (mmx::sign(lval) == mmx::sign(rval)) { 00146 if (z == 0) status_ = -1; else status_ = -2; 00147 return approx; 00148 } 00149 00150 #ifdef INEWTONDEBUG 00151 std::cerr << "solving " << p << std::endl; 00152 #endif 00153 00154 /* INTERVAL NEWTON MAIN LOOP */ 00155 while (1) { 00156 // should we always do simple bisection with intervals of length >= 1 ? 00157 // if (oldw >= 1 || z == 0) { 00158 00159 // slope contains zero 00160 if (z == 0) { 00161 approx = bisection_iteration(p, approx, status_); 00162 neww = width(approx); 00163 #ifdef INEWTONDEBUG 00164 i++; 00165 std::cerr << '(' << i << ") bisect: " << approx << ' ' << neww << std::endl; 00166 #endif 00167 if (status_ < 0) return approx; 00168 if (neww < delta_) break; 00169 assert(neww <= oldw); 00170 if (neww == oldw) break; 00171 oldw = neww; 00172 00173 } else { 00174 // slope doesn't contain zero 00175 00176 BT x = median(approx); 00177 if (x == mmx::lower(approx) || x == mmx::upper(approx)) break; 00178 #ifdef INEWTONDEBUG 00179 std::cerr << "mid = " << x << " mid_INT = " << INT(x) << std::endl; 00180 std::cerr << "p(x) = " << p(INT(x)) << std::endl; 00181 std::cerr << "p'(approx) = " << dp(approx) << std::endl; 00182 std::cerr << "p/p' = " << INT(p(INT(x))/dp(approx)) << std::endl; 00183 #endif 00184 testapprox = INT(x) - INT(p(INT(x))/dp(approx)); // Newton approximation 00185 00186 #ifdef INEWTONDEBUG 00187 i++; 00188 #endif 00189 if (mmx::intersect(newapprox, testapprox, approx)) { 00190 #ifdef INEWTONDEBUG 00191 std::cerr <<"new =" << newapprox << " test = " << testapprox << " old = " << approx << std::endl; 00192 #endif 00193 neww = width(newapprox); 00194 00195 // new approximation SHOULD contain the root 00196 // if you want to verify, use the debugging code below 00197 #ifdef INEWTONDEBUG 00198 INT ver = p(newapprox); 00199 if (ver != INT(0)) { 00200 assert (mmx::sign(ver) == 0); 00201 } 00202 00203 std::cerr << '(' << i << ") Newton: " << newapprox << ' ' << neww << std::endl; 00204 #endif 00205 //if (Float(2)*neww >= oldw) break; 00206 if (neww >= oldw) break; 00207 approx = newapprox; 00208 if (neww < delta_) break; 00209 oldw = neww; 00210 } else { 00211 // no convergence 00212 status_ = -4; 00213 return approx; 00214 } 00215 } 00216 slope = dp(approx); 00217 #ifdef INEWTONDEBUG 00218 std::cerr << "slope = " << slope << width(slope) << std::endl; 00219 #endif 00220 z = mmx::sign(slope); 00221 } 00222 if (mmx::sign(dp(approx)) != 0) status_ = 1; 00223 #ifdef INEWTONDEBUG 00224 std::cerr << "iters = " << i << std::endl; 00225 #endif 00226 return approx; 00227 } 00228 00229 }; 00230 00231 template<class INT, class C> struct IntervalNewton 00232 { 00233 typedef interval_newton<INT, C> solver_t; 00234 typedef typename INT::boundary_type BT; 00235 typedef Seq<INT> solutions_t; 00236 00237 INT initial_approx; 00238 BT delta; 00239 int status; 00240 00241 IntervalNewton(INT approx_, BT delta_ = 0): initial_approx(approx_), delta(delta_), status(0) { }; 00242 00243 }; 00244 00245 template<class C, class M> struct solver_of; 00246 00247 template<class C, class IT> 00248 struct solver_of<C,IntervalNewton<IT,C> > { 00249 typedef IT Cell; 00250 typedef IntervalNewton<IT,C> Strategy; 00251 typedef Seq<IT> Solutions; 00252 }; 00253 00254 template<class POL, class IT> 00255 typename IntervalNewton<IT,typename POL::coeff_t>::solutions_t 00256 solve(const POL& f, IntervalNewton<IT,typename POL::coeff_t>& s) 00257 { 00258 typedef typename POL::coeff_t coeff_t; 00259 typedef typename IntervalNewton<IT,coeff_t>::solver_t solver_t; 00260 typename IntervalNewton<IT,coeff_t>::solutions_t sol; 00261 00262 POL df = diff(f); 00263 sol << solver_t::newton_iteration(f, df, s.initial_approx, s.delta, s.status); 00264 00265 return sol; 00266 } 00267 00268 00269 } //namespace mmx 00270 #endif // _realroot_solve_inewton_hpp_
