solver_uv_inewton_test.cpp

# include <realroot/realroot-config.hpp>
# ifdef REALROOT_WITH_NUMERIX

#include <iostream>
#include <iomanip>
#include <numerix/kernel.hpp>
#include <numerix/interval.hpp>
#include <realroot/Interval.hpp>
#include <realroot/ring_tens.hpp>
#include <realroot/solver_uv_interval_newton.hpp>

using namespace mmx;

typedef Numerix::integer ZZ;
typedef Numerix::floating RR;

//typedef interval<double> INT_double;
//typedef interval<RR> INT_mpfr;


typedef interval<double> INT_double;
typedef interval<Numerix::floating> INT_mpfr;

typedef polynomial< ring<ZZ,Monomials> > upolz_t;

upolz_t wilk(int n) {
    upolz_t w("x-1");
    upolz_t c(w);
    for (int j = 2; j <= n; j++) {
        c[0] = -j;
        w = w*c;
    }
    return w;
}


template<class INT, class POL>
void polytest(POL& p, INT approx, typename INT::boundary_type delta = typename INT::boundary_type(0))
{
  IntervalNewton<INT, typename POL::coeff_t> in(approx, delta);
  Seq<INT> sol;
  sol = solve(p, in);
  std::cout << sol << std::endl;
  std::cout << "sol status = " << in.status << std::endl;    
}

int main(int argc, char** argv) 
{
  std::cout.precision(16);
  std::cerr.precision(16);
//  mmx_bit_precision = 53;
  polynomial< ring<double, Monomials> > p("x^2+3.01*x+1.002"); // -2.628843850, -0.3811561496
  std::cout <<"p: "<<p<< std::endl;
  std::cout << "IntervalNewton: interval = double, poly = double" << std::endl;
  polytest<INT_double, polynomial< ring<double, Monomials> > >(p, INT_double(-1,0));
    
  std::cout << "*********************" << std::endl; 
  std::cout << "IntervalNewton: interval = double, poly = integer" << std::endl;
//  upoldse<ZZ> f( "x^10+21*x^5-3*x^2+3*x+1"); // -1.850230737, -0.2584591940
  upolz_t f = wilk(10); // try 20 for fun! (and define INEWTONDEBUG)

  std::cout << "f: " << f <<  std::endl; 
  polytest<INT_double, upolz_t>(f, INT_double(8.5,9.7));
  std::cout << "IntervalNewton: interval = mpfr, poly = integer" << std::endl;
//  mmx_bit_precision = 256;
//  mpfr_set_default_prec(256);
  polytest<INT_mpfr, upolz_t>(f, INT_mpfr(8.5,9.7));
                                      
//  std::cout << "IntervalNewton: interval = rational, poly = integer" << std::endl;
//  polytest<INT_rat, upoldse<integer> >(f, INT_rat(-1,0), rational(1)/rational(1E2));
  return 0;                                
}

# else
 int main(int argc, char** argv) { return 0; }
# endif //REALROOT_WITH_NUMERIX