|
This is the second part of the source from [mantepse.spad.pamphlet]
\begin{axiom}
)lib RECOP FAMR2 FFFG FFFGF NEWTON GOPT GOPT0 UFPS UFPS1
\end{axiom}
\begin{spad}
)abbrev package GUESS Guess
++ Author: Martin Rubey
++ Description: This package implements guessing of sequences. Packages for the
++ most common cases are provided as \spadtype{GuessInteger},
++ \spadtype{GuessPolynomial}, etc.
Guess(F, S, EXPRR, R, retract, coerce): Exports == Implementation where
F: Field -- zB.: FRAC POLY PF 5
-- in F wird interpoliert und gerechnet S: GcdDomain
- -- in guessExpRat möchte ich die Wurzeln von Polynomen über F bestimmen. Wenn F
- ein Quotientenkörper ist, dann kann ich den Nenner loswerden.
-- F wäre dann also QFCAT S
- R: Join(OrderedSet, IntegralDomain)
- zB.: FRAC POLY INT
- -- die Ergebnisse werden aber in EXPRR dargestellt
- EXPRR: Join(ExpressionSpace, IntegralDomain,
EXPRR: Join(FunctionSpace Integer, IntegralDomain,
RetractableTo R, RetractableTo Symbol,
RetractableTo Integer, CombinatorialOpsCategory,
PartialDifferentialRing Symbol) with
_ : (%,%) -> %
_/ : (%,%) -> %
__* : (%,%) -> %
numerator : % -> %
denominator : % -> %
ground? : % -> Boolean
-- zB.: EXPR INT
-- EXPR FRAC POLY INT ist ja verboten. Deshalb kann ich nicht einfach EXPR R
-- verwenden
-- EXPRR existiert, falls irgendwann einmal EXPR PF 5 unterstützt wird.
-- das folgende möchte ich eigentlich loswerden.
- retract: R -> F
- zB.: i+->i
coerce: F -> EXPRR -- zB.: i+->i
-- Achtung EXPRR ~= EXPR R
LGOPT ==> List GuessOption
GOPT0 ==> GuessOptionFunctions0
NNI ==> NonNegativeInteger
PI ==> PositiveInteger
EXPRI ==> Expression Integer
GUESSRESULT ==> List Record(function: EXPRR, order: NNI)
UFPSF ==> UnivariateFormalPowerSeries F
UFPS1 ==> UnivariateFormalPowerSeriesFunctions
UFPSS ==> UnivariateFormalPowerSeries S
SUP ==> SparseUnivariatePolynomial
UFPSSUPF ==> UnivariateFormalPowerSeries SUP F
FFFG ==> FractionFreeFastGaussian
FFFGF ==> FractionFreeFastGaussianFractions
-- CoeffAction
DIFFSPECA ==> (NNI, NNI, SUP S) -> S DIFFSPECAF ==> (NNI, NNI, UFPSSUPF) -> SUP F
DIFFSPECAX ==> (NNI, Symbol, EXPRR) -> EXPRR
-- the diagonal of the C-matrix
DIFFSPECC ==> NNI -> List S HPSPEC ==> Record(guessStream: UFPSF -> Stream UFPSF,
degreeStream: Stream NNI,
testStream: UFPSSUPF -> Stream UFPSSUPF,
exprStream: (EXPRR, Symbol) -> Stream EXPRR,
A: DIFFSPECA,
AF: DIFFSPECAF,
AX: DIFFSPECAX,
C: DIFFSPECC)
- -- note that empty?(guessStream.o) has to return always. In other words, if the
- stream is finite, empty? should recognize it.
- DIFFSPECN ==> EXPRR -> EXPRR
- eg.: i+->q^i
GUESSER ==> (List F, LGOPT) -> GUESSRESULT
FSUPS ==> Fraction SUP S
FSUPF ==> Fraction SUP F
V ==> OrderedVariableList(['a1, 'A])
POLYF ==> SparseMultivariatePolynomial(F, V)
FPOLYF ==> Fraction POLYF
FSUPFPOLYF ==> Fraction SUP FPOLYF
POLYS ==> SparseMultivariatePolynomial(S, V)
FPOLYS ==> Fraction POLYS
FSUPFPOLYS ==> Fraction SUP FPOLYS
--<>
Exports == with guess: List F -> GUESSRESULT
++ \spad{guess l} applies recursively \spadfun{guessRec} and
++ \spadfun{guessADE} to the successive differences and quotients of
++ the list. Default options as described in
++ \spadtype{GuessOptionFunctions0} are used.
guess: (List F, LGOPT) -> GUESSRESULT
++ \spad{guess(l, options)} applies recursively \spadfun{guessRec}
++ and \spadfun{guessADE} to the successive differences and quotients
++ of the list. The given options are used.
guess: (List F, List GUESSER, List Symbol) -> GUESSRESULT
++ \spad{guess(l, guessers, ops)} applies recursively the given
++ guessers to the successive differences if ops contains the symbol
++ guessSum and quotients if ops contains the symbol guessProduct to
++ the list. Default options as described in
++ \spadtype{GuessOptionFunctions0} are used.
guess: (List F, List GUESSER, List Symbol, LGOPT) -> GUESSRESULT
++ \spad{guess(l, guessers, ops)} applies recursively the given
++ guessers to the successive differences if ops contains the symbol
++ \spad{guessSum} and quotients if ops contains the symbol
++ \spad{guessProduct} to the list. The given options are used.
guessExpRat: List F -> GUESSRESULT
++ \spad{guessExpRat l} tries to find a function of the form
++ n+->(a+b n)^n r(n), where r(n) is a rational function, that fits
++ l.
guessExpRat: (List F, LGOPT) -> GUESSRESULT
++ \spad{guessExpRat(l, options)} tries to find a function of the
++ form n+->(a+b n)^n r(n), where r(n) is a rational function, that
++ fits l.
if F has RetractableTo Symbol and S has RetractableTo Symbol then
guessExpRat: Symbol -> GUESSER
++ \spad{guessExpRat q} returns a guesser that tries to find a
++ function of the form n+->(a+b q^n)^n r(q^n), where r(n) is a
++ rational function, that fits l.
guessHP: (LGOPT -> HPSPEC) -> GUESSER
++ \spad{guessHP f} constructs an operation that applies Hermite-Pade
++ approximation to the series generated by the given function f.
guessADE: List F -> GUESSRESULT
++ \spad{guessADE l} tries to find an algebraic differential equation
++ for a generating function whose first Taylor coefficients are
++ given by l, using the default options described in
++ \spadtype{GuessOptionFunctions0}.
guessADE: (List F, LGOPT) -> GUESSRESULT
++ \spad{guessADE(l, options)} tries to find an algebraic
++ differential equation for a generating function whose first Taylor
++ coefficients are given by l, using the given options.
guessAlg: List F -> GUESSRESULT
++ \spad{guessAlg l} tries to find an algebraic equation for a
++ generating function whose first Taylor coefficients are given by
++ l, using the default options described in
++ \spadtype{GuessOptionFunctions0}. It is equivalent to
++ \spadfun{guessADE}(l, maxDerivative == 0).
guessAlg: (List F, LGOPT) -> GUESSRESULT
++ \spad{guessAlg(l, options)} tries to find an algebraic equation
++ for a generating function whose first Taylor coefficients are
++ given by l, using the given options. It is equivalent to
++ \spadfun{guessADE}(l, options) with \spad{maxDerivative == 0}.
guessHolo: List F -> GUESSRESULT
++ \spad{guessHolo l} tries to find an ordinary linear differential
++ equation for a generating function whose first Taylor coefficients
++ are given by l, using the default options described in
++ \spadtype{GuessOptionFunctions0}. It is equivalent to
++ \spadfun{guessADE}\spad{(l, maxPower == 1)}.
guessHolo: (List F, LGOPT) -> GUESSRESULT
++ \spad{guessHolo(l, options)} tries to find an ordinary linear
++ differential equation for a generating function whose first Taylor
++ coefficients are given by l, using the given options. It is
++ equivalent to \spadfun{guessADE}\spad{(l, options)} with
++ \spad{maxPower == 1}.
guessPade: (List F, LGOPT) -> GUESSRESULT
++ \spad{guessPade(l, options)} tries to find a rational function
++ whose first Taylor coefficients are given by l, using the given
++ options. It is equivalent to \spadfun{guessADE}\spad{(l,
++ maxDerivative == 0, maxPower == 1, allDegrees == true)}.
guessPade: List F -> GUESSRESULT
++ \spad{guessPade(l, options)} tries to find a rational function
++ whose first Taylor coefficients are given by l, using the default
++ options described in \spadtype{GuessOptionFunctions0}. It is
++ equivalent to \spadfun{guessADE}\spad{(l, options)} with
++ \spad{maxDerivative == 0, maxPower == 1, allDegrees == true}.
guessRec: List F -> GUESSRESULT
++ \spad{guessRec l} tries to find an ordinary difference equation
++ whose first values are given by l, using the default options
++ described in \spadtype{GuessOptionFunctions0}.
guessRec: (List F, LGOPT) -> GUESSRESULT
++ \spad{guessRec(l, options)} tries to find an ordinary difference
++ equation whose first values are given by l, using the given
++ options.
guessPRec: (List F, LGOPT) -> GUESSRESULT
++ \spad{guessPRec(l, options)} tries to find a linear recurrence
++ with polynomial coefficients whose first values are given by l,
++ using the given options. It is equivalent to
++ \spadfun{guessRec}\spad{(l, options)} with \spad{maxPower == 1}.
guessPRec: List F -> GUESSRESULT
++ \spad{guessPRec l} tries to find a linear recurrence with
++ polynomial coefficients whose first values are given by l, using
++ the default options described in
++ \spadtype{GuessOptionFunctions0}. It is equivalent to
++ \spadfun{guessRec}\spad{(l, maxPower == 1)}.
guessRat: (List F, LGOPT) -> GUESSRESULT
++ \spad{guessRat(l, options)} tries to find a rational function
++ whose first values are given by l, using the given options. It is
++ equivalent to \spadfun{guessRec}\spad{(l, maxShift == 0, maxPower
++ == 1, allDegrees == true)}.
guessRat: List F -> GUESSRESULT
++ \spad{guessRat l} tries to find a rational function whose first
++ values are given by l, using the default options described in
++ \spadtype{GuessOptionFunctions0}. It is equivalent to
++ \spadfun{guessRec}\spad{(l, maxShift == 0, maxPower == 1,
++ allDegrees == true)}.
diffHP: LGOPT -> HPSPEC
++ \spad{diffHP options} returns a specification for Hermite-Pade
++ approximation with the differential operator
shiftHP: LGOPT -> HPSPEC
++ \spad{shiftHP options} returns a specification for Hermite-Pade
++ approximation with the shift operator
if F has RetractableTo Symbol and S has RetractableTo Symbol then
shiftHP: Symbol -> (LGOPT -> HPSPEC)
++ \spad{shiftHP options} returns a specification for
++ Hermite-Pade approximation with the $q$-shift operator diffHP: Symbol -> (LGOPT -> HPSPEC)
++ \spad{diffHP options} returns a specification for Hermite-Pade
++ approximation with the $q$-dilation operator
guessRec: Symbol -> GUESSER
++ \spad{guessRec q} returns a guesser that finds an ordinary
++ q-difference equation whose first values are given by l, using
++ the given options.
guessPRec: Symbol -> GUESSER
++ \spad{guessPRec q} returns a guesser that tries to find
++ a linear q-recurrence with polynomial coefficients whose first
++ values are given by l, using the given options. It is
++ equivalent to \spadfun{guessRec}\spad{(q)} with
++ \spad{maxPower == 1}.
guessRat: Symbol -> GUESSER
++ \spad{guessRat q} returns a guesser that tries to find a
++ q-rational function whose first values are given by l, using
++ the given options. It is equivalent to \spadfun{guessRec} with
++ \spad{(l, maxShift == 0, maxPower == 1, allDegrees == true)}.
guessADE: Symbol -> GUESSER
++ \spad{guessADE q} returns a guesser that tries to find an
++ algebraic differential equation for a generating function whose
++ first Taylor coefficients are given by l, using the given
++ options.
--<>
Implementation == add
-- We have to put this chunk at the beginning, because otherwise it will take
-- very long to compile.
ord1(x: List Integer, i: Integer): Integer ==
n := #x - 3 - i
x.(n+1)*reduce(_+, [x.j for j in 1..n], 0) + _
2*reduce(_+, [reduce(_+, [x.k*x.j for k in 1..j-1], 0) _
for j in 1..n], 0)
ord2(x: List Integer, i: Integer): Integer ==
if zero? i then
n := #x - 3 - i
ord1(x, i) + reduce(_+, [x.j for j in 1..n], 0)*(x.(n+2)-x.(n+1))
else
ord1(x, i)
deg1(x: List Integer, i: Integer): Integer ==
m := #x - 3
(x.(m+3)+x.(m+1)+x.(1+i))*reduce(_+, [x.j for j in 2+i..m], 0) + _
x.(m+3)x.(m+1) + _
2reduce(_+, [reduce(_+, [x.k*x.j for k in 2+i..j-1], 0) _
for j in 2+i..m], 0)
deg2(x: List Integer, i: Integer): Integer ==
m := #x - 3
deg1(x, i) + _
(x.(m+3) + reduce(_+, [x.j for j in 2+i..m], 0)) * _
(x.(m+2)-x.(m+1))
checkResult(res: EXPRR, n: Symbol, l: Integer, list: List F,
options: LGOPT): NNI ==
for i in l..1 by -1 repeat
den := eval(denominator res, n::EXPRR, (i-1)::EXPRR)
if den = 0 then return i::NNI
num := eval(numerator res, n::EXPRR, (i-1)::EXPRR)
if list.i ~= retract(retract(num/den)@R)
then return i::NNI
0$NNI
SUPS2SUPF(p: SUP S): SUP F ==
if F is S then
p pretend SUP(F)
else if F is Fraction S then
map(coerce(#1)$Fraction(S), p)
$SparseUnivariatePolynomialFunctions2(S, F)
else error "Type parameter F should be either equal to S or equal _
to Fraction S"
F2FPOLYS(p: F): FPOLYS ==
if F is S then
p::POLYF::FPOLYF pretend FPOLYS
else if F is Fraction S then
numer(p)$Fraction(S)::POLYS/denom(p)$Fraction(S)::POLYS
else error "Type parameter F should be either equal to S or equal _
to Fraction S"
MPCSF ==> MPolyCatFunctions2(V, IndexedExponents V,
IndexedExponents V, S, F,
POLYS, POLYF)
SUPF2EXPRR(xx: Symbol, p: SUP F): EXPRR ==
zero? p => 0
(coerce(leadingCoefficient p))::EXPRR (xx::EXPRR)*degree p
+ SUPF2EXPRR(xx, reductum p)
FSUPF2EXPRR(xx: Symbol, p: FSUPF): EXPRR ==
(SUPF2EXPRR(xx, numer p)) / (SUPF2EXPRR(xx, denom p))
POLYS2POLYF(p: POLYS): POLYF ==
if F is S then
p pretend POLYF
else if F is Fraction S then
map(coerce(#1)$Fraction(S), p)$MPCSF
else error "Type parameter F should be either equal to S or equal _
to Fraction S"
SUPPOLYS2SUPF(p: SUP POLYS, a1v: F, Av: F): SUP F ==
zero? p => 0
lc: POLYF := POLYS2POLYF leadingCoefficient(p)
monomial(retract(eval(lc, [index(1)$V, index(2)$V]::List V,
[a1v, Av])),
degree p)
+ SUPPOLYS2SUPF(reductum p, a1v, Av)
SUPFPOLYS2FSUPPOLYS(p: SUP FPOLYS): Fraction SUP POLYS ==
cden := splitDenominator(p)
$UnivariatePolynomialCommonDenominator(POLYS, FPOLYS,SUP FPOLYS) pnum: SUP POLYS
:= map(retract(#1 * cden.den)$FPOLYS, p)
$SparseUnivariatePolynomialFunctions2(FPOLYS, POLYS)
pden: SUP POLYS := (cden.den)::SUP POLYS
pnum/pden
POLYF2EXPRR(p: POLYF): EXPRR ==
map(convert(#1)@Symbol::EXPRR, coerce(#1)@EXPRR, p)
$PolynomialCategoryLifting(IndexedExponents V, V,
F, POLYF, EXPRR)
defaultD: DIFFSPECN
defaultD(expr: EXPRR): EXPRR == expr
-- applies n+->q^n or whatever DN is to i
DN2DL: (DIFFSPECN, Integer) -> F
DN2DL(DN, i) == retract(retract(DN(i::EXPRR))@R)
evalResultant(p1: POLYS, p2: POLYS, o: Integer, d: Integer, va1: V, vA: V)_
: List S ==
res: List S := []
d1 := degree(p1, va1)
d2 := degree(p2, va1)
lead: S
for k in 1..d-o+1 repeat
p1atk := univariate eval(p1, vA, k::S)
p2atk := univariate eval(p2, vA, k::S) d1atk := degree p1atk
d2atk := degree p2atk
-- output("k: " string(k))$OutputPackage
-- output("d1: " string(d1) " d1atk: " string(d1atk))$OutputPackage
-- output("d2: " string(d2) " d2atk: " string(d2atk))$OutputPackage if d2atk < d2 then
if d1atk < d1
then lead := 0$S
else lead := (leadingCoefficient p1atk)((d2-d2atk)::NNI)
else
if d1atk < d1
then lead := (-1$S)d2 (leadingCoefficient p2atk)*((d1-d1atk)::NNI)
else lead := 1$S
if zero? lead
then res := cons(0, res)
else res := cons(lead (resultant(p1atk, p2atk)$SUP(S) exquo _
(k::S)*(o::NNI))::S,
res)
reverse res
p(xm: Integer, i: Integer, va1: V, vA: V, basis: DIFFSPECN): FPOLYS ==
vA::POLYS::FPOLYS + va1::POLYS::FPOLYS _
* F2FPOLYS(DN2DL(basis, i) - DN2DL(basis, xm))
p2(xm: Integer, i: Symbol, a1v: F, Av: F, basis: DIFFSPECN): EXPRR ==
coerce(Av) + coerce(a1v)*(basis(i::EXPRR) - basis(xm::EXPRR))
GF ==> GeneralizedMultivariateFactorize(SingletonAsOrderedSet,
IndexedExponents V, F, F,
SUP F)
guessExpRatAux(xx: Symbol, list: List F, basis: DIFFSPECN,
xValues: List Integer, options: LGOPT): List EXPRR == a1: V := index(1)$V
A: V := index(2)$V
len: NNI := #list
if len < 4 then return []
else len := (len-3)::NNI
xlist := [F2FPOLYS DN2DL(basis, xValues.i) for i in 1..len]
x1 := F2FPOLYS DN2DL(basis, xValues.(len+1))
x2 := F2FPOLYS DN2DL(basis, xValues.(len+2))
x3 := F2FPOLYS DN2DL(basis, xValues.(len+3))
y: NNI -> FPOLYS :=
F2FPOLYS(list.#1) _
p(last xValues, (xValues.#1)::Integer, a1, A, basis)*_
(-(xValues.#1)::Integer)
ylist: List FPOLYS := [y i for i in 1..len]
y1 := y(len+1)
y2 := y(len+2)
y3 := y(len+3)
res := []::List EXPRR
if maxDegree(options)$GOPT0 = -1
then maxDeg := len-1
else maxDeg := min(maxDegree(options)$GOPT0, len-1)
for i in 0..maxDeg repeat
if debug(options)$GOPT0 then
output(hconcat(degree ExpRat , i::OutputForm))
$OutputPackage if debug(options)$GOPT0 then
systemCommand("sys date +%s")$MoreSystemCommands
output(interpolating...)$OutputPackage
ri: FSUPFPOLYS
:= interpolate(xlist, ylist, (len-1-i)::NNI) _
$FFFG(FPOLYS, SUP FPOLYS)
-- for experimental fraction free interpolation
-- ri: Fraction SUP POLYS
-- := interpolate(xlist, ylist, (len-1-i)::NNI) _
-- $FFFG(POLYS, SUP POLYS)
- if debug(options)$GOPT0 then
- output(hconcat("xlist: ", xlist::OutputForm))$OutputPackage
-- output(hconcat("ylist: ", ylist::OutputForm))$OutputPackage
-- output(hconcat("ri: ", ri::OutputForm))$OutputPackage
systemCommand("sys date +%s")$MoreSystemCommands
output(polynomials...)$OutputPackage
poly1: POLYS := numer(elt(ri, x1)$SUP(FPOLYS) - y1)
poly2: POLYS := numer(elt(ri, x2)$SUP(FPOLYS) - y2)
poly3: POLYS := numer(elt(ri, x3)$SUP(FPOLYS) - y3)
-- for experimental fraction free interpolation
-- ri2: FSUPFPOLYS := map(#1::FPOLYS, numer ri)
-- $SparseUnivariatePolynomialFunctions2(POLYS, FPOLYS)
-- /map(#1::FPOLYS, denom ri) _
-- $SparseUnivariatePolynomialFunctions2(POLYS, FPOLYS)
--
-- poly1: POLYS := numer(elt(ri2, x1)$SUP(FPOLYS) - y1)
-- poly2: POLYS := numer(elt(ri2, x2)$SUP(FPOLYS) - y2)
-- poly3: POLYS := numer(elt(ri2, x3)$SUP(FPOLYS) - y3) n:Integer := len - i
o1: Integer := ord1(xValues, i)
d1: Integer := deg1(xValues, i)
o2: Integer := ord2(xValues, i)
d2: Integer := deg2(xValues, i)
-- another compiler bug: using i as iterator here makes the loop break
if debug(options)$GOPT0 then
systemCommand("sys date +%s")$MoreSystemCommands
output(interpolating resultants...)$OutputPackage
res1: SUP S
:= newton(evalResultant(poly1, poly3, o1, d1, a1, A))
$NewtonInterpolation(S)
res2: SUP S
:= newton(evalResultant(poly2, poly3, o2, d2, a1, A))
$NewtonInterpolation(S)
- if debug(options)$GOPT0 then
- res1: SUP S := univariate(resultant(poly1, poly3, a1))
-- res2: SUP S := univariate(resultant(poly2, poly3, a1))
-- if res1 ~= res1res or res2 ~= res2res then
-- output(hconcat("poly1 ", poly1::OutputForm))$OutputPackage
-- output(hconcat("poly2 ", poly2::OutputForm))$OutputPackage
-- output(hconcat("poly3 ", poly3::OutputForm))$OutputPackage
-- output(hconcat("res1 ", res1::OutputForm))$OutputPackage
-- output(hconcat("res2 ", res2::OutputForm))$OutputPackage
output("n/i: " string(n) " " string(i))$OutputPackage
output("res1 ord: " string(o1) " " string(minimumDegree res1))
$OutputPackage
output("res1 deg: " string(d1) " " string(degree res1))
$OutputPackage
output("res2 ord: " string(o2) " " string(minimumDegree res2))
$OutputPackage
output("res2 deg: " string(d2) " " string(degree res2))
$OutputPackage
if debug(options)$GOPT0 then
systemCommand("sys date +%s")$MoreSystemCommands
output(computing gcd...)$OutputPackage
-- we want to solve over F
res3: SUP F := SUPS2SUPF(primitivePart(gcd(res1, res2))) if debug(options)$GOPT0 then
systemCommand("sys date +%s")$MoreSystemCommands
output(solving...)$OutputPackage
-- res3 is a polynomial in A=a0+(len+3)*a1
-- now we have to find the roots of res3
- for f in factors factor(res3)$GF | degree f.factor = 1 repeat
- we are only interested in the linear factors
if debug(options)$GOPT0 then
output(hconcat("f: ", f::OutputForm))$OutputPackage
Av: F := -coefficient(f.factor, 0)
/ leadingCoefficient f.factor
-- FIXME: in an earlier version, we disregarded vanishing Av
-- maybe we intended to disregard vanishing a1v? Either doesn't really
-- make sense to me right now. evalPoly := eval(POLYS2POLYF poly3, A, Av)
if zero? evalPoly
then evalPoly := eval(POLYS2POLYF poly1, A, Av)
-- Note that it really may happen that poly3 vanishes when specializing
-- A. Consider for example guessExpRat([1,1,1,1]).
-- FIXME: We check poly1 below, too. I should work out in what cases poly3
-- vanishes. for g in factors factor(univariate evalPoly)$GF
| degree g.factor = 1 repeat
if debug(options)$GOPT0 then
output(hconcat("g: ", g::OutputForm))$OutputPackage a1v: F := -coefficient(g.factor, 0)
/ leadingCoefficient g.factor
-- check whether poly1 and poly2 really vanish. Note that we could have found
-- an extraneous solution, since we only computed the gcd of the two
-- resultants. t1 := eval(POLYS2POLYF poly1, [a1, A]::List V,
[a1v, Av]::List F)
if zero? t1 then
t2 := eval(POLYS2POLYF poly2, [a1, A]::List V,
[a1v, Av]::List F)
if zero? t2 then ri1: Fraction SUP POLYS
:= SUPFPOLYS2FSUPPOLYS(numer ri)
/ SUPFPOLYS2FSUPPOLYS(denom ri)
-- for experimental fraction free interpolation
-- ri1: Fraction SUP POLYS := ri numr: SUP F := SUPPOLYS2SUPF(numer ri1, a1v, Av)
denr: SUP F := SUPPOLYS2SUPF(denom ri1, a1v, Av)
if not zero? denr then
res4: EXPRR := eval(FSUPF2EXPRR(xx, numr/denr),
kernel(xx),
basis(xx::EXPRR))
p2(last xValues,
xx, a1v, Av, basis)
*xx::EXPRR
res := cons(res4, res)
else if zero? numr and debug(options)$GOPT0 then
output("numerator and denominator vanish!")
$OutputPackage
-- If we are only interested in one solution, we do not try other degrees if we
-- have found already some solutions. I.e., the indentation here is correct. if not null(res) and one(options)$GOPT0 then return res
res
guessExpRatAux0(list: List F, basis: DIFFSPECN, options: LGOPT): GUESSRESULT ==
if zero? safety(options)$GOPT0 then
error "Guess: guessExpRat does not support zero safety"
-- guesses Functions of the Form (a1n+a0)^nrat(n)
xx := indexName(options)$GOPT0
-- restrict to safety
len: Integer := #list
if len-safety(options)$GOPT0+1 < 0 then return []
shortlist: List F := first(list, (len-safety(options)$GOPT0+1)::NNI)
-- remove zeros from list
zeros: EXPRR := 1
newlist: List F
xValues: List Integer
i: Integer := -1
for x in shortlist repeat
i := i+1
if x = 0 then
zeros := zeros * (basis(xx::EXPRR) - basis(i::EXPRR))
i := -1
for x in shortlist repeat
i := i+1
if x ~= 0 then
newlist := cons(x/retract(retract(eval(zeros, xx::EXPRR,
i::EXPRR))@R),
newlist)
xValues := cons(i, xValues)
newlist := reverse newlist
xValues := reverse xValues
res: List EXPRR
:= [eval(zeros * f, xx::EXPRR, xx::EXPRR) _
for f in guessExpRatAux(xx, newlist, basis, xValues, options)]
reslist := map([#1, checkResult(#1, xx, len, list, options)], res)
$ListFunctions2(EXPRR, Record(function: EXPRR, order: NNI))
select(#1.order < len-safety(options)$GOPT0, reslist)
guessExpRat(list : List F): GUESSRESULT ==
guessExpRatAux0(list, defaultD, [])
guessExpRat(list: List F, options: LGOPT): GUESSRESULT ==
guessExpRatAux0(list, defaultD, options)
if F has RetractableTo Symbol and S has RetractableTo Symbol then
guessExpRat(q: Symbol): GUESSER ==
guessExpRatAux0(#1, q::EXPRR**#1, #2)
-- some useful types for Ore operators that work on series
-- the differentiation operator
DIFFSPECX ==> (EXPRR, Symbol, NonNegativeInteger) -> EXPRR
-- eg.: f(x)+->f(qx)
-- f(x)+->D(f, x)
DIFFSPECS ==> (UFPSF, NonNegativeInteger) -> UFPSF
-- eg.: f(x)+->f(qx)
- DIFFSPECSF ==> (UFPSSUPF, NonNegativeInteger) -> UFPSSUPF
- eg.: f(x)+->f(q*x)
-- the constant term for the inhomogeneous case
DIFFSPEC1 ==> UFPSF
DIFFSPEC1F ==> UFPSSUPF
DIFFSPEC1X ==> Symbol -> EXPRR
termAsEXPRR(f: EXPRR, xx: Symbol, l: List Integer,
DX: DIFFSPECX, D1X: DIFFSPEC1X): EXPRR ==
if empty? l then D1X(xx)
else
ll: List List Integer := powers(l)$Partition fl: List EXPRR := [DX(f, xx, (first part-1)::NonNegativeInteger)
* second(part)::NNI for part in ll]
reduce(_, fl)
termAsUFPSF(f: UFPSF, l: List Integer, DS: DIFFSPECS, D1: DIFFSPEC1): UFPSF ==
if empty? l then D1
else
ll: List List Integer := powers(l)$Partition
-- first of each element of ll is the derivative, second is the power
fl: List UFPSF := [DS(f, (first part -1)::NonNegativeInteger) _
** second(part)::NNI for part in ll]
reduce(_*, fl)
-- returns \prod f^(l.i), but using the Hadamard product
termAsUFPSF2(f: UFPSF, l: List Integer,
DS: DIFFSPECS, D1: DIFFSPEC1): UFPSF ==
if empty? l then D1
else
ll: List List Integer := powers(l)$Partition
-- first of each element of ll is the derivative, second is the power
fl: List UFPSF
:= [map(#1** second(part)::NNI, DS(f, (first part -1)::NNI)) _
for part in ll]
reduce(hadamard$UFPS1(F), fl)
termAsUFPSSUPF(f: UFPSSUPF, l: List Integer,
DSF: DIFFSPECSF, D1F: DIFFSPEC1F): UFPSSUPF ==
if empty? l then D1F
else
ll: List List Integer := powers(l)$Partition
-- first of each element of ll is the derivative, second is the power
fl: List UFPSSUPF
:= [DSF(f, (first part -1)::NonNegativeInteger)
** second(part)::NNI for part in ll]
reduce(_*, fl)
-- returns \prod f^(l.i), but using the Hadamard product
termAsUFPSSUPF2(f: UFPSSUPF, l: List Integer,
DSF: DIFFSPECSF, D1F: DIFFSPEC1F): UFPSSUPF ==
if empty? l then D1F
else
ll: List List Integer := powers(l)$Partition
-- first of each element of ll is the derivative, second is the power
fl: List UFPSSUPF
:= [map(#1 ** second(part)::NNI, DSF(f, (first part -1)::NNI)) _
for part in ll]
reduce(hadamard$UFPS1(SUP F), fl)
FilteredPartitionStream(options: LGOPT): Stream List Integer ==
maxD := 1+maxDerivative(options)$GOPT0
maxP := maxPower(options)$GOPT0 if maxD > 0 and maxP > -1 then
s := partitions(maxD, maxP)$PartitionsAndPermutations
else
s1: Stream Integer := generate(inc, 1)$Stream(Integer)
s2: Stream Stream List Integer
:= map(partitions(#1)$PartitionsAndPermutations, s1)
$StreamFunctions2(Integer, Stream List Integer)
s3: Stream List Integer
:= concat(s2)$StreamFunctions1(List Integer)
-- s := cons([],
-- select(((maxD = 0) or (first #1 <= maxD)) _
-- and ((maxP = -1) or (# #1 <= maxP)), s3)) s := cons([],
select(((maxD = 0) or (# #1 <= maxD)) _
and ((maxP = -1) or (first #1 <= maxP)), s3))
s := conjugates(s)$PartitionsAndPermutations
if homogeneous(options)$GOPT0 then rest s else s -- for functions
ADEguessStream(f: UFPSF, partitions: Stream List Integer,
DS: DIFFSPECS, D1: DIFFSPEC1): Stream UFPSF ==
map(termAsUFPSF(f, #1, DS, D1), partitions)
$StreamFunctions2(List Integer, UFPSF)
-- for coefficients, i.e., using the Hadamard product
ADEguessStream2(f: UFPSF, partitions: Stream List Integer,
DS: DIFFSPECS, D1: DIFFSPEC1): Stream UFPSF ==
map(termAsUFPSF2(f, #1, DS, D1), partitions)
$StreamFunctions2(List Integer, UFPSF)
ADEdegreeStream(partitions: Stream List Integer): Stream NNI ==
scan(0, max((if empty? #1 then 0 else (first #1 - 1)::NNI), #2),
partitions)$StreamFunctions2(List Integer, NNI)
ADEtestStream(f: UFPSSUPF, partitions: Stream List Integer,
DSF: DIFFSPECSF, D1F: DIFFSPEC1F): Stream UFPSSUPF ==
map(termAsUFPSSUPF(f, #1, DSF, D1F), partitions)
$StreamFunctions2(List Integer, UFPSSUPF)
ADEtestStream2(f: UFPSSUPF, partitions: Stream List Integer,
DSF: DIFFSPECSF, D1F: DIFFSPEC1F): Stream UFPSSUPF ==
map(termAsUFPSSUPF2(f, #1, DSF, D1F), partitions)
$StreamFunctions2(List Integer, UFPSSUPF)
ADEEXPRRStream(f: EXPRR, xx: Symbol, partitions: Stream List Integer,
DX: DIFFSPECX, D1X: DIFFSPEC1X): Stream EXPRR ==
map(termAsEXPRR(f, xx, #1, DX, D1X), partitions)
$StreamFunctions2(List Integer, EXPRR)
diffDX: DIFFSPECX
diffDX(expr, x, n) == D(expr, x, n)
diffDS: DIFFSPECS
diffDS(s, n) == D(s, n)
diffDSF: DIFFSPECSF
diffDSF(s, n) ==
-- I have to help the compiler here a little to choose the right signature...
if SUP F has _*: (NonNegativeInteger, SUP F) -> SUP F
then D(s, n)
diffAX: DIFFSPECAX
diffAX(l: NNI, x: Symbol, f: EXPRR): EXPRR ==
(x::EXPRR)*l f
diffA: DIFFSPECA
diffA(k: NNI, l: NNI, f: SUP S): S ==
DiffAction(k, l, f)$FFFG(S, SUP S)
diffAF: DIFFSPECAF
diffAF(k: NNI, l: NNI, f: UFPSSUPF): SUP F ==
DiffAction(k, l, f)$FFFG(SUP F, UFPSSUPF)
diffC: DIFFSPECC
diffC(total: NNI): List S == DiffC(total)$FFFG(S, SUP S)
diff1X: DIFFSPEC1X
diff1X(x: Symbol)== 1$EXPRR
diffHP options ==
if displayAsGF(options)$GOPT0 then
partitions := FilteredPartitionStream options
[ADEguessStream(#1, partitions, diffDS, 1$UFPSF),
ADEdegreeStream partitions,
ADEtestStream(#1, partitions, diffDSF, 1$UFPSSUPF),
ADEEXPRRStream(#1, #2, partitions, diffDX, diff1X),
diffA, diffAF, diffAX, diffC]$HPSPEC
else
error "Guess: guessADE supports only displayAsGF"
if F has RetractableTo Symbol and S has RetractableTo Symbol then
qDiffDX(q: Symbol, expr: EXPRR, x: Symbol, n: NonNegativeInteger): EXPRR ==
eval(expr, x::EXPRR, (q::EXPRR)*nx::EXPRR)
qDiffDS(q: Symbol, s: UFPSF, n: NonNegativeInteger): UFPSF ==
multiplyCoefficients((q::F)*((n#1)::NonNegativeInteger), s)
qDiffDSF(q: Symbol, s: UFPSSUPF, n: NonNegativeInteger): UFPSSUPF ==
multiplyCoefficients((q::F::SUP F)*((n#1)::NonNegativeInteger), s)
diffHP(q: Symbol): (LGOPT -> HPSPEC) ==
if displayAsGF(#1)$GOPT0 then
partitions := FilteredPartitionStream #1
[ADEguessStream(#1, partitions, qDiffDS(q, #1, #2), 1$UFPSF), _
repeating([0$NNI])$Stream(NNI),
ADEtestStream(#1, partitions, qDiffDSF(q, #1, #2), 1$UFPSSUPF),
ADEEXPRRStream(#1, #2, partitions, qDiffDX(q, #1, #2, #3), diff1X), _
diffA, diffAF, diffAX, diffC]$HPSPEC
else
error "Guess: guessADE supports only displayAsGF"
ShiftSX(expr: EXPRR, x: Symbol, n: NNI): EXPRR ==
eval(expr, x::EXPRR, x::EXPRR+n::EXPRR)
ShiftSXGF(expr: EXPRR, x: Symbol, n: NNI): EXPRR ==
if zero? n then expr
else
l := [eval(D(expr, x, i)/factorial(i)::EXPRR, x::EXPRR, 0$EXPRR)_
*(x::EXPRR)i for i in 0..n-1]
(expr-reduce(_+, l))/(x::EXPRRn)
ShiftSS(s:UFPSF, n:NNI): UFPSF ==
((quoByVar #1)**n)$MappingPackage1(UFPSF) (s)
ShiftSF(s:UFPSSUPF, n: NNI):UFPSSUPF ==
((quoByVar #1)**n)$MappingPackage1(UFPSSUPF) (s)
ShiftAX(l: NNI, n: Symbol, f: EXPRR): EXPRR ==
n::EXPRR*l f
- ShiftAXGF(l: NNI, x: Symbol, f: EXPRR): EXPRR ==
- I need to help the compiler here, unfortunately
if zero? l then f
else
s := [stirling2(l, i)$IntegerCombinatoricFunctions(Integer)::EXPRR _
(x::EXPRR)iD(f, x, i) for i in 1..l]
reduce(_+, s)
ShiftA(k: NNI, l: NNI, f: SUP S): S ==
ShiftAction(k, l, f)$FFFG(S, SUP S)
ShiftAF(k: NNI, l: NNI, f: UFPSSUPF): SUP F ==
ShiftAction(k, l, f)$FFFG(SUP F, UFPSSUPF)
ShiftC(total: NNI): List S ==
ShiftC(total)$FFFG(S, SUP S)
shiftHP options ==
partitions := FilteredPartitionStream options
if displayAsGF(options)$GOPT0 then
if maxPower(options)$GOPT0 = 1 then
[ADEguessStream(#1, partitions, ShiftSS, (1-monomial(1,1))(-1)),
ADEdegreeStream partitions,
ADEtestStream(#1, partitions, ShiftSF, (1-monomial(1,1))(-1)),
ADEEXPRRStream(#1, #2, partitions, ShiftSXGF, 1/(1-#1::EXPRR)),
ShiftA, ShiftAF, ShiftAXGF, ShiftC]$HPSPEC
else
error "Guess: no support for the Shift operator with displayAsGF _
and maxPower>1"
else
[ADEguessStream2(#1, partitions, ShiftSS, (1-monomial(1,1))(-1)),
ADEdegreeStream partitions,
ADEtestStream2(#1, partitions, ShiftSF, (1-monomial(1,1))(-1)),
ADEEXPRRStream(#1, #2, partitions, ShiftSX, diff1X),
ShiftA, ShiftAF, ShiftAX, ShiftC]$HPSPEC
if F has RetractableTo Symbol and S has RetractableTo Symbol then
qShiftAX(q: Symbol, l: NNI, n: Symbol, f: EXPRR): EXPRR ==
(q::EXPRR)*(ln::EXPRR) * f
qShiftA(q: Symbol, k: NNI, l: NNI, f: SUP S): S ==
qShiftAction(q::S, k, l, f)$FFFG(S, SUP S)
qShiftAF(q: Symbol, k: NNI, l: NNI, f: UFPSSUPF): SUP F ==
qShiftAction(q::F::SUP(F), k, l, f)$FFFG(SUP F, UFPSSUPF)
qShiftC(q: Symbol, total: NNI): List S ==
qShiftC(q::S, total)$FFFG(S, SUP S)
shiftHP(q: Symbol): (LGOPT -> HPSPEC) ==
partitions := FilteredPartitionStream #1
if displayAsGF(#1)$GOPT0 then
error "Guess: no support for the qShift operator with displayAsGF"
else
[ADEguessStream2(#1, partitions, ShiftSS, _
(1-monomial(1,1))(-1)),
ADEdegreeStream partitions,
ADEtestStream2(#1, partitions, ShiftSF, _
(1-monomial(1,1))(-1)),
ADEEXPRRStream(#1, #2, partitions, ShiftSX, diff1X),
qShiftA(q, #1, #2, #3), qShiftAF(q, #1, #2, #3), _
qShiftAX(q, #1, #2, #3), qShiftC(q, #1)]$HPSPEC
makeEXPRR(DAX: DIFFSPECAX, x: Symbol, p: SUP F, expr: EXPRR): EXPRR ==
if zero? p then 0$EXPRR
else
coerce(leadingCoefficient p)::EXPRR * DAX(degree p, x, expr) _
+ makeEXPRR(DAX, x, reductum p, expr)
list2UFPSF(list: List F): UFPSF == series(list::Stream F)$UFPSF
list2UFPSSUPF(list: List F): UFPSSUPF ==
l := [e::SUP(F) for e in list for i in 0..]::Stream SUP F
series(l)$UFPSSUPF + monomial(monomial(1,1)$SUP(F), #list)$UFPSSUPF
SUPF2SUPSUPF(p: SUP F): SUP SUP F ==
map(#1::SUP F, p)$SparseUnivariatePolynomialFunctions2(F, SUP F)
UFPSF2SUPF(f: UFPSF, deg: NNI): SUP F ==
makeSUP univariatePolynomial(f, deg)
getListSUPF(s: Stream UFPSF, o: NNI, deg: NNI): List SUP F ==
map(UFPSF2SUPF(#1, deg), entries complete first(s, o))
$ListFunctions2(UFPSF, SUP F)
S2EXPRR(s: S): EXPRR ==
if F is S then
coerce(s pretend F)@EXPRR
else if F is Fraction S then
coerce(s::Fraction(S))@EXPRR
else error "Type parameter F should be either equal to S or equal _
to Fraction S"
guessInterpolate(guessList: List SUP F, eta: List NNI, D: HPSPEC)
: Matrix SUP S ==
if F is S then
vguessList: Vector SUP S := vector(guessList pretend List(SUP(S)))
generalInterpolation((D.C)(reduce(_+, eta)), D.A,
vguessList, eta)$FFFG(S, SUP S)
else if F is Fraction S then
vguessListF: Vector SUP F := vector(guessList)
generalInterpolation((D.C)(reduce(_+, eta)), D.A,
vguessListF, eta)$FFFGF(S, SUP S, SUP F) else error "Type parameter F should be either equal to S or equal _
to Fraction S"
guessInterpolate2(guessList: List SUP F,
sumEta: NNI, maxEta: NNI,
D: HPSPEC): Stream Matrix SUP S ==
if F is S then
vguessList: Vector SUP S := vector(guessList pretend List(SUP(S)))
generalInterpolation((D.C)(sumEta), D.A,
vguessList, sumEta, maxEta)
$FFFG(S, SUP S)
else if F is Fraction S then
vguessListF: Vector SUP F := vector(guessList)
generalInterpolation((D.C)(sumEta), D.A,
vguessListF, sumEta, maxEta)
$FFFGF(S, SUP S, SUP F)
else error "Type parameter F should be either equal to S or equal _
to Fraction S"
testInterpolant(resi: List SUP S,
list: List F,
testList: List UFPSSUPF,
exprList: List EXPRR,
initials: List EXPRR,
guessDegree: NNI,
D: HPSPEC,
dummy: Symbol, op: BasicOperator, options: LGOPT,
list: List F)
: Union("failed", Record(function: EXPRR, order: NNI)) ==
((maxDegree(options)$GOPT0 = -1) and
(allDegrees(options)$GOPT0 = false) and
zero?(last resi))
=> return "failed"
nonZeroCoefficient: Integer := 0
for i in 1..#resi repeat
if not zero? resi.i then
if zero? nonZeroCoefficient then
nonZeroCoefficient := i
else
nonZeroCoefficient := 0
break
if not zero? nonZeroCoefficient then
(freeOf?(exprList.nonZeroCoefficient, name op)) => return "failed" for e in list repeat
if not zero? e then return "failed"
else
resiSUPF := map(SUPF2SUPSUPF SUPS2SUPF #1, resi)
$ListFunctions2(SUP S, SUP SUP F)
iterate? := true;
for d in guessDegree+1.. repeat
c: SUP F := generalCoefficient(D.AF, vector testList,
d, vector resiSUPF)
$FFFG(SUP F, UFPSSUPF) if not zero? c then
iterate? := ground? c
break
iterate? => return "failed"
g: SUP S := gcd resi
resiF := map(SUPS2SUPF((#1 exquo g)::SUP(S)), resi)
$ListFunctions2(SUP S, SUP F)
if debug(options)$GOPT0 then
output(hconcat("trying possible solution ", resiF::OutputForm))
$OutputPackage
-- transform each term into an expression
ex: List EXPRR := [makeEXPRR(D.AX, dummy, p, e) _
for p in resiF for e in exprList]
-- transform the list of expressions into a sum of expressions
res: EXPRR
if displayAsGF(options)$GOPT0 then
res := evalADE(op, dummy, variableName(options)$GOPT0::EXPRR,
indexName(options)$GOPT0::EXPRR,
numerator reduce(_+, ex),
reverse initials)
$RecurrenceOperator(Integer, EXPRR)
ord: NNI := 0
-- FIXME: checkResult doesn't really work yet for generating functions
else
res := evalRec(op, dummy, indexName(options)$GOPT0::EXPRR,
indexName(options)$GOPT0::EXPRR,
numerator reduce(+, ex),
reverse initials)
$RecurrenceOperator(Integer, EXPRR)
ord: NNI := checkResult(res, indexName(options)$GOPT0,
#list, list, options)
[res, ord]$Record(function: EXPRR, order: NNI)
guessHPaux(list: List F, D: HPSPEC, options: LGOPT): GUESSRESULT ==
reslist: GUESSRESULT := [] listDegree := #list-1-safety(options)$GOPT0
if listDegree < 0 then return reslist
a := functionName(options)$GOPT0
op := operator a
x := variableName(options)$GOPT0
dummy := new$Symbol
initials: List EXPRR
if displayAsGF(options)$GOPT0
then initials := [coerce(e)@EXPRR*factorial(k::EXPRR) _
for e in list for k in 0..]
else initials := [coerce(e)@EXPRR for e in list]
guessS := (D.guessStream)(list2UFPSF list)
degreeS := D.degreeStream
testS := (D.testStream)(list2UFPSSUPF list)
exprS := (D.exprStream)(op(dummy::EXPRR)::EXPRR, dummy)
iterate?: Boolean := false -- this is necessary because the compiler
-- doesn't understand => "iterate" properly
-- the latter just leaves the current block, it
-- seems
for o in 2.. repeat
empty? rest(guessS, (o-1)::NNI) => break
guessDegree: Integer := listDegree-(degreeS.o)::Integer
guessDegree < 0 => break
if debug(options)$GOPT0 then
output(hconcat("Trying order ", o::OutputForm))$OutputPackage
output(hconcat("guessDegree is ", guessDegree::OutputForm))
$OutputPackage
if allDegrees(options)$GOPT0 then
(o > guessDegree+2) => return reslist maxEta: Integer := 1+maxDegree(options)$GOPT0
if maxEta = 0
then maxEta := guessDegree+1
else
maxParams := divide(guessDegree::NNI+1, o)
if debug(options)$GOPT0
then output(hconcat("maxParams: ", maxParams::OutputForm))
$OutputPackage
if maxParams.quotient = 0 and maxParams.remainder < o-1
then return reslist
if ((maxDegree(options)$GOPT0 ~= -1) and
(maxDegree(options)$GOPT0 < maxParams.quotient)) and
not (empty? rest(guessS, o) or
((newGuessDegree := listDegree-(degreeS.(o+1))::Integer)
< 0) or
(((newMaxParams := divide(newGuessDegree::NNI+1, o+1))
.quotient = 0) and
(newMaxParams.remainder < o)))
then iterate? := true
else eta: List NNI
:= [(if i <= maxParams.remainder
then maxParams.quotient + 1
else maxParams.quotient)::NNI for i in 1..o]
if iterate?
then
iterate? := false
if debug(options)$GOPT0 then output("iterating")$OutputPackage
else
guessList: List SUP F := getListSUPF(guessS, o, guessDegree::NNI)
testList: List UFPSSUPF := entries complete first(testS, o)
exprList: List EXPRR := entries complete first(exprS, o)
if debug(options)$GOPT0 then
output("The list of expressions is")$OutputPackage
output(exprList::OutputForm)$OutputPackage
if allDegrees(options)$GOPT0 then
MS: Stream Matrix SUP S := guessInterpolate2(guessList,
guessDegree::NNI+1,
maxEta::NNI, D)
repeat
(empty? MS) => break
M := first MS for i in 1..o repeat
res := testInterpolant(entries column(M, i),
list,
testList,
exprList,
initials,
guessDegree::NNI,
D, dummy, op, options, list) (res case "failed") => "iterate"
if not member?(res, reslist)
then reslist := cons(res, reslist)
if one(options)$GOPT0 then return reslist
MS := rest MS
else
M: Matrix SUP S := guessInterpolate(guessList, eta, D)
for i in 1..o repeat
res := testInterpolant(entries column(M, i),
list,
testList,
exprList,
initials,
guessDegree::NNI,
D, dummy, op, options, list)
(res case "failed") => "iterate" if not member?(res, reslist)
then reslist := cons(res, reslist)
if one(options)$GOPT0 then return reslist
reslist
guessHP(D: LGOPT -> HPSPEC): GUESSER == guessHPaux(#1, D #2, #2)
guessADE(list: List F, options: LGOPT): GUESSRESULT ==
opts: LGOPT := cons(displayAsGF(true)$GuessOption, options)
guessHPaux(list, diffHP opts, opts)
guessADE(list: List F): GUESSRESULT == guessADE(list, [])
guessAlg(list: List F, options: LGOPT) ==
guessADE(list, cons(maxDerivative(0)$GuessOption, options))
guessAlg(list: List F): GUESSRESULT == guessAlg(list, [])
guessHolo(list: List F, options: LGOPT): GUESSRESULT ==
guessADE(list, cons(maxPower(1)$GuessOption, options))
guessHolo(list: List F): GUESSRESULT == guessHolo(list, [])
guessPade(list: List F, options: LGOPT): GUESSRESULT ==
opts := append(options, [maxDerivative(0)$GuessOption,
maxPower(1)$GuessOption,
allDegrees(true)$GuessOption])
guessADE(list, opts)
guessPade(list: List F): GUESSRESULT == guessPade(list, [])
if F has RetractableTo Symbol and S has RetractableTo Symbol then guessADE(q: Symbol): GUESSER ==
opts: LGOPT := cons(displayAsGF(true)$GuessOption, #2)
guessHPaux(#1, (diffHP q)(opts), opts)
guessRec(list: List F, options: LGOPT): GUESSRESULT ==
opts: LGOPT := cons(displayAsGF(false)$GuessOption, options)
guessHPaux(list, shiftHP opts, opts)
guessRec(list: List F): GUESSRESULT == guessRec(list, [])
guessPRec(list: List F, options: LGOPT): GUESSRESULT ==
guessRec(list, cons(maxPower(1)$GuessOption, options))
guessPRec(list: List F): GUESSRESULT == guessPRec(list, [])
guessRat(list: List F, options: LGOPT): GUESSRESULT ==
opts := append(options, [maxShift(0)$GuessOption,
maxPower(1)$GuessOption,
allDegrees(true)$GuessOption])
guessRec(list, opts)
guessRat(list: List F): GUESSRESULT == guessRat(list, [])
if F has RetractableTo Symbol and S has RetractableTo Symbol then
guessRec(q: Symbol): GUESSER ==
opts: LGOPT := cons(displayAsGF(false)$GuessOption, #2)
guessHPaux(#1, (shiftHP q)(opts), opts)
guessPRec(q: Symbol): GUESSER ==
opts: LGOPT := append([displayAsGF(false)$GuessOption,
maxPower(1)$GuessOption], #2)
guessHPaux(#1, (shiftHP q)(opts), opts)
guessRat(q: Symbol): GUESSER ==
opts := append(#2, [displayAsGF(false)$GuessOption,
maxShift(0)$GuessOption,
maxPower(1)$GuessOption,
allDegrees(true)$GuessOption])
guessHPaux(#1, (shiftHP q)(opts), opts)
guess(list: List F, guessers: List GUESSER, ops: List Symbol,
options: LGOPT): GUESSRESULT ==
maxLevel := maxLevel(options)$GOPT0
xx := indexName(options)$GOPT0
if debug(options)$GOPT0 then
output(hconcat("guessing level ", maxLevel::OutputForm))
$OutputPackage
output(hconcat("guessing ", list::OutputForm))$OutputPackage
res: GUESSRESULT := []
len := #list :: PositiveInteger
if len <= 1 then return res for guesser in guessers repeat
res := append(guesser(list, options), res) if debug(options)$GOPT0 then
output(hconcat("res ", res::OutputForm))$OutputPackage
if one(options)$GOPT0 and not empty? res then return res
if (maxLevel = 0) then return res
if member?('guessProduct, ops) and not member?(0$F, list) then
prodList: List F := [(list.(i+1))/(list.i) for i in 1..(len-1)] if not every?(one?, prodList) then
var: Symbol := subscript('p, [len::OutputForm])
prodGuess :=
[[coerce(list.(guess.order+1))
* product(guess.function,
equation(var,
(guess.order)::EXPRR..xx::EXPRR-1)), _
guess.order] _
for guess in guess(prodList, guessers, ops,
append([indexName(var)$GuessOption,
maxLevel(maxLevel-1)
$GuessOption],
options))$%] if debug(options)$GOPT0 then
output(hconcat(prodGuess ,
prodGuess::OutputForm))
$OutputPackage
for guess in prodGuess
| not any?(guess.function = #1.function, res) repeat
res := cons(guess, res)
if one(options)$GOPT0 and not empty? res then return res
if member?('guessSum, ops) then
sumList: List F := [(list.(i+1))-(list.i) for i in 1..(len-1)] if not every?(#1=sumList.1, sumList) then
var: Symbol := subscript('s, [len::OutputForm])
sumGuess :=
[[coerce(list.(guess.order+1))
+ summation(guess.function,
equation(var,
(guess.order)::EXPRR..xx::EXPRR-1)),
guess.order] _
for guess in guess(sumList, guessers, ops,
append([indexName(var)$GuessOption,
maxLevel(maxLevel-1)
$GuessOption],
options))$%] for guess in sumGuess
| not any?(guess.function = #1.function, res) repeat
res := cons(guess, res)
res
guess(list: List F): GUESSRESULT ==
guess(list, [guessADE$%, guessRec$%], ['guessProduct, 'guessSum], [])
guess(list: List F, options: LGOPT): GUESSRESULT ==
guess(list, [guessADE$%, guessRat$%], ['guessProduct, 'guessSum],
options)
guess(list: List F, guessers: List GUESSER, ops: List Symbol)
: GUESSRESULT ==
guess(list, guessers, ops, [])
\end{spad}
- \begin{spad}
- concerning algebraic functions, it may make sense to look at A.J.van der
-- Poorten, Power Series Representing Algebraic Functions, Section 6.
)abbrev package GUESSINT GuessInteger
++ Description:
++ This package exports guessing of sequences of rational numbers
GuessInteger() == Guess(Fraction Integer, Integer, Expression Integer,
Fraction Integer,
id$MappingPackage1(Fraction Integer),
coerce$Expression(Integer))
)abbrev package GUESSF1 GuessFiniteFunctions
++ Description:
++ This package exports guessing of sequences of numbers in a finite field
GuessFiniteFunctions(F:Join(FiniteFieldCategory, ConvertibleTo Integer)):
Exports == Implementation where EXPRR ==> Expression Integer
Exports == with
F2EXPRR: F -> EXPRR
Implementation == add
F2EXPRR(p: F): EXPRR == convert(p)@Integer::EXPRR
)abbrev package GUESSF GuessFinite
++ Description:
++ This package exports guessing of sequences of numbers in a finite field
GuessFinite(F:Join(FiniteFieldCategory, ConvertibleTo Integer)) ==
Guess(F, F, Expression Integer, Integer, coerce$F,
F2EXPRR$GuessFiniteFunctions(F))
)abbrev package GUESSP GuessPolynomial
++ Description:
++ This package exports guessing of sequences of rational functions
GuessPolynomial() == Guess(Fraction Polynomial Integer, Polynomial Integer,
Expression Integer,
Fraction Polynomial Integer,
id$MappingPackage1(Fraction Polynomial Integer),
coerce$Expression(Integer))
\end{spad}
Some or all expressions may not have rendered properly,
because Axiom returned the following error:
Error: export AXIOM=/usr/local/lib/axiom/target/x86_64-unknown-linux; export ALDORROOT=/usr/local/aldor/linux/1.1.0; export PATH=$ALDORROOT/bin:$PATH; export HOME=/var/zope2/var/LatexWiki; ulimit -t 240; $AXIOM/bin/AXIOMsys < /var/zope2/var/LatexWiki/7164797305649977448-25px.axm
GCL (GNU Common Lisp) 2.6.8 CLtL1 Nov 9 2007 07:47:56
Source License: LGPL(gcl,gmp), GPL(unexec,bfd,xgcl)
Binary License: GPL due to GPL'ed components: (READLINE BFD UNEXEC)
Modifications of this banner must retain notice of a compatible license
Dedicated to the memory of W. Schelter
Use (help) to get some basic information on how to use GCL.
Temporary directory for compiler files set to /tmp/
FriCAS (AXIOM fork) Computer Algebra System
Version: FriCAS 2007-10-02
Timestamp: Friday November 9, 2007 at 19:35:06
-----------------------------------------------------------------------------
Issue )copyright to view copyright notices.
Issue )summary for a summary of useful system commands.
Issue )quit to leave FriCAS and return to shell.
-----------------------------------------------------------------------------
(1) -> (1) -> (1) -> (1) -> (1) -> )lib RECOP FAMR2 FFFG FFFGF NEWTON GOPT GOPT0 UFPS UFPS1
RecurrenceOperator is now explicitly exposed in frame initial
RecurrenceOperator will be automatically loaded when needed from
/var/zope2/var/LatexWiki/RECOP.NRLIB/code
FiniteAbelianMonoidRingFunctions2 is now explicitly exposed in frame
initial
FiniteAbelianMonoidRingFunctions2 will be automatically loaded when
needed from /var/zope2/var/LatexWiki/FAMR2.NRLIB/code
FractionFreeFastGaussian is now explicitly exposed in frame initial
FractionFreeFastGaussian will be automatically loaded when needed
from /var/zope2/var/LatexWiki/FFFG.NRLIB/code
FractionFreeFastGaussianFractions is now explicitly exposed in frame
initial
FractionFreeFastGaussianFractions will be automatically loaded when
needed from /var/zope2/var/LatexWiki/FFFGF.NRLIB/code
NewtonInterpolation is now explicitly exposed in frame initial
NewtonInterpolation will be automatically loaded when needed from
/var/zope2/var/LatexWiki/NEWTON.NRLIB/code
GuessOption is now explicitly exposed in frame initial
GuessOption will be automatically loaded when needed from
/var/zope2/var/LatexWiki/GOPT.NRLIB/code
GuessOptionFunctions0 is now explicitly exposed in frame initial
GuessOptionFunctions0 will be automatically loaded when needed from
/var/zope2/var/LatexWiki/GOPT0.NRLIB/code
UnivariateFormalPowerSeries is now explicitly exposed in frame
initial
UnivariateFormalPowerSeries will be automatically loaded when needed
from /var/zope2/var/LatexWiki/UFPS.NRLIB/code
UnivariateFormalPowerSeriesFunctions is now explicitly exposed in
frame initial
UnivariateFormalPowerSeriesFunctions will be automatically loaded
when needed from /var/zope2/var/LatexWiki/UFPS1.NRLIB/code
(1) -> <spad>
)abbrev package GUESS Guess
++ Author: Martin Rubey
++ Description: This package implements guessing of sequences. Packages for the
++ most common cases are provided as \spadtype{GuessInteger},
++ \spadtype{GuessPolynomial}, etc.
Guess(F, S, EXPRR, R, retract, coerce): Exports == Implementation where
F: Field -- zB.: FRAC POLY PF 5
-- in F wird interpoliert und gerechnet S: GcdDomain
- -- in guessExpRat möchte ich die Wurzeln von Polynomen über F bestimmen. Wenn F
- ein Quotientenkörper ist, dann kann ich den Nenner loswerden.
-- F wäre dann also QFCAT S
- R: Join(OrderedSet, IntegralDomain)
- zB.: FRAC POLY INT
- -- die Ergebnisse werden aber in EXPRR dargestellt
- EXPRR: Join(ExpressionSpace, IntegralDomain,
EXPRR: Join(FunctionSpace Integer, IntegralDomain,
RetractableTo R, RetractableTo Symbol,
RetractableTo Integer, CombinatorialOpsCategory,
PartialDifferentialRing Symbol) with
_ : (%,%) -> %
/ : (%,%) -> %
_* : (%,%) -> %
numerator : % -> %
denominator : % -> %
ground? : % -> Boolean
-- zB.: EXPR INT
-- EXPR FRAC POLY INT ist ja verboten. Deshalb kann ich nicht einfach EXPR R
-- verwenden
-- EXPRR existiert, falls irgendwann einmal EXPR PF 5 unterstützt wird.
-- das folgende möchte ich eigentlich loswerden.
- retract: R -> F
- zB.: i+->i
coerce: F -> EXPRR -- zB.: i+->i
-- Achtung EXPRR ~= EXPR R
LGOPT ==> List GuessOption
GOPT0 ==> GuessOptionFunctions0
NNI ==> NonNegativeInteger
PI ==> PositiveInteger
EXPRI ==> Expression Integer
GUESSRESULT ==> List Record(function: EXPRR, order: NNI)
UFPSF ==> UnivariateFormalPowerSeries F
UFPS1 ==> UnivariateFormalPowerSeriesFunctions
UFPSS ==> UnivariateFormalPowerSeries S
SUP ==> SparseUnivariatePolynomial
UFPSSUPF ==> UnivariateFormalPowerSeries SUP F
FFFG ==> FractionFreeFastGaussian
FFFGF ==> FractionFreeFastGaussianFractions
-- CoeffAction
DIFFSPECA ==> (NNI, NNI, SUP S) -> S DIFFSPECAF ==> (NNI, NNI, UFPSSUPF) -> SUP F
DIFFSPECAX ==> (NNI, Symbol, EXPRR) -> EXPRR
-- the diagonal of the C-matrix
DIFFSPECC ==> NNI -> List S HPSPEC ==> Record(guessStream: UFPSF -> Stream UFPSF,
degreeStream: Stream NNI,
testStream: UFPSSUPF -> Stream UFPSSUPF,
exprStream: (EXPRR, Symbol) -> Stream EXPRR,
A: DIFFSPECA,
AF: DIFFSPECAF,
AX: DIFFSPECAX,
C: DIFFSPECC)
- -- note that empty?(guessStream.o) has to return always. In other words, if the
- stream is finite, empty? should recognize it.
- DIFFSPECN ==> EXPRR -> EXPRR
- eg.: i+->q^i
GUESSER ==> (List F, LGOPT) -> GUESSRESULT
FSUPS ==> Fraction SUP S
FSUPF ==> Fraction SUP F
V ==> OrderedVariableList(['a1, 'A])
POLYF ==> SparseMultivariatePolynomial(F, V)
FPOLYF ==> Fraction POLYF
FSUPFPOLYF ==> Fraction SUP FPOLYF
POLYS ==> SparseMultivariatePolynomial(S, V)
FPOLYS ==> Fraction POLYS
FSUPFPOLYS ==> Fraction SUP FPOLYS
--<<implementation: Guess - Hermite-Pade - Types for Operators>>
Exports == with guess: List F -> GUESSRESULT
++ \spad{guess l} applies recursively \spadfun{guessRec} and
++ \spadfun{guessADE} to the successive differences and quotients of
++ the list. Default options as described in
++ \spadtype{GuessOptionFunctions0} are used.
guess: (List F, LGOPT) -> GUESSRESULT
++ \spad{guess(l, options)} applies recursively \spadfun{guessRec}
++ and \spadfun{guessADE} to the successive differences and quotients
++ of the list. The given options are used.
guess: (List F, List GUESSER, List Symbol) -> GUESSRESULT
++ \spad{guess(l, guessers, ops)} applies recursively the given
++ guessers to the successive differences if ops contains the symbol
++ guessSum and quotients if ops contains the symbol guessProduct to
++ the list. Default options as described in
++ \spadtype{GuessOptionFunctions0} are used.
guess: (List F, List GUESSER, List Symbol, LGOPT) -> GUESSRESULT
++ \spad{guess(l, guessers, ops)} applies recursively the given
++ guessers to the successive differences if ops contains the symbol
++ \spad{guessSum} and quotients if ops contains the symbol
++ \spad{guessProduct} to the list. The given options are used.
guessExpRat: List F -> GUESSRESULT
++ \spad{guessExpRat l} tries to find a function of the form
++ n+->(a+b n)^n r(n), where r(n) is a rational function, that fits
++ l.
guessExpRat: (List F, LGOPT) -> GUESSRESULT
++ \spad{guessExpRat(l, options)} tries to find a function of the
++ form n+->(a+b n)^n r(n), where r(n) is a rational function, that
++ fits l.
if F has RetractableTo Symbol and S has RetractableTo Symbol then
guessExpRat: Symbol -> GUESSER
++ \spad{guessExpRat q} returns a guesser that tries to find a
++ function of the form n+->(a+b q^n)^n r(q^n), where r(n) is a
++ rational function, that fits l.
guessHP: (LGOPT -> HPSPEC) -> GUESSER
++ \spad{guessHP f} constructs an operation that applies Hermite-Pade
++ approximation to the series generated by the given function f.
guessADE: List F -> GUESSRESULT
++ \spad{guessADE l} tries to find an algebraic differential equation
++ for a generating function whose first Taylor coefficients are
++ given by l, using the default options described in
++ \spadtype{GuessOptionFunctions0}.
guessADE: (List F, LGOPT) -> GUESSRESULT
++ \spad{guessADE(l, options)} tries to find an algebraic
++ differential equation for a generating function whose first Taylor
++ coefficients are given by l, using the given options.
guessAlg: List F -> GUESSRESULT
++ \spad{guessAlg l} tries to find an algebraic equation for a
++ generating function whose first Taylor coefficients are given by
++ l, using the default options described in
++ \spadtype{GuessOptionFunctions0}. It is equivalent to
++ \spadfun{guessADE}(l, maxDerivative == 0).
guessAlg: (List F, LGOPT) -> GUESSRESULT
++ \spad{guessAlg(l, options)} tries to find an algebraic equation
++ for a generating function whose first Taylor coefficients are
++ given by l, using the given options. It is equivalent to
++ \spadfun{guessADE}(l, options) with \spad{maxDerivative == 0}.
guessHolo: List F -> GUESSRESULT
++ \spad{guessHolo l} tries to find an ordinary linear differential
++ equation for a generating function whose first Taylor coefficients
++ are given by l, using the default options described in
++ \spadtype{GuessOptionFunctions0}. It is equivalent to
++ \spadfun{guessADE}\spad{(l, maxPower == 1)}.
guessHolo: (List F, LGOPT) -> GUESSRESULT
++ \spad{guessHolo(l, options)} tries to find an ordinary linear
++ differential equation for a generating function whose first Taylor
++ coefficients are given by l, using the given options. It is
++ equivalent to \spadfun{guessADE}\spad{(l, options)} with
++ \spad{maxPower == 1}.
guessPade: (List F, LGOPT) -> GUESSRESULT
++ \spad{guessPade(l, options)} tries to find a rational function
++ whose first Taylor coefficients are given by l, using the given
++ options. It is equivalent to \spadfun{guessADE}\spad{(l,
++ maxDerivative == 0, maxPower == 1, allDegrees == true)}.
guessPade: List F -> GUESSRESULT
++ \spad{guessPade(l, options)} tries to find a rational function
++ whose first Taylor coefficients are given by l, using the default
++ options described in \spadtype{GuessOptionFunctions0}. It is
++ equivalent to \spadfun{guessADE}\spad{(l, options)} with
++ \spad{maxDerivative == 0, maxPower == 1, allDegrees == true}.
guessRec: List F -> GUESSRESULT
++ \spad{guessRec l} tries to find an ordinary difference equation
++ whose first values are given by l, using the default options
++ described in \spadtype{GuessOptionFunctions0}.
guessRec: (List F, LGOPT) -> GUESSRESULT
++ \spad{guessRec(l, options)} tries to find an ordinary difference
++ equation whose first values are given by l, using the given
++ options.
guessPRec: (List F, LGOPT) -> GUESSRESULT
++ \spad{guessPRec(l, options)} tries to find a linear recurrence
++ with polynomial coefficients whose first values are given by l,
++ using the given options. It is equivalent to
++ \spadfun{guessRec}\spad{(l, options)} with \spad{maxPower == 1}.
guessPRec: List F -> GUESSRESULT
++ \spad{guessPRec l} tries to find a linear recurrence with
++ polynomial coefficients whose first values are given by l, using
++ the default options described in
++ \spadtype{GuessOptionFunctions0}. It is equivalent to
++ \spadfun{guessRec}\spad{(l, maxPower == 1)}.
guessRat: (List F, LGOPT) -> GUESSRESULT
++ \spad{guessRat(l, options)} tries to find a rational function
++ whose first values are given by l, using the given options. It is
++ equivalent to \spadfun{guessRec}\spad{(l, maxShift == 0, maxPower
++ == 1, allDegrees == true)}.
guessRat: List F -> GUESSRESULT
++ \spad{guessRat l} tries to find a rational function whose first
++ values are given by l, using the default options described in
++ \spadtype{GuessOptionFunctions0}. It is equivalent to
++ \spadfun{guessRec}\spad{(l, maxShift == 0, maxPower == 1,
++ allDegrees == true)}.
diffHP: LGOPT -> HPSPEC
++ \spad{diffHP options} returns a specification for Hermite-Pade
++ approximation with the differential operator
shiftHP: LGOPT -> HPSPEC
++ \spad{shiftHP options} returns a specification for Hermite-Pade
++ approximation with the shift operator
if F has RetractableTo Symbol and S has RetractableTo Symbol then
shiftHP: Symbol -> (LGOPT -> HPSPEC)
++ \spad{shiftHP options} returns a specification for
++ Hermite-Pade approximation with the $q$-shift operator diffHP: Symbol -> (LGOPT -> HPSPEC)
++ \spad{diffHP options} returns a specification for Hermite-Pade
++ approximation with the $q$-dilation operator
guessRec: Symbol -> GUESSER
++ \spad{guessRec q} returns a guesser that finds an ordinary
++ q-difference equation whose first values are given by l, using
++ the given options.
guessPRec: Symbol -> GUESSER
++ \spad{guessPRec q} returns a guesser that tries to find
++ a linear q-recurrence with polynomial coefficients whose first
++ values are given by l, using the given options. It is
++ equivalent to \spadfun{guessRec}\spad{(q)} with
++ \spad{maxPower == 1}.
guessRat: Symbol -> GUESSER
++ \spad{guessRat q} returns a guesser that tries to find a
++ q-rational function whose first values are given by l, using
++ the given options. It is equivalent to \spadfun{guessRec} with
++ \spad{(l, maxShift == 0, maxPower == 1, allDegrees == true)}.
guessADE: Symbol -> GUESSER
++ \spad{guessADE q} returns a guesser that tries to find an
++ algebraic differential equation for a generating function whose
++ first Taylor coefficients are given by l, using the given
++ options.
--<<debug exports: Guess>>
Implementation == add
-- We have to put this chunk at the beginning, because otherwise it will take
-- very long to compile.
ord1(x: List Integer, i: Integer): Integer ==
n := #x - 3 - i
x.(n+1)*reduce(_+, [x.j for j in 1..n], 0) + _
2*reduce(_+, [reduce(_+, [x.k*x.j for k in 1..j-1], 0) _
for j in 1..n], 0)
ord2(x: List Integer, i: Integer): Integer ==
if zero? i then
n := #x - 3 - i
ord1(x, i) + reduce(_+, [x.j for j in 1..n], 0)*(x.(n+2)-x.(n+1))
else
ord1(x, i)
deg1(x: List Integer, i: Integer): Integer ==
m := #x - 3
(x.(m+3)+x.(m+1)+x.(1+i))*reduce(_+, [x.j for j in 2+i..m], 0) + _
x.(m+3)x.(m+1) + _
2reduce(_+, [reduce(_+, [x.k*x.j for k in 2+i..j-1], 0) _
for j in 2+i..m], 0)
deg2(x: List Integer, i: Integer): Integer ==
m := #x - 3
deg1(x, i) + _
(x.(m+3) + reduce(_+, [x.j for j in 2+i..m], 0)) * _
(x.(m+2)-x.(m+1))
checkResult(res: EXPRR, n: Symbol, l: Integer, list: List F,
options: LGOPT): NNI ==
for i in l..1 by -1 repeat
den := eval(denominator res, n::EXPRR, (i-1)::EXPRR)
if den = 0 then return i::NNI
num := eval(numerator res, n::EXPRR, (i-1)::EXPRR)
if list.i ~= retract(retract(num/den)@R)
then return i::NNI
0$NNI
SUPS2SUPF(p: SUP S): SUP F ==
if F is S then
p pretend SUP(F)
else if F is Fraction S then
map(coerce(#1)$Fraction(S), p)
$SparseUnivariatePolynomialFunctions2(S, F)
else error "Type parameter F should be either equal to S or equal _
to Fraction S"
F2FPOLYS(p: F): FPOLYS ==
if F is S then
p::POLYF::FPOLYF pretend FPOLYS
else if F is Fraction S then
numer(p)$Fraction(S)::POLYS/denom(p)$Fraction(S)::POLYS
else error "Type parameter F should be either equal to S or equal _
to Fraction S"
MPCSF ==> MPolyCatFunctions2(V, IndexedExponents V,
IndexedExponents V, S, F,
POLYS, POLYF)
SUPF2EXPRR(xx: Symbol, p: SUP F): EXPRR ==
zero? p => 0
(coerce(leadingCoefficient p))::EXPRR (xx::EXPRR)*degree p
+ SUPF2EXPRR(xx, reductum p)
FSUPF2EXPRR(xx: Symbol, p: FSUPF): EXPRR ==
(SUPF2EXPRR(xx, numer p)) / (SUPF2EXPRR(xx, denom p))
POLYS2POLYF(p: POLYS): POLYF ==
if F is S then
p pretend POLYF
else if F is Fraction S then
map(coerce(#1)$Fraction(S), p)$MPCSF
else error "Type parameter F should be either equal to S or equal _
to Fraction S"
SUPPOLYS2SUPF(p: SUP POLYS, a1v: F, Av: F): SUP F ==
zero? p => 0
lc: POLYF := POLYS2POLYF leadingCoefficient(p)
monomial(retract(eval(lc, [index(1)$V, index(2)$V]::List V,
[a1v, Av])),
degree p)
+ SUPPOLYS2SUPF(reductum p, a1v, Av)
SUPFPOLYS2FSUPPOLYS(p: SUP FPOLYS): Fraction SUP POLYS ==
cden := splitDenominator(p)
$UnivariatePolynomialCommonDenominator(POLYS, FPOLYS,SUP FPOLYS) pnum: SUP POLYS
:= map(retract(#1 * cden.den)$FPOLYS, p)
$SparseUnivariatePolynomialFunctions2(FPOLYS, POLYS)
pden: SUP POLYS := (cden.den)::SUP POLYS
pnum/pden
POLYF2EXPRR(p: POLYF): EXPRR ==
map(convert(#1)@Symbol::EXPRR, coerce(#1)@EXPRR, p)
$PolynomialCategoryLifting(IndexedExponents V, V,
F, POLYF, EXPRR)
defaultD: DIFFSPECN
defaultD(expr: EXPRR): EXPRR == expr
-- applies n+->q^n or whatever DN is to i
DN2DL: (DIFFSPECN, Integer) -> F
DN2DL(DN, i) == retract(retract(DN(i::EXPRR))@R)
evalResultant(p1: POLYS, p2: POLYS, o: Integer, d: Integer, va1: V, vA: V)_
: List S ==
res: List S := []
d1 := degree(p1, va1)
d2 := degree(p2, va1)
lead: S
for k in 1..d-o+1 repeat
p1atk := univariate eval(p1, vA, k::S)
p2atk := univariate eval(p2, vA, k::S) d1atk := degree p1atk
d2atk := degree p2atk
-- output("k: " string(k))$OutputPackage
-- output("d1: " string(d1) " d1atk: " string(d1atk))$OutputPackage
-- output("d2: " string(d2) " d2atk: " string(d2atk))$OutputPackage if d2atk < d2 then
if d1atk < d1
then lead := 0$S
else lead := (leadingCoefficient p1atk)((d2-d2atk)::NNI)
else
if d1atk < d1
then lead := (-1$S)d2 (leadingCoefficient p2atk)*((d1-d1atk)::NNI)
else lead := 1$S
if zero? lead
then res := cons(0, res)
else res := cons(lead (resultant(p1atk, p2atk)$SUP(S) exquo _
(k::S)*(o::NNI))::S,
res)
reverse res
p(xm: Integer, i: Integer, va1: V, vA: V, basis: DIFFSPECN): FPOLYS ==
vA::POLYS::FPOLYS + va1::POLYS::FPOLYS _
* F2FPOLYS(DN2DL(basis, i) - DN2DL(basis, xm))
p2(xm: Integer, i: Symbol, a1v: F, Av: F, basis: DIFFSPECN): EXPRR ==
coerce(Av) + coerce(a1v)*(basis(i::EXPRR) - basis(xm::EXPRR))
GF ==> GeneralizedMultivariateFactorize(SingletonAsOrderedSet,
IndexedExponents V, F, F,
SUP F)
guessExpRatAux(xx: Symbol, list: List F, basis: DIFFSPECN,
xValues: List Integer, options: LGOPT): List EXPRR == a1: V := index(1)$V
A: V := index(2)$V
len: NNI := #list
if len < 4 then return []
else len := (len-3)::NNI
xlist := [F2FPOLYS DN2DL(basis, xValues.i) for i in 1..len]
x1 := F2FPOLYS DN2DL(basis, xValues.(len+1))
x2 := F2FPOLYS DN2DL(basis, xValues.(len+2))
x3 := F2FPOLYS DN2DL(basis, xValues.(len+3))
y: NNI -> FPOLYS :=
F2FPOLYS(list.#1) _
p(last xValues, (xValues.#1)::Integer, a1, A, basis)*_
(-(xValues.#1)::Integer)
ylist: List FPOLYS := [y i for i in 1..len]
y1 := y(len+1)
y2 := y(len+2)
y3 := y(len+3)
res := []::List EXPRR
if maxDegree(options)$GOPT0 = -1
then maxDeg := len-1
else maxDeg := min(maxDegree(options)$GOPT0, len-1)
for i in 0..maxDeg repeat
if debug(options)$GOPT0 then
output(hconcat(degree ExpRat , i::OutputForm))
$OutputPackage if debug(options)$GOPT0 then
systemCommand("sys date +%s")$MoreSystemCommands
output(interpolating...)$OutputPackage
ri: FSUPFPOLYS
:= interpolate(xlist, ylist, (len-1-i)::NNI) _
$FFFG(FPOLYS, SUP FPOLYS)
-- for experimental fraction free interpolation
-- ri: Fraction SUP POLYS
-- := interpolate(xlist, ylist, (len-1-i)::NNI) _
-- $FFFG(POLYS, SUP POLYS)
- if debug(options)$GOPT0 then
- output(hconcat("xlist: ", xlist::OutputForm))$OutputPackage
-- output(hconcat("ylist: ", ylist::OutputForm))$OutputPackage
-- output(hconcat("ri: ", ri::OutputForm))$OutputPackage
systemCommand("sys date +%s")$MoreSystemCommands
output(polynomials...)$OutputPackage
poly1: POLYS := numer(elt(ri, x1)$SUP(FPOLYS) - y1)
poly2: POLYS := numer(elt(ri, x2)$SUP(FPOLYS) - y2)
poly3: POLYS := numer(elt(ri, x3)$SUP(FPOLYS) - y3)
-- for experimental fraction free interpolation
-- ri2: FSUPFPOLYS := map(#1::FPOLYS, numer ri)
-- $SparseUnivariatePolynomialFunctions2(POLYS, FPOLYS)
-- /map(#1::FPOLYS, denom ri) _
-- $SparseUnivariatePolynomialFunctions2(POLYS, FPOLYS)
--
-- poly1: POLYS := numer(elt(ri2, x1)$SUP(FPOLYS) - y1)
-- poly2: POLYS := numer(elt(ri2, x2)$SUP(FPOLYS) - y2)
-- poly3: POLYS := numer(elt(ri2, x3)$SUP(FPOLYS) - y3) n:Integer := len - i
o1: Integer := ord1(xValues, i)
d1: Integer := deg1(xValues, i)
o2: Integer := ord2(xValues, i)
d2: Integer := deg2(xValues, i)
-- another compiler bug: using i as iterator here makes the loop break
if debug(options)$GOPT0 then
systemCommand("sys date +%s")$MoreSystemCommands
output(interpolating resultants...)$OutputPackage
res1: SUP S
:= newton(evalResultant(poly1, poly3, o1, d1, a1, A))
$NewtonInterpolation(S)
res2: SUP S
:= newton(evalResultant(poly2, poly3, o2, d2, a1, A))
$NewtonInterpolation(S)
- if debug(options)$GOPT0 then
- res1: SUP S := univariate(resultant(poly1, poly3, a1))
-- res2: SUP S := univariate(resultant(poly2, poly3, a1))
-- if res1 ~= res1res or res2 ~= res2res then
-- output(hconcat("poly1 ", poly1::OutputForm))$OutputPackage
-- output(hconcat("poly2 ", poly2::OutputForm))$OutputPackage
-- output(hconcat("poly3 ", poly3::OutputForm))$OutputPackage
-- output(hconcat("res1 ", res1::OutputForm))$OutputPackage
-- output(hconcat("res2 ", res2::OutputForm))$OutputPackage
output("n/i: " string(n) " " string(i))$OutputPackage
output("res1 ord: " string(o1) " " string(minimumDegree res1))
$OutputPackage
output("res1 deg: " string(d1) " " string(degree res1))
$OutputPackage
output("res2 ord: " string(o2) " " string(minimumDegree res2))
$OutputPackage
output("res2 deg: " string(d2) " " string(degree res2))
$OutputPackage
if debug(options)$GOPT0 then
systemCommand("sys date +%s")$MoreSystemCommands
output(computing gcd...)$OutputPackage
-- we want to solve over F
res3: SUP F := SUPS2SUPF(primitivePart(gcd(res1, res2))) if debug(options)$GOPT0 then
systemCommand("sys date +%s")$MoreSystemCommands
output(solving...)$OutputPackage
-- res3 is a polynomial in A=a0+(len+3)*a1
-- now we have to find the roots of res3
- for f in factors factor(res3)$GF | degree f.factor = 1 repeat
- we are only interested in the linear factors
if debug(options)$GOPT0 then
output(hconcat("f: ", f::OutputForm))$OutputPackage
Av: F := -coefficient(f.factor, 0)
/ leadingCoefficient f.factor
-- FIXME: in an earlier version, we disregarded vanishing Av
-- maybe we intended to disregard vanishing a1v? Either doesn't really
-- make sense to me right now. evalPoly := eval(POLYS2POLYF poly3, A, Av)
if zero? evalPoly
then evalPoly := eval(POLYS2POLYF poly1, A, Av)
-- Note that it really may happen that poly3 vanishes when specializing
-- A. Consider for example guessExpRat([1,1,1,1]).
-- FIXME: We check poly1 below, too. I should work out in what cases poly3
-- vanishes. for g in factors factor(univariate evalPoly)$GF
| degree g.factor = 1 repeat
if debug(options)$GOPT0 then
output(hconcat("g: ", g::OutputForm))$OutputPackage a1v: F := -coefficient(g.factor, 0)
/ leadingCoefficient g.factor
-- check whether poly1 and poly2 really vanish. Note that we could have found
-- an extraneous solution, since we only computed the gcd of the two
-- resultants. t1 := eval(POLYS2POLYF poly1, [a1, A]::List V,
[a1v, Av]::List F)
if zero? t1 then
t2 := eval(POLYS2POLYF poly2, [a1, A]::List V,
[a1v, Av]::List F)
if zero? t2 then ri1: Fraction SUP POLYS
:= SUPFPOLYS2FSUPPOLYS(numer ri)
/ SUPFPOLYS2FSUPPOLYS(denom ri)
-- for experimental fraction free interpolation
-- ri1: Fraction SUP POLYS := ri numr: SUP F := SUPPOLYS2SUPF(numer ri1, a1v, Av)
denr: SUP F := SUPPOLYS2SUPF(denom ri1, a1v, Av)
if not zero? denr then
res4: EXPRR := eval(FSUPF2EXPRR(xx, numr/denr),
kernel(xx),
basis(xx::EXPRR))
p2(last xValues,
xx, a1v, Av, basis)
*xx::EXPRR
res := cons(res4, res)
else if zero? numr and debug(options)$GOPT0 then
output("numerator and denominator vanish!")
$OutputPackage
-- If we are only interested in one solution, we do not try other degrees if we
-- have found already some solutions. I.e., the indentation here is correct. if not null(res) and one(options)$GOPT0 then return res
res
guessExpRatAux0(list: List F, basis: DIFFSPECN, options: LGOPT): GUESSRESULT ==
if zero? safety(options)$GOPT0 then
error "Guess: guessExpRat does not support zero safety"
-- guesses Functions of the Form (a1n+a0)^nrat(n)
xx := indexName(options)$GOPT0
-- restrict to safety
len: Integer := #list
if len-safety(options)$GOPT0+1 < 0 then return []
shortlist: List F := first(list, (len-safety(options)$GOPT0+1)::NNI)
-- remove zeros from list
zeros: EXPRR := 1
newlist: List F
xValues: List Integer
i: Integer := -1
for x in shortlist repeat
i := i+1
if x = 0 then
zeros := zeros * (basis(xx::EXPRR) - basis(i::EXPRR))
i := -1
for x in shortlist repeat
i := i+1
if x ~= 0 then
newlist := cons(x/retract(retract(eval(zeros, xx::EXPRR,
i::EXPRR))@R),
newlist)
xValues := cons(i, xValues)
newlist := reverse newlist
xValues := reverse xValues
res: List EXPRR
:= [eval(zeros * f, xx::EXPRR, xx::EXPRR) _
for f in guessExpRatAux(xx, newlist, basis, xValues, options)]
reslist := map([#1, checkResult(#1, xx, len, list, options)], res)
$ListFunctions2(EXPRR, Record(function: EXPRR, order: NNI))
select(#1.order < len-safety(options)$GOPT0, reslist)
guessExpRat(list : List F): GUESSRESULT ==
guessExpRatAux0(list, defaultD, [])
guessExpRat(list: List F, options: LGOPT): GUESSRESULT ==
guessExpRatAux0(list, defaultD, options)
if F has RetractableTo Symbol and S has RetractableTo Symbol then
guessExpRat(q: Symbol): GUESSER ==
guessExpRatAux0(#1, q::EXPRR**#1, #2)
-- some useful types for Ore operators that work on series
-- the differentiation operator
DIFFSPECX ==> (EXPRR, Symbol, NonNegativeInteger) -> EXPRR
-- eg.: f(x)+->f(qx)
-- f(x)+->D(f, x)
DIFFSPECS ==> (UFPSF, NonNegativeInteger) -> UFPSF
-- eg.: f(x)+->f(qx)
- DIFFSPECSF ==> (UFPSSUPF, NonNegativeInteger) -> UFPSSUPF
- eg.: f(x)+->f(q*x)
-- the constant term for the inhomogeneous case
DIFFSPEC1 ==> UFPSF
DIFFSPEC1F ==> UFPSSUPF
DIFFSPEC1X ==> Symbol -> EXPRR
termAsEXPRR(f: EXPRR, xx: Symbol, l: List Integer,
DX: DIFFSPECX, D1X: DIFFSPEC1X): EXPRR ==
if empty? l then D1X(xx)
else
ll: List List Integer := powers(l)$Partition fl: List EXPRR := [DX(f, xx, (first part-1)::NonNegativeInteger)
* second(part)::NNI for part in ll]
reduce(_, fl)
termAsUFPSF(f: UFPSF, l: List Integer, DS: DIFFSPECS, D1: DIFFSPEC1): UFPSF ==
if empty? l then D1
else
ll: List List Integer := powers(l)$Partition
-- first of each element of ll is the derivative, second is the power
fl: List UFPSF := [DS(f, (first part -1)::NonNegativeInteger) _
** second(part)::NNI for part in ll]
reduce(_*, fl)
-- returns \prod f^(l.i), but using the Hadamard product
termAsUFPSF2(f: UFPSF, l: List Integer,
DS: DIFFSPECS, D1: DIFFSPEC1): UFPSF ==
if empty? l then D1
else
ll: List List Integer := powers(l)$Partition
-- first of each element of ll is the derivative, second is the power
fl: List UFPSF
:= [map(#1** second(part)::NNI, DS(f, (first part -1)::NNI)) _
for part in ll]
reduce(hadamard$UFPS1(F), fl)
termAsUFPSSUPF(f: UFPSSUPF, l: List Integer,
DSF: DIFFSPECSF, D1F: DIFFSPEC1F): UFPSSUPF ==
if empty? l then D1F
else
ll: List List Integer := powers(l)$Partition
-- first of each element of ll is the derivative, second is the power
fl: List UFPSSUPF
:= [DSF(f, (first part -1)::NonNegativeInteger)
** second(part)::NNI for part in ll]
reduce(_*, fl)
-- returns \prod f^(l.i), but using the Hadamard product
termAsUFPSSUPF2(f: UFPSSUPF, l: List Integer,
DSF: DIFFSPECSF, D1F: DIFFSPEC1F): UFPSSUPF ==
if empty? l then D1F
else
ll: List List Integer := powers(l)$Partition
-- first of each element of ll is the derivative, second is the power
fl: List UFPSSUPF
:= [map(#1 ** second(part)::NNI, DSF(f, (first part -1)::NNI)) _
for part in ll]
reduce(hadamard$UFPS1(SUP F), fl)
FilteredPartitionStream(options: LGOPT): Stream List Integer ==
maxD := 1+maxDerivative(options)$GOPT0
maxP := maxPower(options)$GOPT0 if maxD > 0 and maxP > -1 then
s := partitions(maxD, maxP)$PartitionsAndPermutations
else
s1: Stream Integer := generate(inc, 1)$Stream(Integer)
s2: Stream Stream List Integer
:= map(partitions(#1)$PartitionsAndPermutations, s1)
$StreamFunctions2(Integer, Stream List Integer)
s3: Stream List Integer
:= concat(s2)$StreamFunctions1(List Integer)
-- s := cons([],
-- select(((maxD = 0) or (first #1 <= maxD)) _
-- and ((maxP = -1) or (# #1 <= maxP)), s3)) s := cons([],
select(((maxD = 0) or (# #1 <= maxD)) _
and ((maxP = -1) or (first #1 <= maxP)), s3))
s := conjugates(s)$PartitionsAndPermutations
if homogeneous(options)$GOPT0 then rest s else s -- for functions
ADEguessStream(f: UFPSF, partitions: Stream List Integer,
DS: DIFFSPECS, D1: DIFFSPEC1): Stream UFPSF ==
map(termAsUFPSF(f, #1, DS, D1), partitions)
$StreamFunctions2(List Integer, UFPSF)
-- for coefficients, i.e., using the Hadamard product
ADEguessStream2(f: UFPSF, partitions: Stream List Integer,
DS: DIFFSPECS, D1: DIFFSPEC1): Stream UFPSF ==
map(termAsUFPSF2(f, #1, DS, D1), partitions)
$StreamFunctions2(List Integer, UFPSF)
ADEdegreeStream(partitions: Stream List Integer): Stream NNI ==
scan(0, max((if empty? #1 then 0 else (first #1 - 1)::NNI), #2),
partitions)$StreamFunctions2(List Integer, NNI)
ADEtestStream(f: UFPSSUPF, partitions: Stream List Integer,
DSF: DIFFSPECSF, D1F: DIFFSPEC1F): Stream UFPSSUPF ==
map(termAsUFPSSUPF(f, #1, DSF, D1F), partitions)
$StreamFunctions2(List Integer, UFPSSUPF)
ADEtestStream2(f: UFPSSUPF, partitions: Stream List Integer,
DSF: DIFFSPECSF, D1F: DIFFSPEC1F): Stream UFPSSUPF ==
map(termAsUFPSSUPF2(f, #1, DSF, D1F), partitions)
$StreamFunctions2(List Integer, UFPSSUPF)
ADEEXPRRStream(f: EXPRR, xx: Symbol, partitions: Stream List Integer,
DX: DIFFSPECX, D1X: DIFFSPEC1X): Stream EXPRR ==
map(termAsEXPRR(f, xx, #1, DX, D1X), partitions)
$StreamFunctions2(List Integer, EXPRR)
diffDX: DIFFSPECX
diffDX(expr, x, n) == D(expr, x, n)
diffDS: DIFFSPECS
diffDS(s, n) == D(s, n)
diffDSF: DIFFSPECSF
diffDSF(s, n) ==
-- I have to help the compiler here a little to choose the right signature...
if SUP F has _*: (NonNegativeInteger, SUP F) -> SUP F
then D(s, n)
diffAX: DIFFSPECAX
diffAX(l: NNI, x: Symbol, f: EXPRR): EXPRR ==
(x::EXPRR)*l f
diffA: DIFFSPECA
diffA(k: NNI, l: NNI, f: SUP S): S ==
DiffAction(k, l, f)$FFFG(S, SUP S)
diffAF: DIFFSPECAF
diffAF(k: NNI, l: NNI, f: UFPSSUPF): SUP F ==
DiffAction(k, l, f)$FFFG(SUP F, UFPSSUPF)
diffC: DIFFSPECC
diffC(total: NNI): List S == DiffC(total)$FFFG(S, SUP S)
diff1X: DIFFSPEC1X
diff1X(x: Symbol)== 1$EXPRR
diffHP options ==
if displayAsGF(options)$GOPT0 then
partitions := FilteredPartitionStream options
[ADEguessStream(#1, partitions, diffDS, 1$UFPSF),
ADEdegreeStream partitions,
ADEtestStream(#1, partitions, diffDSF, 1$UFPSSUPF),
ADEEXPRRStream(#1, #2, partitions, diffDX, diff1X),
diffA, diffAF, diffAX, diffC]$HPSPEC
else
error "Guess: guessADE supports only displayAsGF"
if F has RetractableTo Symbol and S has RetractableTo Symbol then
qDiffDX(q: Symbol, expr: EXPRR, x: Symbol, n: NonNegativeInteger): EXPRR ==
eval(expr, x::EXPRR, (q::EXPRR)*nx::EXPRR)
qDiffDS(q: Symbol, s: UFPSF, n: NonNegativeInteger): UFPSF ==
multiplyCoefficients((q::F)*((n#1)::NonNegativeInteger), s)
qDiffDSF(q: Symbol, s: UFPSSUPF, n: NonNegativeInteger): UFPSSUPF ==
multiplyCoefficients((q::F::SUP F)*((n#1)::NonNegativeInteger), s)
diffHP(q: Symbol): (LGOPT -> HPSPEC) ==
if displayAsGF(#1)$GOPT0 then
partitions := FilteredPartitionStream #1
[ADEguessStream(#1, partitions, qDiffDS(q, #1, #2), 1$UFPSF), _
repeating([0$NNI])$Stream(NNI),
ADEtestStream(#1, partitions, qDiffDSF(q, #1, #2), 1$UFPSSUPF),
ADEEXPRRStream(#1, #2, partitions, qDiffDX(q, #1, #2, #3), diff1X), _
diffA, diffAF, diffAX, diffC]$HPSPEC
else
error "Guess: guessADE supports only displayAsGF"
ShiftSX(expr: EXPRR, x: Symbol, n: NNI): EXPRR ==
eval(expr, x::EXPRR, x::EXPRR+n::EXPRR)
ShiftSXGF(expr: EXPRR, x: Symbol, n: NNI): EXPRR ==
if zero? n then expr
else
l := [eval(D(expr, x, i)/factorial(i)::EXPRR, x::EXPRR, 0$EXPRR)_
*(x::EXPRR)i for i in 0..n-1]
(expr-reduce(_+, l))/(x::EXPRRn)
ShiftSS(s:UFPSF, n:NNI): UFPSF ==
((quoByVar #1)**n)$MappingPackage1(UFPSF) (s)
ShiftSF(s:UFPSSUPF, n: NNI):UFPSSUPF ==
((quoByVar #1)**n)$MappingPackage1(UFPSSUPF) (s)
ShiftAX(l: NNI, n: Symbol, f: EXPRR): EXPRR ==
n::EXPRR*l f
- ShiftAXGF(l: NNI, x: Symbol, f: EXPRR): EXPRR ==
- I need to help the compiler here, unfortunately
if zero? l then f
else
s := [stirling2(l, i)$IntegerCombinatoricFunctions(Integer)::EXPRR _
(x::EXPRR)iD(f, x, i) for i in 1..l]
reduce(_+, s)
ShiftA(k: NNI, l: NNI, f: SUP S): S ==
ShiftAction(k, l, f)$FFFG(S, SUP S)
ShiftAF(k: NNI, l: NNI, f: UFPSSUPF): SUP F ==
ShiftAction(k, l, f)$FFFG(SUP F, UFPSSUPF)
ShiftC(total: NNI): List S ==
ShiftC(total)$FFFG(S, SUP S)
shiftHP options ==
partitions := FilteredPartitionStream options
if displayAsGF(options)$GOPT0 then
if maxPower(options)$GOPT0 = 1 then
[ADEguessStream(#1, partitions, ShiftSS, (1-monomial(1,1))(-1)),
ADEdegreeStream partitions,
ADEtestStream(#1, partitions, ShiftSF, (1-monomial(1,1))(-1)),
ADEEXPRRStream(#1, #2, partitions, ShiftSXGF, 1/(1-#1::EXPRR)),
ShiftA, ShiftAF, ShiftAXGF, ShiftC]$HPSPEC
else
error "Guess: no support for the Shift operator with displayAsGF _
and maxPower>1"
else
[ADEguessStream2(#1, partitions, ShiftSS, (1-monomial(1,1))(-1)),
ADEdegreeStream partitions,
ADEtestStream2(#1, partitions, ShiftSF, (1-monomial(1,1))(-1)),
ADEEXPRRStream(#1, #2, partitions, ShiftSX, diff1X),
ShiftA, ShiftAF, ShiftAX, ShiftC]$HPSPEC
if F has RetractableTo Symbol and S has RetractableTo Symbol then
qShiftAX(q: Symbol, l: NNI, n: Symbol, f: EXPRR): EXPRR ==
(q::EXPRR)*(ln::EXPRR) * f
qShiftA(q: Symbol, k: NNI, l: NNI, f: SUP S): S ==
qShiftAction(q::S, k, l, f)$FFFG(S, SUP S)
qShiftAF(q: Symbol, k: NNI, l: NNI, f: UFPSSUPF): SUP F ==
qShiftAction(q::F::SUP(F), k, l, f)$FFFG(SUP F, UFPSSUPF)
qShiftC(q: Symbol, total: NNI): List S ==
qShiftC(q::S, total)$FFFG(S, SUP S)
shiftHP(q: Symbol): (LGOPT -> HPSPEC) ==
partitions := FilteredPartitionStream #1
if displayAsGF(#1)$GOPT0 then
error "Guess: no support for the qShift operator with displayAsGF"
else
[ADEguessStream2(#1, partitions, ShiftSS, _
(1-monomial(1,1))(-1)),
ADEdegreeStream partitions,
ADEtestStream2(#1, partitions, ShiftSF, _
(1-monomial(1,1))(-1)),
ADEEXPRRStream(#1, #2, partitions, ShiftSX, diff1X),
qShiftA(q, #1, #2, #3), qShiftAF(q, #1, #2, #3), _
qShiftAX(q, #1, #2, #3), qShiftC(q, #1)]$HPSPEC
makeEXPRR(DAX: DIFFSPECAX, x: Symbol, p: SUP F, expr: EXPRR): EXPRR ==
if zero? p then 0$EXPRR
else
coerce(leadingCoefficient p)::EXPRR * DAX(degree p, x, expr) _
+ makeEXPRR(DAX, x, reductum p, expr)
list2UFPSF(list: List F): UFPSF == series(list::Stream F)$UFPSF
list2UFPSSUPF(list: List F): UFPSSUPF ==
l := [e::SUP(F) for e in list for i in 0..]::Stream SUP F
series(l)$UFPSSUPF + monomial(monomial(1,1)$SUP(F), #list)$UFPSSUPF
SUPF2SUPSUPF(p: SUP F): SUP SUP F ==
map(#1::SUP F, p)$SparseUnivariatePolynomialFunctions2(F, SUP F)
UFPSF2SUPF(f: UFPSF, deg: NNI): SUP F ==
makeSUP univariatePolynomial(f, deg)
getListSUPF(s: Stream UFPSF, o: NNI, deg: NNI): List SUP F ==
map(UFPSF2SUPF(#1, deg), entries complete first(s, o))
$ListFunctions2(UFPSF, SUP F)
S2EXPRR(s: S): EXPRR ==
if F is S then
coerce(s pretend F)@EXPRR
else if F is Fraction S then
coerce(s::Fraction(S))@EXPRR
else error "Type parameter F should be either equal to S or equal _
to Fraction S"
guessInterpolate(guessList: List SUP F, eta: List NNI, D: HPSPEC)
: Matrix SUP S ==
if F is S then
vguessList: Vector SUP S := vector(guessList pretend List(SUP(S)))
generalInterpolation((D.C)(reduce(_+, eta)), D.A,
vguessList, eta)$FFFG(S, SUP S)
else if F is Fraction S then
vguessListF: Vector SUP F := vector(guessList)
generalInterpolation((D.C)(reduce(_+, eta)), D.A,
vguessListF, eta)$FFFGF(S, SUP S, SUP F) else error "Type parameter F should be either equal to S or equal _
to Fraction S"
guessInterpolate2(guessList: List SUP F,
sumEta: NNI, maxEta: NNI,
D: HPSPEC): Stream Matrix SUP S ==
if F is S then
vguessList: Vector SUP S := vector(guessList pretend List(SUP(S)))
generalInterpolation((D.C)(sumEta), D.A,
vguessList, sumEta, maxEta)
$FFFG(S, SUP S)
else if F is Fraction S then
vguessListF: Vector SUP F := vector(guessList)
generalInterpolation((D.C)(sumEta), D.A,
vguessListF, sumEta, maxEta)
$FFFGF(S, SUP S, SUP F)
else error "Type parameter F should be either equal to S or equal _
to Fraction S"
testInterpolant(resi: List SUP S,
list: List F,
testList: List UFPSSUPF,
exprList: List EXPRR,
initials: List EXPRR,
guessDegree: NNI,
D: HPSPEC,
dummy: Symbol, op: BasicOperator, options: LGOPT,
list: List F)
: Union("failed", Record(function: EXPRR, order: NNI)) ==
((maxDegree(options)$GOPT0 = -1) and
(allDegrees(options)$GOPT0 = false) and
zero?(last resi))
=> return "failed"
nonZeroCoefficient: Integer := 0
for i in 1..#resi repeat
if not zero? resi.i then
if zero? nonZeroCoefficient then
nonZeroCoefficient := i
else
nonZeroCoefficient := 0
break
if not zero? nonZeroCoefficient then
(freeOf?(exprList.nonZeroCoefficient, name op)) => return "failed" for e in list repeat
if not zero? e then return "failed"
else
resiSUPF := map(SUPF2SUPSUPF SUPS2SUPF #1, resi)
$ListFunctions2(SUP S, SUP SUP F)
iterate? := true;
for d in guessDegree+1.. repeat
c: SUP F := generalCoefficient(D.AF, vector testList,
d, vector resiSUPF)
$FFFG(SUP F, UFPSSUPF) if not zero? c then
iterate? := ground? c
break
iterate? => return "failed"
g: SUP S := gcd resi
resiF := map(SUPS2SUPF((#1 exquo g)::SUP(S)), resi)
$ListFunctions2(SUP S, SUP F)
if debug(options)$GOPT0 then
output(hconcat("trying possible solution ", resiF::OutputForm))
$OutputPackage
-- transform each term into an expression
ex: List EXPRR := [makeEXPRR(D.AX, dummy, p, e) _
for p in resiF for e in exprList]
-- transform the list of expressions into a sum of expressions
res: EXPRR
if displayAsGF(options)$GOPT0 then
res := evalADE(op, dummy, variableName(options)$GOPT0::EXPRR,
indexName(options)$GOPT0::EXPRR,
numerator reduce(_+, ex),
rev |