\documentclass{article} \usepackage{open-axiom} \begin{document} \title{\$SPAD/src/algebra syssolp.spad} \author{Patrizia Gianni} \maketitle \begin{abstract} \end{abstract} \eject \tableofcontents \eject \section{package SYSSOLP SystemSolvePackage} <>= )abbrev package SYSSOLP SystemSolvePackage ++ Author: P. Gianni ++ Date Created: summer 1988 ++ Date Last Updated: summer 1990 ++ Basic Functions: ++ Related Constructors: Fraction, Polynomial, FloatingRealPackage, ++ FloatingComplexPackage, RadicalSolvePackage ++ Also See: LinearSystemMatrixPackage, GroebnerSolve ++ AMS Classifications: ++ Keywords: ++ References: ++ Description: ++ Symbolic solver for systems of rational functions with coefficients ++ in an integral domain R. ++ The systems are solved in the field of rational functions over R. ++ Solutions are exact of the form variable = value when the value is ++ a member of the coefficient domain R. Otherwise the solutions ++ are implicitly expressed as roots of univariate polynomial equations over R. ++ Care is taken to guarantee that the denominators of the input ++ equations do not vanish on the solution sets. ++ The arguments to solve can either be given as equations or ++ as rational functions interpreted as equal ++ to zero. The user can specify an explicit list of symbols to ++ be solved for, treating all other symbols appearing as parameters ++ or omit the list of symbols in which case the system tries to ++ solve with respect to all symbols appearing in the input. NNI ==> NonNegativeInteger P ==> Polynomial EQ ==> Equation L ==> List V ==> Vector M ==> Matrix UP ==> SparseUnivariatePolynomial SE ==> Symbol IE ==> IndexedExponents Symbol SUP ==> SparseUnivariatePolynomial SystemSolvePackage(R): Cat == Cap where R : IntegralDomain F ==> Fraction Polynomial R PP2 ==> PolynomialFunctions2(R,F) PPR ==> Polynomial Polynomial R Cat == with solve: (L F, L SE) -> L L EQ F ++ solve(lp,lv) finds the solutions of the list lp of ++ rational functions with respect to the list of symbols lv. solve: (L EQ F, L SE) -> L L EQ F ++ solve(le,lv) finds the solutions of the ++ list le of equations of rational functions ++ with respect to the list of symbols lv. solve: L F -> L L EQ F ++ solve(lp) finds the solutions of the list lp of rational ++ functions with respect to all symbols appearing in lp. solve: L EQ F -> L L EQ F ++ solve(le) finds the solutions of the list le of equations of ++ rational functions with respect to all symbols appearing in le. solve: (F, SE) -> L EQ F ++ solve(p,v) solves the equation p=0, where p is a rational function ++ with respect to the variable v. solve: (EQ F,SE) -> L EQ F ++ solve(eq,v) finds the solutions of the equation ++ eq with respect to the variable v. solve: F -> L EQ F ++ solve(p) finds the solution of a rational function p = 0 ++ with respect to the unique variable appearing in p. solve: EQ F -> L EQ F ++ solve(eq) finds the solutions of the equation eq ++ with respect to the unique variable appearing in eq. triangularSystems: (L F, L SE) -> L L P R ++ triangularSystems(lf,lv) solves the system of equations ++ defined by lf with respect to the list of symbols lv; ++ the system of equations is obtaining ++ by equating to zero the list of rational functions lf. ++ The output is a list of solutions where ++ each solution is expressed as a "reduced" triangular system of ++ polynomials. Cap == add import MPolyCatRationalFunctionFactorizer(IE,SE,R,P F) ---- Local Functions ---- linSolve: (L F, L SE) -> Union(L EQ F, "failed") makePolys : L EQ F -> L F makeR2F(r : R) : F == r :: (P R) :: F makeP2F(p:P F):F == lv:=variables p lv = [] => retract p for v in lv repeat p:=pushdown(p,v) retract p ---- Local Functions ---- makeEq(p:P F,lv:L SE): EQ F == z:=last lv np:=numer makeP2F p lx:=variables np x : SE for x: free in lv repeat if member?(x,lx) then leave x up:=univariate(np,x) (degree up)=1 => equation(x::P(R)::F,-coefficient(up,0)/leadingCoefficient up) equation(np::F,0$F) varInF(v: SE): F == v::P(R) :: F newInF(n: Integer):F==varInF new()$SE testDegree(f :P R , lv :L SE) : Boolean == "or"/[positive? degree(f,vv) for vv in lv] ---- Exported Functions ---- -- solve a system of rational functions triangularSystems(lf: L F,lv:L SE) : L L P R == empty? lv => empty() empty? lf => empty() #lf = 1 => p:= numer(first lf) fp:=(factor p)$GeneralizedMultivariateFactorize(SE,IE,R,R,P R) [[ff.factor] for ff in factors fp | testDegree(ff.factor,lv)] dmp:=DistributedMultivariatePolynomial(lv,P R) OV:=OrderedVariableList(lv) DP:=DirectProduct(#lv, NonNegativeInteger) push:=PushVariables(R,DP,OV,dmp) lq : L dmp lvv:L OV:=[variable(vv)::OV for vv in lv] lq:=[pushup(df::dmp,lvv)$push for f in lf|not one?(df:=denom f)] lp:=[pushup(numer(f)::dmp,lvv)$push for f in lf] parRes:=groebSolve(lp,lvv)$GroebnerSolve(lv,P R,R) if lq~=[] then gb:=GroebnerInternalPackage(P R,DirectProduct(#lv,NNI),OV,dmp) parRes:=[pr for pr in parRes| and/[not zero?(redPol(fq,pr pretend List(dmp))$gb) for fq in lq]] [[retract pushdown(pf,lvv)$push for pf in pr] for pr in parRes] -- One polynomial. Implicit variable -- solve(pol : F) == zero? pol => error "equation is always satisfied" lv:=removeDuplicates concat(variables numer pol, variables denom pol) empty? lv => error "inconsistent equation" #lv>1 => error "too many variables" solve(pol,first lv) -- general solver. Input in equation style. Implicit variables -- solve(eq : EQ F) == pol:= lhs eq - rhs eq zero? pol => error "equation is always satisfied" lv:=removeDuplicates concat(variables numer pol, variables denom pol) empty? lv => error "inconsistent equation" #lv>1 => error "too many variables" solve(pol,first lv) -- general solver. Input in equation style -- solve(eq:EQ F,var:SE) == solve(lhs eq - rhs eq,var) -- general solver. Input in polynomial style -- solve(pol:F,var:SE) == if R has GcdDomain then p:=primitivePart(numer pol,var) fp:=(factor p)$GeneralizedMultivariateFactorize(SE,IE,R,R,P R) [makeEq(map(makeR2F,ff.factor)$PP2,[var]) for ff in factors fp] else empty() -- Convert a list of Equations in a list of Polynomials makePolys(l: L EQ F):L F == [lhs e - rhs e for e in l] -- linear systems solver. Input as list of polynomials -- linSolve(lp:L F,lv:L SE) == rec:Record(particular:Union(V F,"failed"),basis:L V F) lr : L P R:=[numer f for f in lp] rec:=linSolve(lr,lv)$LinearSystemPolynomialPackage(R,IE,SE,P R) rec.particular case "failed" => "failed" rhs := rec.particular :: V F zeron:V F:=zero(#lv) for p in rec.basis | p ~= zeron repeat sym := newInF(1) for i in 1..#lv repeat rhs.i := rhs.i + sym*p.i eqs: L EQ F := [] for i in 1..#lv repeat eqs := append(eqs,[(lv.i)::(P R)::F = rhs.i]) eqs -- general solver. Input in polynomial style. Implicit variables -- solve(lr : L F) == lv :="setUnion"/[setUnion(variables numer p, variables denom p) for p in lr] solve(lr,lv) -- general solver. Input in equation style. Implicit variables -- solve(le : L EQ F) == lr:=makePolys le lv :="setUnion"/[setUnion(variables numer p, variables denom p) for p in lr] solve(lr,lv) -- general solver. Input in equation style -- solve(le:L EQ F,lv:L SE) == solve(makePolys le, lv) checkLinear(lr:L F,vl:L SE):Boolean == ld:=[denom pol for pol in lr] for f in ld repeat if (or/[member?(x,vl) for x in variables f]) then return false and/[totalDegree(numer pol,vl) < 2 for pol in lr] -- general solver. Input in polynomial style -- solve(lr:L F,vl:L SE) == empty? vl => empty() checkLinear(lr,vl) => -- linear system -- soln := linSolve(lr, vl) soln case "failed" => [] eqns: L EQ F := [] for i in 1..#vl repeat lhs := (vl.i::(P R))::F rhs := rhs soln.i eqns := append(eqns, [lhs = rhs]) [eqns] -- polynomial system -- if R has GcdDomain then parRes:=triangularSystems(lr,vl) [[makeEq(map(makeR2F,f)$PP2,vl) for f in pr] for pr in parRes] else [[]] @ \section{License} <>= --Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd. --All rights reserved. -- --Redistribution and use in source and binary forms, with or without --modification, are permitted provided that the following conditions are --met: -- -- - Redistributions of source code must retain the above copyright -- notice, this list of conditions and the following disclaimer. -- -- - Redistributions in binary form must reproduce the above copyright -- notice, this list of conditions and the following disclaimer in -- the documentation and/or other materials provided with the -- distribution. -- -- - Neither the name of The Numerical ALgorithms Group Ltd. nor the -- names of its contributors may be used to endorse or promote products -- derived from this software without specific prior written permission. -- --THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS --IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED --TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A --PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER --OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, --EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, --PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR --PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF --LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING --NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS --SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. @ <<*>>= <> <> @ \eject \begin{thebibliography}{99} \bibitem{1} nothing \end{thebibliography} \end{document}