realroot: include/realroot/solver_uv_interval

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