calc: Decomposing Polynomials

 
 11.6.3 Decomposing Polynomials
 ------------------------------
 
 The ‘poly’ function takes a polynomial and a variable as arguments, and
 returns a vector of polynomial coefficients (constant coefficient
 first).  For example, ‘poly(x^3 + 2 x, x)’ returns ‘[0, 2, 0, 1]’.  If
 the input is not a polynomial in ‘x’, the call to ‘poly’ is left in
 symbolic form.  If the input does not involve the variable ‘x’, the
 input is returned in a list of length one, representing a polynomial
 with only a constant coefficient.  The call ‘poly(x, x)’ returns the
 vector ‘[0, 1]’.  The last element of the returned vector is guaranteed
 to be nonzero; note that ‘poly(0, x)’ returns the empty vector ‘[]’.
 Note also that ‘x’ may actually be any formula; for example,
 ‘poly(sin(x)^2 - sin(x) + 3, sin(x))’ returns ‘[3, -1, 1]’.
 
    To get the ‘x^k’ coefficient of polynomial ‘p’, use ‘poly(p,
 x)_(k+1)’.  To get the degree of polynomial ‘p’, use ‘vlen(poly(p, x)) -
 1’.  For example, ‘poly((x+1)^4, x)’ returns ‘[1, 4, 6, 4, 1]’, so
 ‘poly((x+1)^4, x)_(2+1)’ gives the ‘x^2’ coefficient of this polynomial,
 6.
 
    One important feature of the solver is its ability to recognize
 formulas which are “essentially” polynomials.  This ability is made
 available to the user through the ‘gpoly’ function, which is used just
 like ‘poly’: ‘gpoly(EXPR, VAR)’.  If EXPR is a polynomial in some term
 which includes VAR, then this function will return a vector ‘[X, C, A]’
 where X is the term that depends on VAR, C is a vector of polynomial
 coefficients (like the one returned by ‘poly’), and A is a multiplier
 which is usually 1.  Basically, ‘EXPR = A*(C_1 + C_2 X + C_3 X^2 +
 ...)’.  The last element of C is guaranteed to be non-zero, and C will
 not equal ‘[1]’ (i.e., the trivial decomposition EXPR = X is not
 considered a polynomial).  One side effect is that ‘gpoly(x, x)’ and
 ‘gpoly(6, x)’, both of which might be expected to recognize their
 arguments as polynomials, will not because the decomposition is
 considered trivial.
 
    For example, ‘gpoly((x-2)^2, x)’ returns ‘[x, [4, -4, 1], 1]’, since
 the expanded form of this polynomial is ‘4 - 4 x + x^2’.
 
    The term X may itself be a polynomial in VAR.  This is done to reduce
 the size of the C vector.  For example, ‘gpoly(x^4 + x^2 - 1, x)’
 returns ‘[x^2, [-1, 1, 1], 1]’, since a quadratic polynomial in ‘x^2’ is
 easier to solve than a quartic polynomial in ‘x’.
 
    A few more examples of the kinds of polynomials ‘gpoly’ can discover:
 
      sin(x) - 1               [sin(x), [-1, 1], 1]
      x + 1/x - 1              [x, [1, -1, 1], 1/x]
      x + 1/x                  [x^2, [1, 1], 1/x]
      x^3 + 2 x                [x^2, [2, 1], x]
      x + x^2:3 + sqrt(x)      [x^1:6, [1, 1, 0, 1], x^1:2]
      x^(2a) + 2 x^a + 5       [x^a, [5, 2, 1], 1]
      (exp(-x) + exp(x)) / 2   [e^(2 x), [0.5, 0.5], e^-x]
 
    The ‘poly’ and ‘gpoly’ functions accept a third integer argument
 which specifies the largest degree of polynomial that is acceptable.  If
 this is ‘n’, then only C vectors of length ‘n+1’ or less will be
 returned.  Otherwise, the ‘poly’ or ‘gpoly’ call will remain in symbolic
 form.  For example, the equation solver can handle quartics and smaller
 polynomials, so it calls ‘gpoly(EXPR, VAR, 4)’ to discover whether EXPR
 can be treated by its linear, quadratic, cubic, or quartic formulas.
 
    The ‘pdeg’ function computes the degree of a polynomial; ‘pdeg(p,x)’
 is the highest power of ‘x’ that appears in ‘p’.  This is the same as
 ‘vlen(poly(p,x))-1’, but is much more efficient.  If ‘p’ is constant
 with respect to ‘x’, then ‘pdeg(p,x) = 0’.  If ‘p’ is not a polynomial
 in ‘x’ (e.g., ‘pdeg(2 cos(x), x)’, the function remains unevaluated.  It
 is possible to omit the second argument ‘x’, in which case ‘pdeg(p)’
 returns the highest total degree of any term of the polynomial, counting
 all variables that appear in ‘p’.  Note that ‘pdeg(c) = pdeg(c,x) = 0’
 for any nonzero constant ‘c’; the degree of the constant zero is
 considered to be ‘-inf’ (minus infinity).
 
    The ‘plead’ function finds the leading term of a polynomial.  Thus
 ‘plead(p,x)’ is equivalent to ‘poly(p,x)_vlen(poly(p,x))’, though again
 more efficient.  In particular, ‘plead((2x+1)^10, x)’ returns 1024
 without expanding out the list of coefficients.  The value of
 ‘plead(p,x)’ will be zero only if ‘p = 0’.
 
    The ‘pcont’ function finds the “content” of a polynomial.  This is
 the greatest common divisor of all the coefficients of the polynomial.
 With two arguments, ‘pcont(p,x)’ effectively uses ‘poly(p,x)’ to get a
 list of coefficients, then uses ‘pgcd’ (the polynomial GCD function) to
 combine these into an answer.  For example, ‘pcont(4 x y^2 + 6 x^2 y,
 x)’ is ‘2 y’.  The content is basically the “biggest” polynomial that
 can be divided into ‘p’ exactly.  The sign of the content is the same as
 the sign of the leading coefficient.
 
    With only one argument, ‘pcont(p)’ computes the numerical content of
 the polynomial, i.e., the ‘gcd’ of the numerical coefficients of all the
 terms in the formula.  Note that ‘gcd’ is defined on rational numbers as
 well as integers; it computes the ‘gcd’ of the numerators and the ‘lcm’
 of the denominators.  Thus ‘pcont(4:3 x y^2 + 6 x^2 y)’ returns 2:3.
 Dividing the polynomial by this number will clear all the denominators,
 as well as dividing by any common content in the numerators.  The
 numerical content of a polynomial is negative only if all the
 coefficients in the polynomial are negative.
 
    The ‘pprim’ function finds the “primitive part” of a polynomial,
 which is simply the polynomial divided (using ‘pdiv’ if necessary) by
 its content.  If the input polynomial has rational coefficients, the
 result will have integer coefficients in simplest terms.