(* DIXON RESULTANTS - The Mathematica Package. Version $1.0$. *) (* :Name: Dixon` *) (* :Context: Dixon` *) (* :Author: George Nakos and Robert M. Williams. *) (* :Mathematica Version: 2.2 *) (* :Package Version: 1.0 *) (* :Title: Elimination with the Dixon Resultant. *) (* :Copyright: *) (* :Summary: This package computes the Dixon resultant, DixonResultant, which can be used to eliminate a number of unknowns out of a system of polynomial equations. For large systems DixonResultant is a much better eliminant than Mathematica's command Eliminate. So, it can be used to enhance Eliminate. DixonResultant also enhances Mathematica's Resultant, which implements the Sylvester resultant. This is up to constant factor(s) a special case of the Dixon resultant. There are two advantages in using DixonResultant over Resultant: a) The end matrix has smaller size; hence it easier to row reduce it or compute its determinant. b) A whole block of variables can be eliminated in one calculation, automatically, instead of the successive elimination of one variable at a time. As a by-product of our program we offer a symbolic Gauss elimination that complements Mathematica's RowReduce command and may also be of general interest. *) (* :Discussion: This program computes a) the CLASSICAL Dixon Resultant (ClassicalDixonResultant) and b) the improvement made by Kapur, Saxena and Yang (DixonResultant). (c.f. "Algebraic and Geometric Reasoning using Dixon Resultants" Proc. ACM ISSAC, Oxford, England, July 1994.) The method is primarily an elimination method that can also be used to solve systems of polynomial equations. It applies to any system of n+1 equations with n unknowns, parametric coefficients allowed. (Thus, it applies to systems with numerical coefficients and at least as many unknowns as equations; the extra unknowns are declared as parameters.) ClassicalDixonResultant is only available as an illustration of Dixon's approach. If it is identically zero it yields no information on the solutions of the polynomial system. In contrast, DixonResultant, which is our main function, can used even when ClassicalDixonResultant is zero; provided a certain precondition holds. *) (* :History: Alpha test version by George Nakos and Robert M. Williams, July 94. Comments by George Nakos and Janet McShane, July 94. Modified by George Nakos on: 8 May and 28 July 1995. Modified by George Nakos on: 9,10 September, 1995. *) (* :Keywords: Resultants, Dixon Resultant, Polynomial Equations, Linear Systems and Determinants. *) (* :Warning: Due to the intensity of the computations some calculations may be slow. *) (* :Limitation: None known. *) BeginPackage["Dixon`"] ; DixonSub::usage = "DixonSub[polys, xlist, alist] forms the matrix obtained by \n successively substituting the auxiliary variable list alist \n into the polynomial list polys with variables xlist. \n See also DixonPolynomial"; DixonPolynomial::usage = "DixonPolynomial[polys, xlist, alist] computes the Dixon polynomial \n of polys in variables xlist with auxiliary variables alist. \n polys is a list of polynomials. xlist and alist are lists of \n variables. If the length of polys is n then the lengths of xlist \n and alist should be n+1. See also DixonSub."; DixonMatrix::usage = "DixonMatrix[polys, xlist, alist] computes the Dixon matrix \n of polys in variables xlist with auxiliary variables alist. \n polys is a list of polynomials. xlist and alist are lists of \n variables. If the length of polys is n then the lengths of xlist \n and alist should be n+1. See also DixonSub and DixonPolynomial."; ClassicalDixonResultant::usage = "ClassicalDixonResultant[polys, xlist, alist] computes the classical \n Dixon resultant of polys in variables xlist with auxiliary variables \n alist. polys is a list of polynomials. xlist and alist are lists of \n variables. If the length of polys is n then the lengths of xlist \n and alist should be n+1. See also DixonSub, DixonPolynomial and \n DixonMatrix."; DixonReduction::usage = "DixonReduction[polys, xlist, alist] performs Gauss elimination \n on the Dixon matrix of polys in variables xlist with auxiliary \n variables alist. polys is a list of polynomials. xlist and alist \n are lists of variables. If the length of polys is n then the lengths \n of xlist and alist should be n+1. See also DixonMatrix and \n GaussElimination."; DixonResultant::usage = "DixonResultant[polys, xlist, alist] computes the Kapur-Saxena-Yang \n Dixon resultant of polys in variables xlist with auxiliary variables \n alist. polys is a list of polynomials. xlist and alist are lists of \n variables. If the length of polys is n then the lengths of xlist and \n alist should be n+1. See also DixonMatrix and ClassicalDixonResultant."; PreconditionQ::usage = "PreconditionQ[polys, xlist, alist] tests for the validity of the \n Kapur-Saxena-Yang precondition for polys in variables xlist with \n auxiliary variables alist. polys is a list of polynomials. xlist and \n alist are lists of variables. If the length of polys is n then the \n lengths of xlist and alist should be n+1. \n The test is probabilistic in the case of success and deterministic \n in the case of failure. See also DixonMatrix and DixonResultant."; GaussElimination::usage = "GaussElimination[m] performs symbolic Gauss Elimination on the \n matrix m (downward pass only)."; ProductLeadingEntries::usage = "ProductLeadingEntries[m] returns the product of the leading entries \n of matrix m."; GeneralizedDeterminant::usage = "GeneralizedDeterminant[m] returns the product of the leading entries \n of an echelon form of the matrix m. See also GaussElimination and \n ProductLeadingEntries"; Begin["`Private`"] ; (* The next two procedures compute the Dixon polynomial *) (* (DixonPolynomial). *) (* Arguments for DixonSub and DixonPolynomial consist of a list of n+1 *) (* polynomials in n variables (polys), a list of the n original variables *) (* (xlist), and a list of the n auxiliary variables (alist). *) DixonSub[polys_List, xlist_List, alist_List] :=Block[ {i,dix={polys}}, For[i=1, i<=Length[xlist], i++, AppendTo[dix, Last[dix] /. xlist[[i]]->alist[[i]] ]] ; If[Length[dix]===Length[Transpose[dix]], dix, "The variable list is wrong."] ] (* Note: The Dixon polynomial is the determinant of DixonSub divided by *) (* (xlist[[1]]-alist[[1]])(xlist[[2]]-alist[[2]]) ... *) (* (xlist[[n]]-alist[[n]]) However, an efficient way of computing it is by *) (* diving the difference of any two consecutive rows, say the (i+1)th minus*) (* the ith row, by (xlist[[i]]-alist[[i]]). This way, the determinant is *) (* kept small and there is no need for the very costly factorization of *) (* Det[DixonSub]. This is exactly what DixonPolynomial does. *) DixonPolynomial[polys_List, xlist_List, alist_List] :=Block[ {i,dix=DixonSub[polys, xlist, alist], lis}, If[Head[dix]===String,dix, lis = {dix[[1]]}; For[i=1, i<=Length[xlist], i++, AppendTo[lis, Together[(dix[[i+1]]-dix[[i]]) / (xlist[[i]]-alist[[i]]) ] ] ]; Det[lis] ] ] (* MaxDegVec takes a list of polynomials, polys, and a list of variables, *) (* xlist and returns a list of the maximum degrees of the variables. *) MaxDegVec[polys_List, xlist_List] := Max[Exponent[polys,#]]& /@ xlist (* Computes all lists with entries less than or equal to the *) (* corresponding entries of a given list of integers. Zeros are included *) (* with the exception of the zero vector. *) LessLists[list_List] := Drop[Flatten[Array[List, list+1]-1, Length[list]-1],1] (* n-dimensional zero vector *) ZeroVector[n_Integer] := 0 Range[n] (* Note: 0 Range[n] is faster than Table[0,{i,1,n}]. *) (* The substitution list1->list2 for equal length lists. *) ListSubstitution[list1_List, list2_List] := Table[list1[[i]]->list2[[i]],{i,1,Length[list1]}] (* The following computes the constant term of the multivariate *) (* polynomial poly in the variables vars. *) ConstantTerm[poly_, vars_List] := poly /. ListSubstitution[vars, ZeroVector[Length[vars]]] (* MonomialBasis lists all monomials out of the list of variables vars *) (* with each variable having max degree from the list maxdegrees. This *) (* function may be of independent interest in other projects. *) MonomialBasis[vars_List, maxdegrees_List] := (Times @@ #&) /@ (vars^# & /@ LessLists[maxdegrees]); (* DixonMatrix computes the classical Dixon Matrix. *) DixonMatrix[polys_List, xlist_List, alist_List] := Block[ {i, dixonpoly = Expand[DixonPolynomial[polys, xlist, alist]], totallist = Join[xlist, alist], lesslists = LessLists[MaxDegVec[{dixonpoly}, totallist]], allsubstitutions, allproducts, coeffmatrix, separator }, If[ dixonpoly === 0,0, If[Head[dixonpoly]===String,dixonpoly, separator = 1 + Length[LessLists[MaxDegVec[{dixonpoly}, alist]]]; allsubstitutions = ListSubstitution[totallist, #]& /@ lesslists; allproducts = (Times @@ #&) /@ (totallist^# & /@ lesslists); coeffmatrix = Coefficient[dixonpoly, #]& /@ allproducts; coeffmatrix = Table[coeffmatrix[[i]] /. allsubstitutions[[i]],{i,1,Length[allproducts]}]; PrependTo[coeffmatrix, ConstantTerm[dixonpoly, totallist]]; Partition[coeffmatrix, separator] // Transpose ] ] ] (* Note: The zero vector was excluded and then added later *) (* since we had to compute the constant term first. Mathematica *) (* does not understand Coefficient[poly, x^0 y^0]. *) (* Note: MonomialBasis was not called for allproducts *) (* in DixonMatrix since LessLists would be computed twice. *) (* ClassicalDixonResultant computes the Dixon Resultant, i.e., *) (* the determinant of the Dixon Matrix. *) (* Note: The number of polynomials, polys, MUST be one more than the *) (* number of variables, vars, otherwise the determinant cannot be computed.*) ClassicalDixonResultant[polys_List, xlist_List, alist_List] := Det[DixonMatrix[polys, xlist, alist]] (* The Kapur - Saxena - Yang approach. *) (* Also available, testing for the precondition, separately. *) (* First we deal with GAUSS ELIMINATION *) (* GAUSSIAN ELIMINATION suited for symbolic calculations. The downward *) (* pass only. *) (* Program computing a row echelon form of a matrix without changing the *) (* determinant. The elementary row operations doing that are one of the *) (* following two types. *) (* 1) For interchanging two rows: The sign of one row is flipped. *) (* 2) For adding a multiple of a row to another: the row that gets *) (* replaced is used with coefficient one. *) (* The program also computes the product of the leading entries *) (* of the thus computed row echelon form. *) (* Testing for zero vector. *) ZeroVectorQ[list_List] := list === Table[0, {Length[list]}] (* Testing for zero matrix. *) ZeroMatrixQ[list_List] := list === Table[0,{Length[list]},{Length[list[[1]]]}] (* Swaps rows Rowi and Rowj . *) SwapRow[a_,i_,j_] := If[i==j,a, ReplacePart[ReplacePart[a,a[[j]],i],-a[[i]],j]] (* Replace Rowj with: Rowj + scal * Rowi . *) AddRow[a_,j_,i_,scal_] := ReplacePart[a,Expand[a[[j]] + scal a[[i]]],j] (* The forward pivot. *) PivotFor[a_,i_,j_] := Block[{k,aa=a}, If[a[[i,j]] =!= 0, For[k=i+1,k<=Length[a],k++, aa = AddRow[aa,k,i,-aa[[k,j]]/aa[[i,j]]] ] ]; aa ] Unprotect[Transpose] Transpose[{}]={} (* Drops Rowi and Columnj . *) MinorMat[a_,i_,j_] := Transpose[Drop[Transpose[Drop[a,{i}]],{j}]] (* Drops block past Rowi and Columnj. *) MinorColBlock[a_,i_,j_] := Transpose[Drop[Transpose[Drop[a,i]],j]] (* Finds current pivot. First pivot in each sub-block. *) FindFirstPivot[a_] := Block[{aa=Transpose[a],j=1,i=1}, If[ZeroMatrixQ[aa], -1, While[ZeroVectorQ[aa[[j]]],j++]; While[a[[i,j]] === 0,i++]; {i,j}]] GaussElimination[a_] := Block[{fp,aa=a,mm=aa,top=0,t1=0,t2=0,bb }, If[Length[a] === 1, a, While[Length[mm] =!= 1 && (fp = FindFirstPivot[mm]) =!= -1, fp +={top,t2}; top++; bb=SwapRow[aa,top,fp[[1]]]; aa = bb; bb=PivotFor[aa,top,fp[[2]]]; aa = bb; mm=MinorColBlock[aa,top,fp[[2]]]; {t1,t2}=fp ]; Together[aa] ] ] (* Sometimes the output of GaussElimination may need expansion, by *) (* Expand. Only the Together command was used because rational expressions *) (* are more readable in this form. *) ProductLeadingEntries[m_List] := Block[{mm=DeleteCases[m,0,2]}, mm=DeleteCases[mm,{}]; Product[mm[[i,1]],{i,1,Length[mm]}] // Expand // Together] GeneralizedDeterminant[m_List] := m // GaussElimination // ProductLeadingEntries (* Removing zero rows and zero columns according to Kapur et al.. *) DeleteZeroRows[mat_List] := DeleteCases[mat, 0 Range[Length[mat[[1]]]] ] DeleteZeroColumns[mat_List] := DeleteZeroRows[Transpose[mat]] // Transpose DixonReduction[polys_List, xlist_List, alist_List] := Block[ {dixonmatrix = DixonMatrix[polys, xlist, alist]}, If[ dixonmatrix === 0,0, If[ Head[dixonmatrix] === String, dixonmatrix, dixonmatrix // GaussElimination ] ] ] DixonResultant[polys_List, xlist_List, alist_List] := Block[ {dixonmatrix = DixonMatrix[polys, xlist, alist]}, If[ dixonmatrix === 0,0, If[ Head[dixonmatrix] === String, dixonmatrix, dixonmatrix // DeleteZeroRows // DeleteZeroColumns // GaussElimination // DeleteZeroRows // DeleteZeroColumns // ProductLeadingEntries ] ] ] PreconditionQ[polys_List, xlist_List, alist_List] := Block[ {i, dm = DixonMatrix[polys, xlist, alist], parlist}, If[ dm === 0,Print["\n The Dixon Matrix is zero."], If[ Head[dm] === String, dm, parlist = Variables[dm]; dm = dm /. (Table[parlist[[i]] -> Random[Integer,{1,100}],{i,Length[parlist]}]); dm = RowReduce[dm]; dm = Drop[dm[[1]],1]; If [ ZeroVectorQ[dm], Print["\n The precondition is probabilistically TRUE."], Print["\n The precondition is FALSE."] ] ] ] ] (* Protect all user-accessible functions. *) Protect[ DixonSub, DixonPolynomial, DixonMatrix, ClassicalDixonResultant, DixonReduction, DixonResultant, PreconditionQ, GaussElimination, ProductLeadingEntries, GeneralizedDeterminant ]; End[] ; EndPackage[] ; Print["\n\nPackage Dixon.m, Version 1.0, September 1995\n\n\n User-accessible functions:\n\n { DixonSub, DixonPolynomial, DixonMatrix, ClassicalDixonResultant,\n DixonReduction, DixonResultant, PreconditionQ, GaussElimination,\n ProductLeadingEntries, GeneralizedDeterminant }\n\n "]