\documentclass{article} \usepackage{axiom,amsthm,amsmath} \newtheorem{ToDo}{ToDo}[section] \newcommand{\Axiom}{{\tt Axiom}} \newcommand{\Rate}{{\tt Rate}} \newcommand{\GFUN}{{\tt GFUN}} \begin{document} \title{mantepse.spad} \author{Martin Rubey} \maketitle \begin{abstract} The packages defined in this file enable {\Axiom} to guess formulas for sequences of, for example, rational numbers or rational functions, given the first few terms. It extends and complements Christian Krattenthaler's program \Rate\ and the relevant parts of Bruno Salvy and Paul Zimmermann's \GFUN. \end{abstract} \tableofcontents \section{domain UFPS UnivariateFormalPowerSeries} <>= )abbrev domain UFPS UnivariateFormalPowerSeries UnivariateFormalPowerSeries(Coef: Ring) == UnivariateTaylorSeries(Coef, 'x, 0$Coef) @ %$ \section{package UFPS1 UnivariateFormalPowerSeriesFunctions} <>= )abbrev package UFPS1 UnivariateFormalPowerSeriesFunctions UnivariateFormalPowerSeriesFunctions(Coef: Ring): Exports == Implementation where UFPS ==> UnivariateFormalPowerSeries Coef Exports == with hadamard: (UFPS, UFPS) -> UFPS Implementation == add hadamard(f, g) == series map(#1*#2, coefficients f, coefficients g) $StreamFunctions3(Coef, Coef, Coef) @ %$ \section{domain GOPT GuessOption} <>= )abbrev domain GOPT GuessOption ++ Author: Martin Rubey ++ Description: GuessOption is a domain whose elements are various options used ++ by \spadtype{Guess}. GuessOption(): Exports == Implementation where Exports == SetCategory with maxDerivative: Integer -> % ++ maxDerivative(d) specifies the maximum derivative in an algebraic ++ differential equation. maxDerivative(-1) specifies that the maximum ++ derivative can be arbitrary. This option is expressed in the form ++ \spad{maxDerivative == d}. maxShift: Integer -> % ++ maxShift(d) specifies the maximum shift in a recurrence ++ equation. maxShift(-1) specifies that the maximum shift can be ++ arbitrary. This option is expressed in the form \spad{maxShift == d}. maxPower: Integer -> % ++ maxPower(d) specifies the maximum degree in an algebraic differential ++ equation. For example, the degree of (f'')^3 f' is 4. maxPower(-1) ++ specifies that the maximum exponent can be arbitrary. This option is ++ expressed in the form \spad{maxPower == d}. homogeneous: Boolean -> % ++ homogeneous(d) specifies whether we allow only homogeneous algebraic ++ differential equations. This option is expressed in the form ++ \spad{homogeneous == d}. maxLevel: Integer -> % ++ maxLevel(d) specifies the maximum number of recursion levels operators ++ guessProduct and guessSum will be applied. maxLevel(-1) specifies ++ that all levels are tried. This option is expressed in the form ++ \spad{maxLevel == d}. maxDegree: Integer -> % ++ maxDegree(d) specifies the maximum degree of the coefficient ++ polynomials in an algebraic differential equation or a recursion with ++ polynomial coefficients. For rational functions with an exponential ++ term, \spad{maxDegree} bounds the degree of the denominator ++ polynomial. ++ maxDegree(-1) specifies that the maximum ++ degree can be arbitrary. This option is expressed in the form ++ \spad{maxDegree == d}. allDegrees: Boolean -> % ++ allDegrees(d) specifies whether all possibilities of the degree vector ++ - taking into account maxDegree - should be tried. This is mainly ++ interesting for rational interpolation. This option is expressed in ++ the form \spad{allDegrees == d}. safety: NonNegativeInteger -> % ++ safety(d) specifies the number of values reserved for testing any ++ solutions found. This option is expressed in the form \spad{safety == ++ d}. one: Boolean -> % ++ one(d) specifies whether we are happy with one solution. This option ++ is expressed in the form \spad{one == d}. debug: Boolean -> % ++ debug(d) specifies whether we want additional output on the ++ progress. This option is expressed in the form \spad{debug == d}. functionName: Symbol -> % ++ functionName(d) specifies the name of the function given by the ++ algebraic differential equation or recurrence. This option is ++ expressed in the form \spad{functionName == d}. variableName: Symbol -> % ++ variableName(d) specifies the variable used in by the algebraic ++ differential equation. This option is expressed in the form ++ \spad{variableName == d}. indexName: Symbol -> % ++ indexName(d) specifies the index variable used for the formulas. This ++ option is expressed in the form \spad{indexName == d}. displayAsGF: Boolean -> % ++ displayAsGF(d) specifies whether the result is a generating function ++ or a recurrence. This option should not be set by the user, but rather ++ by the HP-specification. option : (List %, Symbol) -> Union(Any, "failed") ++ option() is not to be used at the top level; ++ option determines internally which drawing options are indicated in ++ a draw command. option?: (List %, Symbol) -> Boolean ++ option?() is not to be used at the top level; ++ option? internally returns true for drawing options which are ++ indicated in a draw command, or false for those which are not. checkOptions: List % -> Void ++ checkOptions checks whether an option is given twice Implementation ==> add import AnyFunctions1(Boolean) import AnyFunctions1(Symbol) import AnyFunctions1(Integer) import AnyFunctions1(NonNegativeInteger) Rep := Record(keyword: Symbol, value: Any) maxLevel d == ["maxLevel"::Symbol, d::Any] maxDerivative d == ["maxDerivative"::Symbol, d::Any] maxShift d == maxDerivative d maxDegree d == ["maxDegree"::Symbol, d::Any] allDegrees d == ["allDegrees"::Symbol, d::Any] maxPower d == ["maxPower"::Symbol, d::Any] safety d == ["safety"::Symbol, d::Any] homogeneous d == ["homogeneous"::Symbol, d::Any] debug d == ["debug"::Symbol, d::Any] one d == ["one"::Symbol, d::Any] functionName d == ["functionName"::Symbol, d::Any] variableName d == ["variableName"::Symbol, d::Any] indexName d == ["indexName"::Symbol, d::Any] displayAsGF d == ["displayAsGF"::Symbol, d::Any] coerce(x:%):OutputForm == x.keyword::OutputForm = x.value::OutputForm x:% = y:% == x.keyword = y.keyword and x.value = y.value option?(l, s) == for x in l repeat x.keyword = s => return true false option(l, s) == for x in l repeat x.keyword = s => return(x.value) "failed" checkOptions l == if not empty? l then if find((first l).keyword = #1.keyword, rest l) case "failed" then checkOptions rest l else error "GuessOption: Option specified twice" @ \section{package GOPT0 GuessOptionFunctions0} <>= )abbrev package GOPT0 GuessOptionFunctions0 ++ Author: Martin Rubey ++ Description: GuessOptionFunctions0 provides operations that extract the ++ values of options for \spadtype{Guess}. GuessOptionFunctions0(): Exports == Implementation where LGOPT ==> List GuessOption Exports == SetCategory with maxLevel: LGOPT -> Integer ++ maxLevel returns the specified maxLevel or -1 as default. maxPower: LGOPT -> Integer ++ maxPower returns the specified maxPower or -1 as default. maxDerivative: LGOPT -> Integer ++ maxDerivative returns the specified maxDerivative or -1 as default. maxShift: LGOPT -> Integer ++ maxShift returns the specified maxShift or -1 as default. maxDegree: LGOPT -> Integer ++ maxDegree returns the specified maxDegree or -1 as default. allDegrees: LGOPT -> Boolean ++ allDegrees returns whether all possibilities of the degree vector ++ should be tried, the default being false. safety: LGOPT -> NonNegativeInteger ++ safety returns the specified safety or 1 as default. one: LGOPT -> Boolean ++ one returns whether we need only one solution, default being true. homogeneous: LGOPT -> Boolean ++ homogeneous returns whether we allow only homogeneous algebraic ++ differential equations, default being false functionName: LGOPT -> Symbol ++ functionName returns the name of the function given by the algebraic ++ differential equation, default being f variableName: LGOPT -> Symbol ++ variableName returns the name of the variable used in by the ++ algebraic differential equation, default being x indexName: LGOPT -> Symbol ++ indexName returns the name of the index variable used for the ++ formulas, default being n displayAsGF: LGOPT -> Boolean ++ displayAsGF specifies whether the result is a generating function ++ or a recurrence. This option should not be set by the user, but rather ++ by the HP-specification, therefore, there is no default. debug: LGOPT -> Boolean ++ debug returns whether we want additional output on the progress, ++ default being false Implementation == add maxLevel l == if (opt := option(l, "maxLevel" :: Symbol)) case "failed" then -1 else retract(opt :: Any)$AnyFunctions1(Integer) maxDerivative l == if (opt := option(l, "maxDerivative" :: Symbol)) case "failed" then -1 else retract(opt :: Any)$AnyFunctions1(Integer) maxShift l == maxDerivative l maxDegree l == if (opt := option(l, "maxDegree" :: Symbol)) case "failed" then -1 else retract(opt :: Any)$AnyFunctions1(Integer) allDegrees l == if (opt := option(l, "allDegrees" :: Symbol)) case "failed" then false else retract(opt :: Any)$AnyFunctions1(Boolean) maxPower l == if (opt := option(l, "maxPower" :: Symbol)) case "failed" then -1 else retract(opt :: Any)$AnyFunctions1(Integer) safety l == if (opt := option(l, "safety" :: Symbol)) case "failed" then 1$NonNegativeInteger else retract(opt :: Any)$AnyFunctions1(Integer)::NonNegativeInteger one l == if (opt := option(l, "one" :: Symbol)) case "failed" then true else retract(opt :: Any)$AnyFunctions1(Boolean) debug l == if (opt := option(l, "debug" :: Symbol)) case "failed" then false else retract(opt :: Any)$AnyFunctions1(Boolean) homogeneous l == if (opt := option(l, "homogeneous" :: Symbol)) case "failed" then false else retract(opt :: Any)$AnyFunctions1(Boolean) variableName l == if (opt := option(l, "variableName" :: Symbol)) case "failed" then "x" :: Symbol else retract(opt :: Any)$AnyFunctions1(Symbol) functionName l == if (opt := option(l, "functionName" :: Symbol)) case "failed" then "f" :: Symbol else retract(opt :: Any)$AnyFunctions1(Symbol) indexName l == if (opt := option(l, "indexName" :: Symbol)) case "failed" then "n" :: Symbol else retract(opt :: Any)$AnyFunctions1(Symbol) displayAsGF l == if (opt := option(l, "displayAsGF" :: Symbol)) case "failed" then error "GuessOption: displayAsGF not set" else retract(opt :: Any)$AnyFunctions1(Boolean) @ \section{package GUESS Guess} <>= )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 @ The following should be commented out when not debugging, see also Section~\ref{sec:Hermite-Pade}. <>= --@<> @ <>= 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. --@<> @ <>= termAsUFPSF: (UFPSF, List Integer, DIFFSPECS, DIFFSPEC1) -> UFPSF termAsUFPSF2: (UFPSF, List Integer, DIFFSPECS, DIFFSPEC1) -> UFPSF termAsEXPRR: (EXPRR, Symbol, List Integer, DIFFSPECX, DIFFSPEC1X) -> EXPRR termAsUFPSSUPF: (UFPSSUPF, List Integer, DIFFSPECSF, DIFFSPEC1F) -> UFPSSUPF termAsUFPSSUPF2: (UFPSSUPF, List Integer, DIFFSPECSF, DIFFSPEC1F) -> UFPSSUPF ShiftSXGF: (EXPRR, Symbol, NNI) -> EXPRR ShiftAXGF: (NNI, Symbol, EXPRR) -> EXPRR FilteredPartitionStream: LGOPT -> Stream List Integer guessInterpolate: (List SUP F, List NNI, HPSPEC) -> Matrix SUP S testInterpolant: (List SUP S, List F, List UFPSSUPF, List EXPRR, List EXPRR, _ NNI, HPSPEC, Symbol, BasicOperator, LGOPT, List F) _ -> Union("failed", Record(function: EXPRR, order: NNI)) @ <>= Implementation == add <> @ <>= -- We have to put this chunk at the beginning, because otherwise it will take -- very long to compile. <> <> <> <> <> @ \subsection{general utilities} <>= 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" @ \subsection{guessing rational functions with an exponential term} <>= <> <> @ \subsubsection{Utilities} \paragraph{conversion routines} <>= 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) @ \paragraph{mimicking $q$-anaologa} <>= 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) <> <> @ \subsubsection{reducing the degree of the polynomials} The degree of [[poly3]] is governed by $(a0+x_m a_1)^{x_m}$. Therefore, we substitute $A-x_m a1$ for $a$, which reduces the degree in $a_1$ by $x_m-x_{i+1}$. <>= 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)) @ <>= p(xm: Integer, i: Integer, va1: V, vA: V, basis: DIFFSPECN): FPOLYS == vA::POLYS::FPOLYS + (i-xm)*va1::POLYS::FPOLYS p2(xm: Integer, i: Symbol, a1v: F, Av: F, basis: DIFFSPECN): EXPRR == coerce(Av) + coerce(a1v)*(i::EXPRR - xm::EXPRR) @ \subsubsection{Order and Degree} The following expressions for order and degree of the resultants [[res1]] and [[res2]] in [[guessExpRatAux]] were first guessed -- partially with the aid of [[guessRat]], and then proven to be correct. <>= 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) + _ 2*reduce(_+, [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)) @ [[evalResultant]] evaluates the resultant of [[p1]] and [[p2]] at [[d-o+1]] points, so that we can recover it by interpolation. <>= 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) @ It may happen, that the leading coefficients of one or both of the polynomials changes, when we evaluate it at $k$. In this case, we need to correct this by multiplying with the corresponding power of the leading coefficient of the other polynomial. Consider the Sylvester matrix of the original polynomials. We want to evaluate it at $A=k$. If the first few leading coefficients of $p2$ vanish, the first few columns of the Sylvester matrix have triangular shape, with the leading coefficient of $p1$ on the diagonal. The same thing happens, if we exchange the roles of $p1$ and $p2$, only that we have to take care of the sign, too. <>= 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) @ %$ Since we also have an lower bound for the order of the resultant, we need to evaluate it only at $d-o+1$ points. Furthermore, we can divide by $k^o$ and still obtain a polynomial. <>= reverse res @ \subsubsection{The main routine} <>= 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)) @ We try to fit the data $(s1,s2,\dots)$ to the model $(a+b n)^n y(n)$, $r$ being a rational function. To obtain $y$, we compute $y(n)=s_n*(a+b n)^{-n}$. <>= 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 "::OutputForm, i::OutputForm)) $OutputPackage @ \begin{ToDo} Shouldn't we use the algorithm over [[POLYS]] here? Strange enough, for polynomial interpolation, it is faster, but for rational interpolation \emph{much} slower. This should be investigated. \end{ToDo} \begin{ToDo} It seems that [[maxDeg]] bounds the degree of the denominator, rather than the numerator? This is now added to the documentation of [[maxDegree]], it should make sense. \end{ToDo} <>= if debug(options)$GOPT0 then systemCommand("sys date +%s")$MoreSystemCommands output("interpolating..."::OutputForm)$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..."::OutputForm)$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..."::OutputForm)$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..."::OutputForm)$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..."::OutputForm)$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 (a1*n+a0)^n*rat(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) @ \subsection{Hermite Pad\'e interpolation}\label{sec:Hermite-Pade} <>= <> <> <> <> <> <> @ \subsubsection{Types for Operators} <>= -- some useful types for Ore operators that work on series -- the differentiation operator DIFFSPECX ==> (EXPRR, Symbol, NonNegativeInteger) -> EXPRR -- eg.: f(x)+->f(q*x) -- f(x)+->D(f, x) DIFFSPECS ==> (UFPSF, NonNegativeInteger) -> UFPSF -- eg.: f(x)+->f(q*x) DIFFSPECSF ==> (UFPSSUPF, NonNegativeInteger) -> UFPSSUPF -- eg.: f(x)+->f(q*x) -- the constant term for the inhomogeneous case DIFFSPEC1 ==> UFPSF DIFFSPEC1F ==> UFPSSUPF DIFFSPEC1X ==> Symbol -> EXPRR @ \subsubsection{Streams}\label{sec:streams} In this section we define some functions that provide streams for [[HermitePade]]. The following three functions transform a partition [[l]] into a product of derivatives of [[f]], using the given operators. We need to provide the same functionality for expressions, series and series with a transcendental element. Only for expressions we do not provide a version using the Hadamard product, although it would be quite easy to define an appropriate operator on expressions. A partition $(\lambda_1^{p_1},\lambda_2^{p_2},\dots)$ is transformed into the expression $(f^{(\lambda_1-1)})^{p_1}(f^{(\lambda_2-1)})^{p_2}\cdots$, i.e., the size of the part is interpreted as derivative, the exponent as power. <>= 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) @ %$ It is not clear whether we should \lq prefer\rq\ shifting and differentiation over powering. Currently, we produce the stream \begin{equation*} \begin{array}{rrrrrrrrr} \emptyset& 1& 11 & 2 & 111& 2 1 & 3 & 1111\\ 1& f& f^2& f'& f^3& f f'& f''& f^4 &\dots \end{array} \end{equation*} Maybe it would be better to produce \begin{equation*} \begin{array}{rrrrrrrrr} \emptyset& 1& 2& 11 & 3 & 21 & 111& 4\\ 1& f& f'& f^2& f''& f f'& f^3& f''' &\dots \end{array} \end{equation*} instead, i.e., to leave the partitions unconjugated. Note however, that shifting and differentiation decrease the number of valid terms, while powering does not. Note that we conjugate all partitions at the very end of the following procedure\dots <>= 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) @ %$ The entries of the following stream indicate how many terms we loose when applying one of the power and shift or differentiation operators. More precisely, the $n$\textsuperscript{th} entry of the stream takes into account all partitions up to index $n$. Thus, the entries of the stream are weakly increasing. <>= 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) @ \subsubsection{Operators} We need to define operators that transform series for differentiation and shifting. We also provide operators for $q$-analoga. The functionality corresponding to powering and taking the Hadamard product if provided by the streams, see Section~\ref{sec:streams}. We have to provide each operator in three versions: \begin{itemize} \item for expressions, \item for series, and \item for series with an additional transcendental element. \end{itemize} The latter makes it possible to detect lazily whether a computed coefficient of a series is valid or not. Furthermore, we have to define for each operator how to extract the coefficient of $x^k$ in $z^l f(x)$, where multiplication with $z$ is defined depending on the operator. Again, it is necessary to provide this functionality for expressions, series and series with a transcendental element. Finally, we define a function that returns the diagonal elements $c_{k,k}$ in the expansion $\langle x^k\rangle z f(x) = \sum_{i=0}^k c_{k,i} \langle x^i\rangle f(x)$, and an expression that represents the constant term for the inhomogeneous case. \paragraph{The Differentiation Setting} In this setting, we have $z f(x) := x f(x)$. <>= 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) @ The next three functions extract the coefficient of $x^k$ in $z^l f(x)$. Only, for expressions, we rather need $\sum_{k\ge0} \langle x^k\rangle z^l f(x)$, i.e., the function itself, which is by definition equal to $x^l f(x)$. <>= 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" @ \paragraph{$q$-dilation} In this setting, we also have $z f(x) := x f(x)$, therefore we can reuse some of the functions of the previous paragraph. Differentiation is defined by $D_q f(x, q) = f(qx, q)$. <>= 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)**n*x::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" @ \paragraph{Shifting} The shift operator transforms a sequence $u(k)$ into $u(k+1)$. We also provide operators [[ShiftSXGF]], [[ShiftAXGF]] that act on the power series, as long as no powering is involved. In this case, shifting transforms $f(x)$ into $\frac{f(x)-f(0)}{x}$. Multiplication with $z$ transforms the coefficients $u(n)$ of the series into $z u(n) := n u(n)$. The description in terms of power series is given by $xDf(x)$. % The coefficients of $x^n$ are $1, f(n), f(n+1), f(n)^2, f(n)f(n+1),\dots$ % What does this remark mean? <>= 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::EXPRR**n) ShiftSS(s:UFPSF, n:NNI): UFPSF == ((quoByVar #1)**n)$MappingPackage1(UFPSF) (s) ShiftSF(s:UFPSSUPF, n: NNI):UFPSSUPF == ((quoByVar #1)**n)$MappingPackage1(UFPSSUPF) (s) @ %$ As before, next three functions extract the coefficient of $x^k$ in $z^l f(x)$. <>= 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)**i*D(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 @ \paragraph{$q$-Shifting} The $q$-shift also transforms $u(n)$ into $u(n+1)$, and we can reuse the corrseponding functions of the previous paragraph. However, this time multiplication with $z$ is defined differently: the coefficient of $x^k$ in $z u(n)$ is $q^n u(n)$. We do not define the corresponding functionality for power series. %The coefficients of $x^n$ are $1, f(n), f(n+1), f(n)^2, f(n)f(n+1),\dots$ % What does this remark mean? <>= if F has RetractableTo Symbol and S has RetractableTo Symbol then qShiftAX(q: Symbol, l: NNI, n: Symbol, f: EXPRR): EXPRR == (q::EXPRR)**(l*n::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 @ %$ \paragraph{Extend the action to polynomials} The following operation uses the given action of $z$ on a function to multiply a $f$ with a polynomial. <>= 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) @ %$ \subsubsection{Utilities} [[list2UFPSF]] and [[list2UFPSSUPF]] transform the list passed to the guessing functions into a series. One might be tempted to transform the list into a polynomial instead, but the present approach makes computing powers and derivatives much cheaper, since, because of laziness, only coefficients that are actually used are computed. The second of the two procedures augments the list with a transcendental element. This way we can easily check whether a coefficient is valid or not: if it contains the transcendental element, it depends on data we do not have. In other words, this transcendental element simply represents the $O(x^d)$, when $d$ is the number of elements in the list. <>= 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]] interprets each coefficient as a univariate polynomial. <>= SUPF2SUPSUPF(p: SUP F): SUP SUP F == map(#1::SUP F, p)$SparseUnivariatePolynomialFunctions2(F, SUP F) @ [[getListSUPF]] returns the first [[o]] elements of the stream, each truncated after degree [[deg]]. <>= 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]] calls the appropriate [[generalInterpolation]] from [[FFFG]], for one vector of degrees, namely [[eta]]. <>= 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]] calls the appropriate [[generalInterpolation]] from [[FFFG]], for all degree vectors with given [[sumEta]] and [[maxEta]]. <>= 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]] checks whether p is really a solution. <>= 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)) == @ First we make sure it is not a solution we should have found already. Note that we cannot check this if [[maxDegree]] is set, in which case some initial solutions may have been overlooked. <>= ((maxDegree(options)$GOPT0 = -1) and (allDegrees(options)$GOPT0 = false) and zero?(last resi)) => return "failed" @ Then we check all the coefficients that should be valid. We want the zero solution only, if the function is really zero. Without this test, every sequence ending with zero is interpreted as the zero sequence, since the factor in front of the only non-vanishing term can cancel everything else. <>= 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 @ We set [[nonZeroCoefficient]] to the only non zero coefficient or, if there are several, it is $0$. It should not happen that all coefficients in [[resi]] vanish. \begin{ToDo} Check that not all coefficients in [[resi]] can vanish simultaneously. \end{ToDo} <>= if not zero? nonZeroCoefficient then (freeOf?(exprList.nonZeroCoefficient, name op)) => return "failed" for e in list repeat if not zero? e then return "failed" @ We first deal with the case that there is only one non-vanishing coefficient in [[resi]]. If the only expression in [[exprList]] whose coefficient does not vanish does not contain the name of the generating function, or if there is a non-zero term in [[list]], we reject the proposed solution. <>= 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" @ Here we check for each degree [[d]] larger than [[guessDegree]] whether the proposed linear combination vanishes. When the first coefficient that does not vanish contains the transcendental element we accept the linear combination as a solution. Finally, it seems that we have found a solution. Now we cancel the greatest common divisor of the equation. Note that this may take quite some time, it seems to be quicker to check [[generalCoefficient]] with the original [[resi]]. <>= 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) @ \subsubsection{The main routine} The following is the main guessing routine, called by all others -- except [[guessExpRat]]. <>= guessHPaux(list: List F, D: HPSPEC, options: LGOPT): GUESSRESULT == reslist: GUESSRESULT := [] listDegree := #list-1-safety(options)$GOPT0 if listDegree < 0 then return reslist @ %$ \sloppypar We do as if we knew only the coefficients up to and including degree [[listDegree-safety]]. Thus, if we have less elements of the sequence than [[safety]], there remain no elements for guessing. Originally we demanded to have at least two elements for guessing. However, we want to be able to induce from $[c,c]$ that the function equals $c$ with [[safety]] one. This is important if we apply, for example, [[guessRat]] recursively. In this case, [[listDegree]] equals zero. <>= 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] @ We need to create several streams. Let $P$ be the univariate power series whose first few elements are given by [[list]]. As an example, consider the differentiation setting. In this case, the elements of [[guessS]] consist of $P$ differentiated and taken to some power. The elements of [[degreeS]] are integers, that tell us how many terms less than in [[list]] are valid in the corresponding element of [[guessS]]. The elements of [[testS]] are very similar to those of [[guessS]], with the difference that they are derived from $P$ with an transcendental element added, which corresponds to $O(x^d)$. Finally, the elements of [[exprS]] contain representations of the transformations applied to $P$ as expressions. \begin{ToDo} I am not sure whether it is better to get rid of denominators in [[list]] here or, as I currently do it, only in [[generalInterpolation$FFFG]]. %$ If we clear them at here, we cannot take advantage of factors that may appear only after differentiating or powering. \end{ToDo} <>= guessS := (D.guessStream)(list2UFPSF list) degreeS := D.degreeStream testS := (D.testStream)(list2UFPSSUPF list) exprS := (D.exprStream)(op(dummy::EXPRR)::EXPRR, dummy) @ We call the number of terms of the linear combination its \emph{order}. We consider linear combinations of at least two terms, since otherwise the function would have to vanish identically\dots When proceeding in the stream [[guessS]], it may happen that we loose valid terms. For example, when trying to find a linear ordinary differential equation, we are looking for a linear combination of $f, f^\prime, f^{\prime\prime}, \dots$. Suppose [[listDegree]] equals $2$, i.e. we have \begin{align*} f &= l_0 + l_1 x + l_2 x^2\\ f^\prime &= l_1 + 2l_2 x\\ f^{\prime\prime} &= 2l_2. \end{align*} Thus, each time we differentiate the series, we loose one term for guessing. Therefore, we cannot use the coefficient of $x^2$ of $f^{\prime\prime}$ for generating a linear combination. [[guessDegree]] contains the degree up to which all the generated series are correct, taking into account [[safety]]. <>= 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 @ %$ We now have to distinguish between the case where we try all combination of degrees and the case where we try only an (nearly) evenly distributed vector of degrees. In the first case, [[guessInterpolate2]] is going to look at all degree vectors with [[o]] elements with total degree [[guessDegree+1]]. We give up as soon as the order [[o]] is greater than the number of available terms plus one. In the extreme case, i.e., when [[o=guessDegree+2]], we allow for example constant coefficients for all terms of the linear combination. It seems that it is a bit arbitrary at what point we stop, however, we would like to be consistent with the evenly distributed case. <>= if allDegrees(options)$GOPT0 then (o > guessDegree+2) => return reslist maxEta: Integer := 1+maxDegree(options)$GOPT0 if maxEta = 0 then maxEta := guessDegree+1 else @ In the second case, we first compute the number of parameters available for determining the coefficient polynomials. We have to take care of the fact, that HermitePade produces solutions with sum of degrees being one more than the sum of elements in [[eta]]. <>= maxParams := divide(guessDegree::NNI+1, o) if debug(options)$GOPT0 then output(hconcat("maxParams: ", maxParams::OutputForm)) $OutputPackage @ If we do not have enough parameters, we give up. We allow for at most one zero entry in the degree vector, because then there is one column that allows at least a constant term in each entry. <>= if maxParams.quotient = 0 and maxParams.remainder < o-1 then return reslist @ If [[maxDegree]] is set, we skip the first few partitions, unless we cannot increase the order anymore. \begin{ToDo} I have no satisfactory way to determine whether we can increase the order or not. \end{ToDo} <>= 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] @ We distribute the parameters as evenly as possible. Is it better to have higher degrees at the end or at the beginning? It remains to prepare [[guessList]], which is the list of [[o]] power series polynomials truncated after degree [[guessDegree]]. Then we can call HermitePade. \begin{ToDo} [[maxDegree]] should be handled differently, maybe: we should only pass as many coefficients to [[FFFG]] as [[maxDegree]] implies! \end{ToDo} <>= 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 @ \begin{ToDo} The duplicated block at the end should really go into [[testInterpolant]], I guess. Furthermore, it would be better to remove duplicates already there. \end{ToDo} \subsubsection{Specialisations} For convenience we provide some specialisations that cover the most common situations. \paragraph{generating functions} <>= 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, []) @ \paragraph{$q$-generating functions} <>= 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) @ %$ \paragraph{coefficients} <>= 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, []) @ %$ \paragraph{$q$-coefficients} <>= 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) @ %$ \subsection{[[guess]] -- applying operators recursively} The main observation made by Christian Krattenthaler in designing his program \Rate\ is the following: it occurs frequently that although a sequence of numbers is not generated by a rational function, the sequence of successive quotients is. We slightly extend upon this idea, and apply recursively one or both of the two following operators: \begin{description} \item[$\Delta_n$] the differencing operator, transforming $f(n)$ into $f(n)-f(n-1)$, and \item[$Q_n$] the operator that transforms $f(n)$ into $f(n)/f(n-1)$. \end{description} <>= 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 "::OutputForm, 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, []) @ %$ \section{package GUESSINT GuessInteger} <>= -- 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)) @ \section{package GUESSF1 GuessFiniteFunctions} <>= )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 @ \section{package GUESSF GuessFinite} <>= )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)) @ \section{package GUESSP GuessPolynomial} <>= )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)) @ \section{License} <>= -- Copyright (C) 2006 Martin Rubey -- -- This program is free software; you can redistribute it and/or -- modify it under the terms of the GNU General Public License as -- published by the Free Software Foundation; either version 2 of -- the License, or (at your option) any later version. -- -- This program is distributed in the hope that it will be -- useful, but WITHOUT ANY WARRANTY; without even the implied -- warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR -- PURPOSE. See the GNU General Public License for more details. @ <<*>>= <> <> <> <> <> <> <> <> <> <> @ \end{document}