Testing whether an expression is a polynomial in one or more variables can be done with the method
bool ex::is_polynomial(const ex & vars) const;
In the case of more than one variable, the variables are given as a list.
(x*y*sin(y)).is_polynomial(x) // Returns true.
(x*y*sin(y)).is_polynomial(lst(x,y)) // Returns false.
A polynomial in one or more variables has many equivalent representations. Some useful ones serve a specific purpose. Consider for example the trivariate polynomial 4*x*y + x*z + 20*y^2 + 21*y*z + 4*z^2 (written down here in output-style). It is equivalent to the factorized polynomial (x + 5*y + 4*z)*(4*y + z). Other representations are the recursive ones where one collects for exponents in one of the three variable. Since the factors are themselves polynomials in the remaining two variables the procedure can be repeated. In our example, two possibilities would be (4*y + z)*x + 20*y^2 + 21*y*z + 4*z^2 and 20*y^2 + (21*z + 4*x)*y + 4*z^2 + x*z.
To bring an expression into expanded form, its method
ex ex::expand(unsigned options = 0);
may be called. In our example above, this corresponds to 4*x*y + x*z + 20*y^2 + 21*y*z + 4*z^2. Again, since the canonical form in GiNaC is not easy to guess you should be prepared to see different orderings of terms in such sums!
Another useful representation of multivariate polynomials is as a
univariate polynomial in one of the variables with the coefficients
being polynomials in the remaining variables. The method
collect() accomplishes this task:
ex ex::collect(const ex & s, bool distributed = false);
The first argument to collect() can also be a list of objects in which
case the result is either a recursively collected polynomial, or a polynomial
in a distributed form with terms like c*x1^e1*...*xn^en, as specified
by the distributed flag.
Note that the original polynomial needs to be in expanded form (for the
variables concerned) in order for collect() to be able to find the
coefficients properly.
The following ginsh transcript shows an application of collect()
together with find():
> a=expand((sin(x)+sin(y))*(1+p+q)*(1+d));
d*p*sin(x)+p*sin(x)+q*d*sin(x)+q*sin(y)+d*sin(x)+q*d*sin(y)+sin(y)+d*sin(y)
+q*sin(x)+d*sin(y)*p+sin(x)+sin(y)*p
> collect(a,{p,q});
d*sin(x)+(d*sin(x)+sin(y)+d*sin(y)+sin(x))*p
+(d*sin(x)+sin(y)+d*sin(y)+sin(x))*q+sin(y)+d*sin(y)+sin(x)
> collect(a,find(a,sin($1)));
(1+q+d+q*d+d*p+p)*sin(y)+(1+q+d+q*d+d*p+p)*sin(x)
> collect(a,{find(a,sin($1)),p,q});
(1+(1+d)*p+d+q*(1+d))*sin(x)+(1+(1+d)*p+d+q*(1+d))*sin(y)
> collect(a,{find(a,sin($1)),d});
(1+q+d*(1+q+p)+p)*sin(y)+(1+q+d*(1+q+p)+p)*sin(x)
Polynomials can often be brought into a more compact form by collecting common factors from the terms of sums. This is accomplished by the function
ex collect_common_factors(const ex & e);
This function doesn't perform a full factorization but only looks for factors which are already explicitly present:
> collect_common_factors(a*x+a*y);
(x+y)*a
> collect_common_factors(a*x^2+2*a*x*y+a*y^2);
a*(2*x*y+y^2+x^2)
> collect_common_factors(a*(b*(a+c)*x+b*((a+c)*x+(a+c)*y)*y));
(c+a)*a*(x*y+y^2+x)*b
The degree and low degree of a polynomial can be obtained using the two methods
int ex::degree(const ex & s);
int ex::ldegree(const ex & s);
which also work reliably on non-expanded input polynomials (they even work on rational functions, returning the asymptotic degree). By definition, the degree of zero is zero. To extract a coefficient with a certain power from an expanded polynomial you use
ex ex::coeff(const ex & s, int n);
You can also obtain the leading and trailing coefficients with the methods
ex ex::lcoeff(const ex & s);
ex ex::tcoeff(const ex & s);
which are equivalent to coeff(s, degree(s)) and coeff(s, ldegree(s)),
respectively.
An application is illustrated in the next example, where a multivariate polynomial is analyzed:
{
symbol x("x"), y("y");
ex PolyInp = 4*pow(x,3)*y + 5*x*pow(y,2) + 3*y
- pow(x+y,2) + 2*pow(y+2,2) - 8;
ex Poly = PolyInp.expand();
for (int i=Poly.ldegree(x); i<=Poly.degree(x); ++i) {
cout << "The x^" << i << "-coefficient is "
<< Poly.coeff(x,i) << endl;
}
cout << "As polynomial in y: "
<< Poly.collect(y) << endl;
}
When run, it returns an output in the following fashion:
The x^0-coefficient is y^2+11*y
The x^1-coefficient is 5*y^2-2*y
The x^2-coefficient is -1
The x^3-coefficient is 4*y
As polynomial in y: -x^2+(5*x+1)*y^2+(-2*x+4*x^3+11)*y
As always, the exact output may vary between different versions of GiNaC or even from run to run since the internal canonical ordering is not within the user's sphere of influence.
degree(), ldegree(), coeff(), lcoeff(),
tcoeff() and collect() can also be used to a certain degree
with non-polynomial expressions as they not only work with symbols but with
constants, functions and indexed objects as well:
{
symbol a("a"), b("b"), c("c"), x("x");
idx i(symbol("i"), 3);
ex e = pow(sin(x) - cos(x), 4);
cout << e.degree(cos(x)) << endl;
// -> 4
cout << e.expand().coeff(sin(x), 3) << endl;
// -> -4*cos(x)
e = indexed(a+b, i) * indexed(b+c, i);
e = e.expand(expand_options::expand_indexed);
cout << e.collect(indexed(b, i)) << endl;
// -> a.i*c.i+(a.i+c.i)*b.i+b.i^2
}
ex quo(const ex & a, const ex & b, const ex & x);
ex rem(const ex & a, const ex & b, const ex & x);
compute the quotient and remainder of univariate polynomials in the variable ‘x’. The results satisfy a = b*quo(a, b, x) + rem(a, b, x).
The additional function
ex prem(const ex & a, const ex & b, const ex & x);
computes the pseudo-remainder of ‘a’ and ‘b’ which satisfies c*a = b*q + prem(a, b, x), where c = b.lcoeff(x) ^ (a.degree(x) - b.degree(x) + 1).
Exact division of multivariate polynomials is performed by the function
bool divide(const ex & a, const ex & b, ex & q);
If ‘b’ divides ‘a’ over the rationals, this function returns true
and returns the quotient in the variable q. Otherwise it returns false
in which case the value of q is undefined.
ex ex::unit(const ex & x);
ex ex::content(const ex & x);
ex ex::primpart(const ex & x);
ex ex::primpart(const ex & x, const ex & c);
return the unit part, content part, and primitive polynomial of a multivariate
polynomial with respect to the variable ‘x’ (the unit part being the sign
of the leading coefficient, the content part being the GCD of the coefficients,
and the primitive polynomial being the input polynomial divided by the unit and
content parts). The second variant of primpart() expects the previously
calculated content part of the polynomial in c, which enables it to
work faster in the case where the content part has already been computed. The
product of unit, content, and primitive part is the original polynomial.
Additionally, the method
void ex::unitcontprim(const ex & x, ex & u, ex & c, ex & p);
computes the unit, content, and primitive parts in one go, returning them
in u, c, and p, respectively.
The functions for polynomial greatest common divisor and least common multiple have the synopsis
ex gcd(const ex & a, const ex & b);
ex lcm(const ex & a, const ex & b);
The functions gcd() and lcm() accept two expressions
a and b as arguments and return a new expression, their
greatest common divisor or least common multiple, respectively. If the
polynomials a and b are coprime gcd(a,b) returns 1
and lcm(a,b) returns the product of a and b. Note that all
the coefficients must be rationals.
#include <ginac/ginac.h>
using namespace GiNaC;
int main()
{
symbol x("x"), y("y"), z("z");
ex P_a = 4*x*y + x*z + 20*pow(y, 2) + 21*y*z + 4*pow(z, 2);
ex P_b = x*y + 3*x*z + 5*pow(y, 2) + 19*y*z + 12*pow(z, 2);
ex P_gcd = gcd(P_a, P_b);
// x + 5*y + 4*z
ex P_lcm = lcm(P_a, P_b);
// 4*x*y^2 + 13*y*x*z + 20*y^3 + 81*y^2*z + 67*y*z^2 + 3*x*z^2 + 12*z^3
}
The resultant of two expressions only makes sense with polynomials. It is always computed with respect to a specific symbol within the expressions. The function has the interface
ex resultant(const ex & a, const ex & b, const ex & s);
Resultants are symmetric in a and b. The following example
computes the resultant of two expressions with respect to x and
y, respectively:
#include <ginac/ginac.h>
using namespace GiNaC;
int main()
{
symbol x("x"), y("y");
ex e1 = x+pow(y,2), e2 = 2*pow(x,3)-1; // x+y^2, 2*x^3-1
ex r;
r = resultant(e1, e2, x);
// -> 1+2*y^6
r = resultant(e1, e2, y);
// -> 1-4*x^3+4*x^6
}
Square-free decomposition is available in GiNaC:
ex sqrfree(const ex & a, const lst & l = lst());
Here is an example that by the way illustrates how the exact form of the result may slightly depend on the order of differentiation, calling for some care with subsequent processing of the result:
...
symbol x("x"), y("y");
ex BiVarPol = expand(pow(2-2*y,3) * pow(1+x*y,2) * pow(x-2*y,2) * (x+y));
cout << sqrfree(BiVarPol, lst(x,y)) << endl;
// -> 8*(1-y)^3*(y*x^2-2*y+x*(1-2*y^2))^2*(y+x)
cout << sqrfree(BiVarPol, lst(y,x)) << endl;
// -> 8*(1-y)^3*(-y*x^2+2*y+x*(-1+2*y^2))^2*(y+x)
cout << sqrfree(BiVarPol) << endl;
// -> depending on luck, any of the above
...
Note also, how factors with the same exponents are not fully factorized with this method.
Polynomials can also be fully factored with a call to the function
ex factor(const ex & a, unsigned int options = 0);
The factorization works for univariate and multivariate polynomials with rational coefficients. The following code snippet shows its capabilities:
...
cout << factor(pow(x,2)-1) << endl;
// -> (1+x)*(-1+x)
cout << factor(expand((x-y*z)*(x-pow(y,2)-pow(z,3))*(x+y+z))) << endl;
// -> (y+z+x)*(y*z-x)*(y^2-x+z^3)
cout << factor(pow(x,2)-1+sin(pow(x,2)-1)) << endl;
// -> -1+sin(-1+x^2)+x^2
...
The results are as expected except for the last one where no factorization seems to have been done. This is due to the default option factor_options::polynomial (equals zero) to factor(), which tells GiNaC to try a factorization only if the expression is a valid polynomial. In the shown example this is not the case, because one term is a function.
There exists a second option factor_options::all, which tells GiNaC to ignore non-polynomial parts of an expression and also to look inside function arguments. With this option the example gives:
...
cout << factor(pow(x,2)-1+sin(pow(x,2)-1), factor_options::all)
<< endl;
// -> (-1+x)*(1+x)+sin((-1+x)*(1+x))
...
GiNaC's factorization functions cannot handle algebraic extensions. Therefore the following example does not factor:
...
cout << factor(pow(x,2)-2) << endl;
// -> -2+x^2 and not (x-sqrt(2))*(x+sqrt(2))
...
Factorization is useful in many applications. A lot of algorithms in computer algebra depend on the ability to factor a polynomial. Of course, factorization can also be used to simplify expressions, but it is costly and applying it to complicated expressions (high degrees or many terms) may consume far too much time. So usually, looking for a GCD at strategic points in a calculation is the cheaper and more appropriate alternative.