/usr/share/axiom-20170501/src/algebra/SYSSOLP.spad is in axiom-source 20170501-3.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 | )abbrev package SYSSOLP SystemSolvePackage
++ Author: P. Gianni
++ Date Created: summer 1988
++ Date Last Updated: summer 1990
++ 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.
SystemSolvePackage(R) : SIG == CODE where
R : IntegralDomain
NNI ==> NonNegativeInteger
P ==> Polynomial
EQ ==> Equation
L ==> List
V ==> Vector
M ==> Matrix
UP ==> SparseUnivariatePolynomial
SE ==> Symbol
IE ==> IndexedExponents Symbol
SUP ==> SparseUnivariatePolynomial
F ==> Fraction Polynomial R
PP2 ==> PolynomialFunctions2(R,F)
PPR ==> Polynomial Polynomial R
SIG ==> 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.
CODE ==> 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
for x 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"/[degree(f,vv)>0 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|(df:=denom f)^=1]
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/[(redPol(fq,pr pretend List(dmp))$gb) ^=0
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 [[]]
|