/* Boost example/newton-raphson.cpp * Newton iteration for intervals * * Copyright 2003 Guillaume Melquiond * * Distributed under the Boost Software License, Version 1.0. * (See accompanying file LICENSE_1_0.txt or * copy at http://www.boost.org/LICENSE_1_0.txt) */#include <boost/numeric/interval.hpp>#include <vector>#include <algorithm>#include <utility>#include <iostream>#include <iomanip>template <class I> I f(const I& x){ return x * (x - 1.) * (x - 2.) * (x - 3.) * (x - 4.); }template <class I> I f_diff(const I& x){ return (((5. * x - 40.) * x + 105.) * x - 100.) * x + 24.; }static const double max_width = 1e-10;static const double alpha = 0.75;using namespace boost;using namespace numeric;using namespace interval_lib;// First method: no empty intervalstypedef interval<double> I1_aux;typedef unprotect<I1_aux>::type I1;std::vector<I1> newton_raphson(const I1& xs) {  std::vector<I1> l, res;  I1 vf, vd, x, x1, x2;  l.push_back(xs);  while (!l.empty()) {    x = l.back();    l.pop_back();    bool x2_used;    double xx = median(x);    vf = f<I1>(xx);    vd = f_diff<I1>(x);    if (zero_in(vf) && zero_in(vd)) {      x1 = I1::whole();      x2_used = false;    } else {      x1 = xx - division_part1(vf, vd, x2_used);      if (x2_used) x2 = xx - division_part2(vf, vd);    }    if (overlap(x1, x)) x1 = intersect(x, x1);    else if (x2_used) { x1 = x2; x2_used = false; }    else continue;    if (x2_used)      if (overlap(x2, x)) x2 = intersect(x, x2);      else x2_used = false;    if (x2_used && width(x2) > width(x1)) std::swap(x1, x2);    if (!zero_in(f(x1)))      if (x2_used) { x1 = x2; x2_used = false; }      else continue;    if (width(x1) < max_width) res.push_back(x1);    else if (width(x1) > alpha * width(x)) {      std::pair<I1, I1> p = bisect(x);      if (zero_in(f(p.first))) l.push_back(p.first);      x2 = p.second;      x2_used = true;    } else l.push_back(x1);    if (x2_used && zero_in(f(x2)))      if (width(x2) < max_width) res.push_back(x2);      else l.push_back(x2);  }  return res;}// Second method: with empty intervalstypedef change_checking<I1_aux, checking_no_nan<double> >::type I2_aux;typedef unprotect<I2_aux>::type I2;std::vector<I2> newton_raphson(const I2& xs) {  std::vector<I2> l, res;  I2 vf, vd, x, x1, x2;  l.push_back(xs);  while (!l.empty()) {    x = l.back();    l.pop_back();    double xx = median(x);    vf = f<I2>(xx);    vd = f_diff<I2>(x);    if (zero_in(vf) && zero_in(vd)) {      x1 = x;      x2 = I2::empty();    } else {      bool x2_used;      x1 = intersect(x, xx - division_part1(vf, vd, x2_used));      x2 = intersect(x, xx - division_part2(vf, vd, x2_used));    }    if (width(x2) > width(x1)) std::swap(x1, x2);    if (empty(x1) || !zero_in(f(x1)))      if (!empty(x2)) { x1 = x2; x2 = I2::empty(); }      else continue;    if (width(x1) < max_width) res.push_back(x1);    else if (width(x1) > alpha * width(x)) {      std::pair<I2, I2> p = bisect(x);      if (zero_in(f(p.first))) l.push_back(p.first);      x2 = p.second;    } else l.push_back(x1);    if (!empty(x2) && zero_in(f(x2)))      if (width(x2) < max_width) res.push_back(x2);      else l.push_back(x2);  }  return res;}template<class T, class Policies>std::ostream &operator<<(std::ostream &os,                         const boost::numeric::interval<T, Policies> &x) {  os << "[" << x.lower() << ", " << x.upper() << "]";  return os;}int main() {  {    I1_aux::traits_type::rounding rnd;    std::vector<I1> res = newton_raphson(I1(-1, 5.1));    std::cout << "Results: " << std::endl << std::setprecision(12);    for(std::vector<I1>::const_iterator i = res.begin(); i != res.end(); ++i)      std::cout << "  " << *i << std::endl;    std::cout << std::endl;  }  {    I2_aux::traits_type::rounding rnd;    std::vector<I2> res = newton_raphson(I2(-1, 5.1));    std::cout << "Results: " << std::endl << std::setprecision(12);    for(std::vector<I2>::const_iterator i = res.begin(); i != res.end(); ++i)      std::cout << "  " << *i << std::endl;    std::cout << std::endl;  }}