solver_uv_interval_newton.hpp

Go to the documentation of this file.

/*******************************************************************
 *   This file is part of the source code of the realroot kernel.
 *   Description: The Interval Newton method to approximate a root of a 
 *                polynomial given an initial interval that contains it.
 *   Author(s): 
 *        G. Tzoumas <geotz@di.uoa.gr>
 *           Department of Informatics and Telecommunications
 *           University of Athens (Greece).
 ********************************************************************/
#ifndef _realroot_solve_inewton_hpp_
#define _realroot_solve_inewton_hpp_

//#include <numerix/interval.hpp>
//#include <realroot/Interval.hpp>
#include <realroot/polynomial_tensor.hpp>
#include <realroot/Seq.hpp>

namespace mmx {
  
/*  // perhaps it might come handy in the future -- I leave it in comments
    namespace inewton {
        
      template<class C> struct IntervalTraits { };
      
      template<> template<class B> struct IntervalTraits<Interval<B> > {
          typedef Interval<B> interval_type;
          typedef B boundary_type;
          
          static inline boundary_type upper(const interval_type& x) { return x.upper(); }
          static inline boundary_type lower(const interval_type& x) { return x.lower(); }
          static inline boundary_type width(const interval_type& x) { return x.width(); }
          static inline boundary_type median(const interval_type& x) { return x; }
          static inline int sign(const interval_type& x) { return ::mmx::sign(x); }
          static inline bool intersect(interval_type &r, const interval_type &a, const interval_type &b) {
              return ::mmx::intersect(r, a, b);
          }
      };
      
      template<> template<class B> struct IntervalTraits<interval<B> > {
          typedef interval<B> interval_type;
          typedef B boundary_type;
          
          static inline boundary_type upper(const interval_type& x) { return ::mmx::upper(x); }
          static inline boundary_type lower(const interval_type& x) { return ::mmx::lower(x); }
          static inline boundary_type width(const interval_type& x) { return upper(x)-lower(x); }
          static inline boundary_type median(const interval_type& x) { return ::mmx::center(x); }
          static inline int sign(const interval_type& x) { return ::mmx::sign(x); }
          static inline bool intersect(interval_type &r, const interval_type &a, const interval_type &b) {
              if ((upper(a) < lower(b)) || (lower(a) > upper(b))) return false;
              r = interval_type(std::max(lower(a), lower(b)), std::min(upper(a), upper(b)));
              return true;
          }
      };
      
  } */

  template<class INT, class CT> struct interval_newton
  {
    typedef typename ring<CT, MonomialTensor>::Polynomial   polynomial;
    typedef typename INT::boundary_type         BT;
      
      static inline BT width(const INT& x) { return mmx::upper(x) - mmx::lower(x); }
      static inline BT median(const INT& x) { return mmx::lower(x) + width(x)/BT(2); }
      static bool singleton(const INT& x) { return mmx::lower(x) == mmx::upper(x); }

      static INT bisection_iteration(const polynomial& p, const INT approx, int& status_)
      {
          BT x = median(approx);
          INT lval, rval, mval;
          lval = p(INT(mmx::lower(approx)));
          rval = p(INT(mmx::upper(approx)));
          mval = p(INT(x));
          //    std::cerr << lval << ' ' << rval << ' ' << mval << ' ' << x << std::endl;
          //    assert( mmx::sign(lval)*mmx::sign(rval) <= 0 );
                   
          status_ = 0;
          if (mval == 0) return INT(x);
          if (lval == 0) return INT(mmx::lower(approx)); 
          if (rval == 0) return INT(mmx::upper(approx));
                   
          if (mmx::sign(lval) == 0 || mmx::sign(rval) == 0) {
              status_ = -3; return approx;
          }

          if (mmx::sign(lval) == mmx::sign(rval)) {
              status_ = -1; return approx;
          }
          
          if (mmx::sign(mval) == 0) return approx;
          
          if (mmx::sign(lval) == mmx::sign(mval)) return INT(x,mmx::upper(approx));
          return INT(mmx::lower(approx), x);
      }
      
      /*
       status:
       -4 = no convergence (should not happen!)
       -3 = precision not enough
       -2 = no root exists
       -1 = no root found (cannot guarantee existence)
        0 = root found (cannot guarantee uniqueness)
        1 = unique root found
       */
      // caution with coeffs including zero!
      static INT newton_iteration(const polynomial& p, const polynomial& dp, 
                                       INT approx, BT delta_, int& status_) 
      { 
          INT slope, testapprox, newapprox, ver;
          BT oldw = width(approx);
          BT neww;
          int z;
          
#ifdef INEWTONDEBUG
          int i = 0;
#endif
          
          INT lval, rval;
          
          lval = p(INT(mmx::lower(approx)));
          rval = p(INT(mmx::upper(approx)));
#ifdef INEWTONDEBUG
          std::cerr << "approx = " << approx << " lval = " << lval << " rval = " << rval << std::endl;
#endif
          if (lval == 0) { status_ = 0; return INT(mmx::lower(approx)); }
          if (rval == 0) { status_ = 0; return INT(mmx::upper(approx)); }
          if (mmx::sign(lval) == 0 || mmx::sign(rval) == 0) {
              status_ = -3; return approx;
          }
          
          status_ = 0;
          slope = dp(approx);
          
#ifdef INEWTONDEBUG
          std::cerr << "slope = " << slope << width(slope) << std::endl;
#endif

          z = mmx::sign(slope);
          if (mmx::sign(lval) == mmx::sign(rval)) {
              if (z == 0) status_ = -1; else status_ = -2;
              return approx;
          }

#ifdef INEWTONDEBUG
          std::cerr << "solving " << p << std::endl;
#endif

          /* INTERVAL NEWTON MAIN LOOP */
          while (1) {
              // should we always do simple bisection with intervals of length >= 1 ?
              // if (oldw >= 1 || z == 0) {
              
              // slope contains zero
              if (z == 0) {
                  approx = bisection_iteration(p, approx, status_);
                  neww = width(approx);
#ifdef INEWTONDEBUG
                  i++;
                  std::cerr << '(' << i << ") bisect: " << approx << ' ' << neww << std::endl;
#endif
                  if (status_ < 0) return approx;
                  if (neww < delta_) break;
                  assert(neww <= oldw);
                  if (neww == oldw) break;
                  oldw = neww;
                  
              } else {
              // slope doesn't contain zero

                  BT x = median(approx);
                  if (x == mmx::lower(approx) || x == mmx::upper(approx)) break;
#ifdef INEWTONDEBUG
                  std::cerr << "mid = " << x << " mid_INT = " << INT(x) << std::endl;
                  std::cerr << "p(x) = " << p(INT(x)) << std::endl;
                  std::cerr << "p'(approx) = " << dp(approx) << std::endl;
                  std::cerr << "p/p' = " << INT(p(INT(x))/dp(approx)) << std::endl;
#endif
                  testapprox = INT(x) - INT(p(INT(x))/dp(approx)); // Newton approximation

#ifdef INEWTONDEBUG
                  i++;
#endif
                  if (mmx::intersect(newapprox, testapprox, approx)) {
#ifdef INEWTONDEBUG
                      std::cerr <<"new =" << newapprox << " test = " << testapprox << " old = " << approx << std::endl;
#endif
                      neww = width(newapprox);
                      
                      // new approximation SHOULD contain the root
                      // if you want to verify, use the debugging code below
#ifdef INEWTONDEBUG
                      INT ver = p(newapprox);
                      if (ver != INT(0)) {
                          assert (mmx::sign(ver) == 0);
                      }
                      
                      std::cerr << '(' << i << ") Newton: " << newapprox << ' ' << neww << std::endl;
#endif
                      //if (Float(2)*neww >= oldw) break;
                      if (neww >= oldw) break;
                      approx = newapprox;
                      if (neww < delta_) break;
                      oldw = neww;
                  } else {
                      // no convergence
                      status_ = -4;
                      return approx;
                  }
              }
              slope = dp(approx);
#ifdef INEWTONDEBUG
              std::cerr << "slope = " << slope << width(slope) << std::endl;
#endif
              z = mmx::sign(slope);
          }
          if (mmx::sign(dp(approx)) != 0) status_ = 1;
#ifdef INEWTONDEBUG
          std::cerr << "iters = " << i << std::endl;
#endif
          return approx;
      }
      
  };
  
  template<class INT, class C> struct IntervalNewton
  {
      typedef interval_newton<INT, C> solver_t;
      typedef typename INT::boundary_type BT;
      typedef Seq<INT> solutions_t;
    
      INT initial_approx;
      BT delta;
      int status;
      
      IntervalNewton(INT approx_, BT delta_ = 0): initial_approx(approx_), delta(delta_), status(0) { };
      
  };
  
  template<class C, class M> struct solver_of;

  template<class C, class IT>  
  struct solver_of<C,IntervalNewton<IT,C> > {
    typedef IT Cell;
    typedef IntervalNewton<IT,C> Strategy;
    typedef Seq<IT> Solutions;
  };

  template<class POL, class IT>
  typename IntervalNewton<IT,typename POL::coeff_t>::solutions_t 
  solve(const POL& f, IntervalNewton<IT,typename POL::coeff_t>& s)
  {
      typedef typename POL::coeff_t coeff_t;
      typedef typename IntervalNewton<IT,coeff_t>::solver_t solver_t;
      typename IntervalNewton<IT,coeff_t>::solutions_t sol;
          
      POL df = diff(f); 
      sol << solver_t::newton_iteration(f, df, s.initial_approx, s.delta, s.status);
          
      return sol;
  }
  

} //namespace mmx 
#endif // _realroot_solve_inewton_hpp_