Treating some variables as constants in a multivariate polynomial

[stla]

Hello,

I'm trying to understand the multivariate polynomials. I have a long polynomial with 5 variables x, y, z, a, b. In fact I consider that a and b are some constants, and then I want the coefficients of the monomials x^i y^j z^k in function of a and b.

I'm able to list the exponents of P(x,y,z,a,b) and their corresponding coefficient. For example, if I want the coefficient of x^1 y^3 z^1 in function of a and b, I can group all the exponents starting with (1, 3, 1). Let's say I found:

.......
exponent: (1,3,1,2,1)
coef: 4
......
exponent: (1,3,1,0,2)
coef: 4
......

Then I know that the coefficient of x^1 y^3 z^1 is 4*a^2*b + 4*b^2.

I could get what I want in this way I think. But isn't there a more direct way?

[stla]

Yes, I get what I wanted:

x^0y^4z^0: (x^2 + (-2)*x + 1)*y^2
x^0y^4z^2: ((-2)*x^2 + 2*x)*y
x^0y^4z^4: (x^2)
x^0y^6z^0: ((-2)*x + (-2))*y^2
x^0y^6z^2: (2*x)*y
x^0y^8z^0: y^2
x^1y^3z^1: ((-4)*x + 4)*y^2 + (4*x^2 + (-4)*x)*y
x^1y^3z^3: (4*x)*y + ((-4)*x^2)
x^1y^5z^1: (-4)*y^2 + (4*x)*y
x^2y^2z^0: (2*x^2 + (-2)*x)*y^2 + ((-2)*x^2 + 2*x)*y
x^2y^2z^2: ((-2)*x + 6)*y^2 + ((-2)*x^2 + (-8)*x)*y + (6*x^2)
x^2y^2z^4: (2*x)*y
x^2y^4z^0: ((-4)*x + (-2))*y^2 + ((-2)*x^2 + (-4)*x)*y
x^2y^4z^2: (-2)*y^2 + (10*x)*y + ((-2)*x^2)
x^2y^6z^0: 2*y^2 + (2*x)*y
x^3y^1z^1: ((-4)*x)*y^2 + (4*x^2 + 4*x)*y + ((-4)*x^2)
x^3y^1z^3: 4*y^2 + ((-4)*x)*y
x^3y^3z^1: (-4)*y^2 + (4*x^2)
x^4y^0z^0: (x^2)*y^2 + ((-2)*x^2)*y + (x^2)
x^4y^0z^2: ((-2)*x)*y^2 + (2*x)*y
x^4y^0z^4: y^2
x^4y^2z^0: ((-2)*x)*y^2 + ((-4)*x^2 + (-4)*x)*y + ((-2)*x^2)
x^4y^2z^2: (-2)*y^2 + (10*x)*y + ((-2)*x^2)
x^4y^4z^0: y^2 + (4*x)*y + (x^2)
x^5y^1z^1: ((-4)*x)*y + (4*x^2)
x^6y^0z^0: ((-2)*x^2)*y + ((-2)*x^2)
x^6y^0z^2: (2*x)*y
x^6y^2z^0: (2*x)*y + (2*x^2)
x^8y^0z^0: (x^2)

The variables of the right members should be a and b.

[stla]

Hello,

CGAL::IO::set_pretty_mode(std::cout); is very nice. But I want to capture the printed polynomial as a string. I do:

    std::stringstream buffer;
    buffer << value << std::endl;
    std::string coefab = buffer.str();

But that does not work. I get "P[1(0,P[2(2,2)])(1,P[1(1,2)])]\n" instead of the pretty polynomial.

[stla]

CGAL::IO::set_pretty_mode(buffer); works! Very nice.

[stla]

I've reproduced this polynomial expansion that I did with Julia. Found the same results :-) The speed is more or less the same (about one minute, while it takes 20 minutes with Python). Now it remains to use gmp coefficients to preserve the fractions (so far I used double in CGAL).

Here is my code. Maybe you will see something that can be improved (I'd like to beat Julia).

typedef CGAL::Polynomial_type_generator<double, 9>::Type                   Poly_9;
typedef CGAL::Polynomial_traits_d<Poly_9>                                  PT_9;
typedef std::pair<CGAL::Exponent_vector, PT_9::Innermost_coefficient_type> Monomial_9;
typedef CGAL::Polynomial_type_generator<double, 6>::Type                   Poly_6;
typedef CGAL::Polynomial_traits_d<Poly_6>                                  PT_6;
typedef std::tuple<int, int, int>                                          Expo3;
typedef std::map<Expo3, Poly_6>                                            XYZ6;

Rcpp::StringMatrix test4(
  Rcpp::IntegerMatrix Powers, Rcpp::NumericVector Coeffs
) {

  PT_9::Construct_polynomial construct_polynomial;

  int nterms = Coeffs.size();
  std::list<Monomial_9> terms;

  for(int i = 0; i < nterms; i++) {
    Rcpp::IntegerVector powers = Powers(Rcpp::_, i);
    terms.push_back(
      std::make_pair(
        CGAL::Exponent_vector(powers.begin(), powers.end()), 
        Coeffs(i)
      )
    );
  }
  Poly_9 P = construct_polynomial(terms.begin(), terms.end());

  std::list<Monomial_9> monoms;
  PT_9::Monomial_representation mrepr;
  mrepr(P, std::back_inserter(monoms));

  XYZ6 Result;
  Poly_6 w0    = PT_6::Shift()(Poly_6(1), 1, 0); // 
  Poly_6 sqrt3 = PT_6::Shift()(Poly_6(1), 1, 1); // 
  Poly_6 A     = PT_6::Shift()(Poly_6(1), 1, 2); // 
  Poly_6 B     = PT_6::Shift()(Poly_6(1), 1, 3); // 
  Poly_6 C     = PT_6::Shift()(Poly_6(1), 1, 4); // 
  Poly_6 D     = PT_6::Shift()(Poly_6(1), 1, 5); // 

  std::list<Monomial_9>::iterator it_monoms;
  for(it_monoms = monoms.begin(); it_monoms != monoms.end(); it_monoms++) {
    CGAL::Exponent_vector allPowers = (*it_monoms).first;
    double coef = (*it_monoms).second;
    Expo3 powersXYZ = {allPowers[0], allPowers[1], allPowers[2]};
    Poly_6 PAB = 
      coef * CGAL::ipower(w0, allPowers[3]) * CGAL::ipower(sqrt3, allPowers[4]) * 
      CGAL::ipower(A, allPowers[5]) * CGAL::ipower(B, allPowers[6]) * 
      CGAL::ipower(C, allPowers[7]) * CGAL::ipower(D, allPowers[8]);
    if(Result.count(powersXYZ) == 0) {
      Result.emplace(powersXYZ, PAB);
    } else {
      Result.at(powersXYZ) += PAB;
    }
  }

  int nmonoms = Result.size();
  Rcpp::StringMatrix POVRay(2, nmonoms);
  int i = 0;
  for(const auto& [key, value] : Result) {
    std::string powers = "xyz(" + 
      std::to_string(std::get<0>(key)) + ", " +
      std::to_string(std::get<1>(key)) + ", " +
      std::to_string(std::get<2>(key)) + "): ";
    std::stringstream buffer;
    CGAL::IO::set_pretty_mode(buffer);
    buffer << value << "," << std::endl;
    std::string coefab = buffer.str();
    Rcpp::StringVector monom = Rcpp::StringVector::create(powers, coefab);
    POVRay(Rcpp::_, i++) = monom;
  }

  return Rcpp::transpose(POVRay);
}

[stla]

Finally I don't need Julia nor CGAL: I found a way to get the result within a minute with R (with a package implemented in C++).