polynomial_mv_sparse_test.cpp

# include <sstream>
# include <iomanip>
# include <realroot/GMP.hpp>
# include <realroot/polynomial_sparse.hpp>

  
int main(int argc, char** argv) 
{
  using namespace mmx;
  using std::cout;
  using std::endl;

  typedef polynomial< double, with<Sparse,LexRevDegree> >  Poly_t;
  typedef Poly_t::Monomial                             Monomial_t;
  //Poly( C, d, i) means C * x_i^d
  Poly_t f = Poly_t( 1,0,0) + Poly_t( 1,1,1) + Poly_t(1,1,2); 
  // + Poly_t(1,1,3);
  
  std::cout << "f   : " << f << std::endl; 
  std::cout << "f^2 : " << (f^2) << std::endl; 
  Poly_t g=f^2;
  Poly_t h=g;
  

  Poly_t k = Poly_t( "x^2+0.00000001*x+100.000000001" ); 

  std::ostringstream os;  os <<std::fixed<< std::setprecision(16)<< k;
  std::cout << "k   : " << os.str() << std::endl; 
  std::cout << "\nk   : " << Poly_t(os.str().data()) << std::endl; 

  Seq<Monomial_t> seq;
  for ( Poly_t::iterator it = g.begin(); it != g.end(); ++it) {
      it->set_coeff( 1);
      seq << *it;      
      //  std::cout << "it: " << *it << std::endl; 
  }


  typedef polynomial< GMP::rational, with<Sparse,DegRevLex> >  Polynomial;
  
  Polynomial P("x^2*y+y^2-1"), Q("x-1");
  cout<<P<<" "<<Q<<endl;
  cout<<"P/Q: "<<P/Q<<endl;
  cout<<"P%Q: "<<P%Q<<endl;
  cout<<(P*P)/P<<endl;

  // std::cout << "seq: " << seq << std::endl; 
  // std::cout << "g[2]: " << g[1] << std::endl; 
  std::cout << "g: " << g << std::endl; 
  std::cout << "h: " << h << std::endl; 

  int dg = degree( h);
  int nbvars = h.nbvar();
  
  std::cout << "dg: " << dg << std::endl; 
  std::cout << "#vars: " << nbvars << std::endl; 

  h=0;  std::cout << "h: " << h << std::endl; 
  dg = degree( h);
  nbvars = h.nbvar() ;
  std::cout << "dg: " << dg << std::endl; 
  std::cout << "#vars: " << nbvars << std::endl; 

  //int dg_fact;
  double dg_fact = factorial<double>(dg);
  
  for (unsigned i = 0; i < seq.size(); ++i) {
      for (unsigned j = 0; j < seq.size(); ++j) {
          Monomial_t tmon = seq[i] * seq[j];
          
          for ( Poly_t::iterator it = h.begin(); it != h.end(); ++it) {
              if( IsComparable( *it, tmon) ) {
                  double c = dg_fact; 
                  
                  if( tmon.size() > 0) {
                      for (int k = 1; k <= nbvars; ++k) {
                        int t = ((unsigned) k < tmon.size()? tmon[k] : 0);
                        c = c / factorial<double>(t);
                      }
                      c = c / factorial<double>( dg - degree( tmon));
                      
                      std::cout << "t: " << tmon << " " 
                                << c << "  "
                                << it->coeff() << " "
                                << ( it->coeff() / c) << " "                      
                                << std::endl; 
                  }
              }
          }
      }
  }
  
  return 0;
  

//   q=p*p; std::cout<<q<<std::endl;
}