PolynomialMV.m Mathsource document CONFTALK.DOC A Package for Manipulating Multivariable Polynomials George Hutchinson, National Institutes of Health, Bethesda, MD E-Mail: GAH@HELIX.NIH.GOV A talk given at the Mathematica Developer Conference (Champaign, IL), April, 1994, with some later additions. A related article "Manipulating Polynomials with Multiple Variables" has been submitted to The Mathematica Journal. Viewgraph 1. The Interface Problem for Symbolic Computation How can interactive symbolic algebra computation be made easier for the occasional user? - Large symbolic expressions cause obstacles to understanding and computation. - Computer specification of pencil-and-paper operations may be laborious or difficult to formulate. Viewgraph 2. An Example of Symbolic Computation In[1]:= mm = {{ka+kb*w, kf, 0, 0}, {-2*kb-2*kc, kc+kd+ke, -3*kf, 0}, {0, -kd, kd+I*xa, -u}, {0, 0, -I*u, kg}} Out[1]= {{ ka + kb w, kf, 0, 0}, {-2 kb - 2 kc, kc + kd + ke, -3 kf, 0}, {0, -kd, kd + I xa, -u}, {0, 0, -I u, kg}} In[2]:= DetFactored[mm] Out[2]:= 2 I (-2 kb - 2 kc) kf u - 3 kd kf kg (ka + kb w) - 2 I (kc + kd + ke) u (ka + kb w) - (-2 kb - 2 kc) kf kg (kd + I xa) + (kc + kd + ke) kg (ka + kb w) (kd + I xa) Commentary: The DetFactored package command provided a more compact representation of a symbolic determinant by replacing matrix elements which were sums by auxiliary variables, computing Det[mm], and then resubstituting the corresponding matrix elements for the auxiliary variables. Direct computation of the determinant leads to a fully expanded polynomial, as seen next. Viewgraph 3. In[3]:= Det[mm] Out[3]= 2 ka kc kd kg + ka kd kg + ka kd ke kg - 3 ka kd kf kg + 2 kb kd kf kg + 2 2 2 kc kd kf kg - I ka kc u - I ka kd u - 2 2 2 I ka ke u - 2 I kb kf u - 2 I kc kf u + 2 kb kc kd kg w + kb kd kg w + kb kd ke kg w - 3 kb kd kf kg w - 2 2 2 I kb kc u w - I kb kd u w - I kb ke u w + I ka kc kg xa + I ka kd kg xa + I ka ke kg xa + 2 I kb kf kg xa + 2 I kc kf kg xa + I kb kc kg w xa + I kb kd kg w xa + I kb ke kg w xa Viewgraph 4. A Package PolynomialMV.m for Interactive Manipulation of Polynomial with Multiple Variables - Oriented towards manipulation of polynomials into equivalent forms involving factorization and expansion of selected parts (applying the distributive law) - Uses subpolynomial lists to increase understanding, control and efficiency of the computational process for polynomials which may be very large - Is best adapted for computation with "low level" polynomials, which do not have deeply nested parentheses (levels 0, 1 and 2). Commentary: Polynomials like Det[mm] above are level zero because they have no parentheses, and DetFactored[mm] is level one because there are no parentheses nested within other parentheses. A polynomial like (2 x (y - a)^3 + 3 y)^2 + 3 a (y + 2 a) has level two because the parentheses around y - a are nested within the parentheses around 2 x (y - a)^3 + 3 y. Polynomials which have parentheses nested more deeply than level two are often difficult to read, and practical interactive polynomial manipulation is not usually directed towards such forms. Viewgraph 5. Polynomials for PolynomialMV.m - Exact number coefficients are recommended (Integer, Rational, or exact Complex) - Variables are restricted to labels or labels with integer indices a, x0A, y[-1], alpha[1,0,-3], ... - Polynomials must be recursively constructed from numbers and variables using Plus, Times and Power operations, with positive integer exponents in all Power operations. a[0]^2 x + (b0 - 3/7 uu^3)(3 x (s - a)^3 - I a) (3 a - 2 I)/5 + (1 + I) uu^2 v Commentary: Use of approximate number coefficients (Real or Complex with Real parts) may lead to inaccurate computations with PolynomialMV.m. Irrational constants like Sqrt[3] or 2 + ArcSin[1/5] can also lead to trouble. It may be possible to replace such numbers by new variables. An alternative approach is to replace p by Rationalize[ N[p], dx ] before the PolynomialMV.m analysis, where the error bound dx is chosen sufficiently small. In Mathematica, Variables[a + Sin[x]] will equal {a, Sin[x]}. In PolynomialMV.m, such generality is not too appropriate, and would add considerable complexity to the algorithms. So, the package uses a narrower definition for variables. Viewgraph 6. Term List Transformations - General Principles The "term list" strategy for PolynomialMV.m: A list of polynomials {p1, p2, ..., pm} is used to represent the sum p = p1 + p2 + ... + pm For large polynomials like matrix determinants, it is easier and more efficient to group summands into smaller subpolynomials for manipulation. Term list transformations are used to separate, regroup and sort subpolynomial lists. These subpolynomials are transformed by separate factorizations and expansions, then recombined. Viewgraph 7. Term List Transformation Arguments "termlist" A list of polynomials adding up to the polynomial under analysis "termspec" An argument specifying a sublist of termlist, as described below. integer k: k-th term of termlist list of distinct integers: term list positions of sublist elements polynomial structure predicate: a True/False function on polynomials, which selects the sublist of all True terms. default specifies all terms of termlist: {1, 2,..., m}, m = Length[termlist] Commentary: The termlist - termspec argument pairs provide a flexible and compact method for describing a sublist of a term list, giving easy operational access to parts of a large polynomial. Viewgraph 8. Gathering and Separating Terms SeparateTerms[termlist, termspec] Specified terms are replaced by a list of all summands of such terms. Unspecified terms are unchanged except for list order. In[4]:= SeparateTerms[ {1 + x, 2 + y, 3 + z}, {1,3} ] Out[4]= {1, x, 3, z, 2 + y} Commentary: The termspec {1,3} specifies terms 1 and 3 of the term list. The summands 1, x, 3, z of these terms begin the output list, and the unspecified term 2 + y is left unchanged at the end of the output. GatherTerms[termlist, termspec] Specified terms are added and become the new first term. Unspecified terms are unchanged. GatherTerms[ { 1, 3 + x, y^2 a, z, 4 a }, ExactDivideQ[#,a]& ] Out[5]= 2 {4 a + a y , 1, 3 + x, z} Commentary: In the above example, the termspec is the polynomial structure predicate ExactDivideQ[#,a]&, which is True iff the term is divisible by a. Since 4 a and y^2 a are the terms of the example divisible by a, these are added to get the new first term, and the other terms 1, 3 + x and z of the term list are left unchanged. Viewgraph 9. User-Supplied Term Transformation TransformTerms[termlist, termspec, transform] Returns termlist with the user-supplied transform function applied to all terms specified by termspec. Polynomial equivalence of transform[t] and t is not tested! In[6]:= TransformTerms[ {-1, b(x-a), (x-a)^2}, 3, Expand ] Out[6]= 2 2 {-1, b (-a + x), a - 2 a x + x } Commentary: The termspec is 3, so the third term of the termlist is expanded. Note that the user is responsible for ensuring that t and transform[t] are equivalent polynomials, so that the sum of terms in termlist will still represent the same large polynomial. Viewgraph 10. Introducing New Polynomial Parts AdjustTerms[termlist, termspec, poly] Returns termlist with poly subtracted from the m specified terms, and one new term m*poly added to the end of the term list to preserve the sum. In[7]:= AdjustTerms[ {1+w^3, 2b, 0}, {1,3}, w^3 ] Out[7]= 3 3 {1, 2 b, -w , 2 w } Commentary: Given a polynomial p, we may want to replace it by p - q + q for some polynomial q, in preparation for further manipulation. AdjustTerms does such operations for term lists, preventing the automatic cancellation of - q + q. In the example, w^3 is subtracted from terms 1 and 3, and a 4-th term 2 w^3 is appended to preserve the polynomial sum. Viewgraph 11. Permuting Order in Term Lists There are two term list transformations which change only the order of terms in the term list: RearrangeTerms[termlist, termspec] Returns termlist re-ordered so that the terms specified by termspec appear first, followed by the terms not specified. PermuteTerms[termlist, termspec] Returns termlist re-ordered by termspec interpreted as a cycle permutation (for example, swaps two terms) Commentary: In the first example, term list elements with Head not equal to Plus are put before those with Head Plus. In the second example, terms 3 and 5 are swapped. In[8]:= RearrangeTerms[ {2 x y^3, 3 + a^2, z, (a + y)^3 + 1, (u + 3 v)^2}, (Head[#]=!=Plus)& ] Out[8]= 3 2 2 3 {2 x y, z, (u + 3 v) , 3 + a , 1 + (a + y) } In[9]:= PermuteTerms[ {1, 4 x, 6 x^2, 4 x^3, x^4}, {3,5} ] Out[9]= 4 3 2 {1, 4 x, x, 4 x , 6 x } Viewgraph 12. Auxiliary Functions for Search and Display ShowTerms[termlist, termspec] Outputs printed list of specified terms. ShowExponents[termlist, termspec] Outputs printed list of exponents, for all variables, for all specified terms. In[10]:= ShowTerms[ {x^2 y, 3 a^2 + x, 2 y + 4}, Head[#]===Plus& ] Term 2 2 3 a + x Term 3 4 + 2 y In[11]:= ShowExponents[{4 a x^2 + a^2 x y^3, 3 a^2 + b}, 1] Exponents, term 1: {a -> 2, b -> 0, x -> 2, y -> 3} Commentary: If desired, ShowExponents can be given a specific list of variables for exponent computation, as a third argument. Viewgraph 13. Accessing Polynomial Parts PolynomialParts[poly, label] Recursively breaks down the polynomial poly into parts, with rules introduced for easy access. In[12]:= PolynomialParts[ 2 a (b - 2 c)^2 d^3 (e + 3), ww] Out[12]= 2 {ww[1] -> 2 a d^3, ww[2] -> (b - 2 c) , ww[3] -> 3 + e, ww[4] -> b - 2 c, ww[5] -> -2 c} Commentary: PolynomialParts breaks sums into their summands, breaks products into a monomial times any non-monomial factors, and returns the base of a power. It recursively breaks down currently listed parts into smaller parts, until no further breaking is possible. It does not retain small terms such as numbers, variables, or variables raised to a power. There is an exception if poly is a monomial, when PolynomialParts breaks poly into its factors. If the label argument is omitted, it produces just the list of polynomial parts, without the easy-access rules. Viewgraph 14. Testing Term Lists LocateTerms[termlist, polypredQ, indexlist] Outputs a list of integers corresponding to specified positions in the term list. Here indexlist is a list of integer positions corresponding to terms of termlist. If t is the k-th term of termlist and k is in indexlist, then k is retained in the output list if polypredQ[t] is True. In[13]:= LocateTerms[{u^2 y, 3 a - u, 2 x}, ExactDivideQ[#,x]&, {1, 2, 3} ] Out[13]= {3} Commentary: In the example, the indexlist is the list {1, 2, 3} of all the terms of termlist, and polypredQ is True for those terms which are exactly divisible by x. Since only the third term is divisible by x, LocateTerms returns the list {3}. This operation may assist users to prepare or examine complicated termspec arguments. Viewgraph 15. Miscellaneous Auxiliary Functions SumTerms[termlist, termspec] Returns the polynomial sum of the specified term list. PolynomialLevel[poly] Returns the polynomial level (maximum depth of nested parentheses) for poly, or returns -1 if poly is not a polynomial. DetFactored[matrix] Returns determinant of a symbolic matrix, the same as Det[matrix] except that products of matrix elements in the determinant are not expanded when matrix elements are sums of terms. Commentary: SumTerms[termlist] is the polynomial represented by termlist. The PolynomialLevel operation may be useful to formulate termspec arguments. An example of the DetFactored operation was shown in Viewgraphs 2-3. Viewgraph 16. Polynomial Structure Predicates - Basic Forms Polynomial structure predicates, evaluating to True or False, can be used as "termspec" arguments for many PolynomialMV.m operations. PolynomialMVQ[expr] True if (strict) polynomial - variables must be labels or indexed labels MonomialFactorQ[poly1] True if poly1 is a number, variable or power of a variable MonomialQ[poly1] True if poly1 is a monomial factor or a product of monomial factors CanonicalQ[poly1] True if poly1 is a monomial or a sum of monomials (level 0) Commentary: Many polynomial structure predicates are easily constructible directly from Mathematica operations. The package adds several that are useful for polynomial manipulation and not immediately definable. Viewgraph 17. Examples In[14]:= PolynomialMVQ[ Sin[x]^2 + 2 ] Out[14]= False In[15]:= MonomialFactorQ[ xa[0]^3 ] Out[15]= True In[16]:= MonomialQ[ 3 a xa^3 y[2]^2 ] Out[16]= True In[17]:= CanonicalQ[ 4 + a^2 + 3 a x^2 ] Out[17]= True In[18]:= CanonicalQ[ (x - 1)(x - 2) ] Out[18]= False Viewgraph 18. Polynomial Exact Division Predicate ExactDivideQ[poly1, poly2] Returns True if there exists a polynomial q such that poly1 equals q * poly2. In[19]:= ExactDivideQ[ (u - 2)(a^2 - b^2), a + b ] Out[19]= True Commentary: ExactDivideQ is implemented by computing PolynomialMVQ[ Cancel[poly1/poly2] ]. Viewgraph 19. Polynomial Part Searching Predicates FindPartQ[poly1, poly2] Returns True if poly1 or any polynomial part of poly1 is a nonzero numerical multiple of poly2 FindSubPartQ[poly1, poly2] Returns True if poly1 or any polynomial part of poly1 is exactly divisible by poly2. In[20]:= FindPartQ[(a + 2 x y^2)(b - 3 x), x y^2 ] Out[20]= True In[21]:= FindSubPartQ[ z (a^2 - 9 u^2) + 2/3 a z, a + 3 u ] Out[21]= True Commentary: In these operations, polynomial parts are computed using PolynomialParts[poly1]. Viewgraph 19. Polynomial Transformations - Exact Division ExactDivide[poly1, poly2] Returns q * poly2 if poly2 exactly divides poly1 with canonical quotient polynomial q. Returns poly1 if poly2 doesn't divide poly1. In[22]:= p1 = 4 a^3 + 12 b + a^4 x + 3 a b x + a^3 b y^2 + 3 b^2 y^2 ; In[23]:= ExactDivide[ p1, a^3 + 3 b ] Out[23]= 3 2 (a + 3 b) (4 + a x + b y ) Commentary: Replaces poly1 by an equivalent polynomial expressible as canonical polynomial q times poly2, when this is possible. Viewgraph 20. Forward and Reverse Distributivity AntiDistribute[poly1, poly2] Returns q * poly2 + r by factoring poly2 from all exactly-divisible summands of poly1 DistributeTerm[poly1, poly2] If poly1 is divisible by poly2 = t1 + t2 + ... + tm, say poly1 = q * poly2, the function returns: q*t1 + q*t2 + ... + q*tm In[24]:= p2 = 2 a (x - u)^2 + (4 x + 3) (x - u)^2 + 7; In[25]:= AntiDistribute[ p2, (x - u)^2 ] Out[25]= 2 7 + (-u + x) (3 + 2 a + 4 x) In[26]:= DistributeTerm[ (y - 1)(x - 2) z^3, x - 2 ] Out[26]= 3 3 -2 (-1 + y) z + x (-1 + y) z Commentary: Corresponds to paper-and-pencil factoring and expanding an expression using the distributivity law. Viewgraph 21. Controlled Expansion of Polynomials TopExpand[poly1] Expands to remove outermost parentheses (nested parentheses are left unchanged) ExpandTerm[poly1, poly2] If poly1 = q * poly2, expand poly2 and returns q unchanged: q * Expand[poly2] In[27]:= TopExpand[ ((s-x0)^2 + k(t - y0 + 1) - w)^2 ] Out[27]= 2 2 4 w - 2 w (s - x0) + (s - x0) - 2 2 k w (1 + t - y0) + 2 k (s - x0) (1 + t - y0) + 2 2 k (1 + t - y0) In[28]:= ExpandTerm[ w (u-w)^2 (u-z)(y-z), w (y-z) ] Out[28]= 2 (u - w) (u - z) (w y - w z) Commentary: These operations allow selective expansion of the arguments. Viewgraph 22. Factoring a Number FactorNumber[poly1, c] If poly1 is m * s^k for m a monomial, s = t1 + t2 + ... + tm, k a positive integer and c a nonzero number, this returns: c^k * m * (t1/c + t2/c + ... + tm/c)^k In[29]:= FactorNumber[ 3 a (2 b - 6)^3, -2 ] Out[29]= 3 -24 a (3 - b) Commentary: FactorNumber works only with the forms m * s^k, m * s, s^k and s. Otherwise, poly1 is returned unchanged. Viewgraph 23. Conclusions - The package PolynomialMV.m demonstrates that addition of simple tools can significantly reduce the problems and effort needed for interactive manipulation of large polynomials. - There are substantial interface problems still remaining for interactive computer algebra, for such tasks as: search for critical parts of lengthy output easy access to expression parts easy description of standard manipulations systematic control of the computation process Feedback? George Hutchinson Bldg. 12A, Room 2041 National Institutes of Health Bethesda, MD 20892 E-mail: GAH@HELIX.NIH.GOV