solver_uv_all_test.cpp

#include <iostream>
#include <iomanip>
#include <realroot/GMP.hpp>
#include <realroot/polynomial_tensor.hpp>
#include <realroot/polynomial_bernstein.hpp>

#include <realroot/solver_uv_bspline.hpp>
#include <realroot/solver_uv_continued_fraction.hpp>
#include <realroot/solver_uv_sleeve.hpp>

using namespace mmx;

typedef GMP::integer                            integer;
typedef GMP::rational                           rational;
typedef polynomial<integer,with<MonomialTensor> >        upol_i;
typedef polynomial<rational,with<MonomialTensor> >       upol_r;

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

upol_i mign(int a, int n, int m) {   
    upol_i f(1,n,0);
    upol_i c(a,1,0);

    c =c+ integer(1);
    c= c^m;
    c+=f;
    return c;
}

upol_i mign2(int a, int n, int m, int k) {   
    upol_i f(a,1,0);
    f=f+integer(1);
    upol_i c(f);
    f = f^k;
    c = c^m;
    c+= upol_i(1,n-k,0);

    return c*f;
}

void mirror(upol_i & f) {   

    for (unsigned i=1;i<f.size();i+=2)
        f[i]=-f[i];
}

template<class S, class P>
void test(P f)
{
  typedef typename P::Ring R;

  unsigned start, total;
  std::cout <<"\n----------- Input:\n"<<f<< std::endl;

  std::cout <<"\n* Using B-Splines: "<<std::endl;
  start = clock();
  std::cout <<"  First root in [2, 4]: "<< 
      solver<ring<double,Bernstein>, Bspline>::first_root(f,2,4)<<std::endl;
  std::cout <<"  First positive root : "<< 
      solver<ring<double,Bernstein>, Bspline>::first_root(f)<<std::endl;
  total= clock () - start;
  std::cout << "Computed in " << total << " ms"<< std::endl; 
  
  start= clock ();
  std::cout <<"\n* Using CF solver: "<<std::endl;
  std::cout <<"  First root in (2, 4) (isolate)    : "<< 
    solver<R,ContFrac<Isolate> >::template solve<S>(f,2,4)<<
    std::endl;
  std::cout <<"  First positive root (isolate)     : "<< 
      solver< R, ContFrac<Isolate> >::template solve<S>(f)<<
      std::endl;
  total= clock () - start;
  std::cout << "Computed in " << total << " ms"<< std::endl; 
  
  std::cout <<"  First root in (2, 4) (floor)      : "<< 
      solver< R, ContFrac<Isolate> >::template solve<S>(f,2,4)<<
      std::endl;
  std::cout <<"  First root in (2, 4) (approximate): "<< 
      solver< R, ContFrac<Isolate> >::template solve<S>(f,2,4)<<
      std::endl; 
  std::cout <<"  First positive root (floor)       : "<< 
      solver< R, ContFrac<Approximate> >::template solve<S>(f)<<
      std::endl;
   std::cout <<"  First positive root (approximate) : "<<
     solver< R, ContFrac<Approximate> >::template solve<S>(f)    <<
       std::endl;

   std::cout <<"\n* Positive real root(s) by CF : "<<
     solver< R, ContFrac<Isolate> >::template solve<S>(f)<<
     std::endl;
   
   std::cout <<"\n* Separate real root(s) by CF : "<<
     solver< R, ContFrac<Isolate> >::template solve<S>(f)<<
     std::endl;

   
   std::cout <<"\n* Real root(s) by Sleeve      : "
             <<solve(f, Sleeve<Approximate>())<<std::endl;
}


int main(int argc, char** argv) 
{
    upol_i f;
    upol_r g("22341334*x^8-2334457888*x^2-71232212*x+34554356");

    f=upol_i("22341334*x^8-2334457888*x^2-71232212*x+34554356");
    /* TODO << Compiling problems -- FIXME, by GL.
    //test<rational>(f);
    //test<rational>(g);

     f= mign(100,3,20);
     mirror(f);
     test<rational>(f);

     f= mign2(100,20,3,3);
     mirror(f);
     test<rational>(f);

     test<rational>( wilk(20) );
    
    //upol_t g("x0^2+x1^2-32");
    //std::cout <<  diff(g,1)   << std::endl;
    */

  return 0;
}