diff --git a/books/bookvol10.3.pamphlet b/books/bookvol10.3.pamphlet index 6bf7891..114599a 100644 --- a/books/bookvol10.3.pamphlet +++ b/books/bookvol10.3.pamphlet @@ -12594,7 +12594,7 @@ string(dX::Symbol) -- succeeds --R --R (23) [dA ,dA ] --R 1 2 ---R Type: List Union(BasicStochasticDifferential,"failed") +--R Type: List(Union(BasicStochasticDifferential,"failed")) --E 27 --S 28 of 30 @@ -12619,7 +12619,7 @@ print copyIto() --R --R (26) [dA ,dA ] --R 1 2 ---R Type: List Union(BasicStochasticDifferential,Integer) +--R Type: List(Union(BasicStochasticDifferential,Integer)) --E 30 @@ -12768,10 +12768,10 @@ o )show BasicStochasticDifferential ++ We allow a separate function convertIfCan which will check whether the ++ argument has previously been declared as a BSD. -BasicStochasticDifferential(): Category == Implementation where +BasicStochasticDifferential(): category == implementation where INT ==> Integer OF ==> OutputForm - Category ==> OrderedSet with + category ==> OrderedSet with ConvertibleTo(Symbol) convertIfCan: Symbol -> Union(%, "failed") @@ -12815,7 +12815,7 @@ BasicStochasticDifferential(): Category == Implementation where ++X introduce!(t,dt) -- dt is a new stochastic differential ++X getSmgl(dt::BSD) - Implementation ==> Symbol add + implementation ==> Symbol add Rep := Symbol @@ -49051,7 +49051,7 @@ FreeNilpotentLie(n:NNI,class:NNI,R: CommutativeRing): Export == Implement where Fx := FRAC UP(x, FRAC INT) --R --R ---R (1) Fraction UnivariatePolynomial(x,Fraction Integer) +--R (1) Fraction(UnivariatePolynomial(x,Fraction Integer)) --R Type: Domain --E 1 @@ -53616,6 +53616,7 @@ heapsort(x) == (empty? x => []; cons(extract!(x),heapsort x)) h1 := heapsort heap [17,-4,9,-11,2,7,-7] --R --R Compiling function heapsort with type Heap(Integer) -> List(Integer) +--R --R --R (12) [17,9,7,2,- 4,- 7,- 11] --R Type: List Integer diff --git a/books/bookvol10.4.pamphlet b/books/bookvol10.4.pamphlet index 79517af..2783623 100644 --- a/books/bookvol10.4.pamphlet +++ b/books/bookvol10.4.pamphlet @@ -4067,14 +4067,16 @@ AnyFunctions1(S:Type): with getDomains 'Collection --R --R (1) ---R {AssociationList, Bits, CharacterClass, DataList, EqTable, FlexibleArray, ---R GeneralPolynomialSet, GeneralSparseTable, GeneralTriangularSet, HashTable, ---R IndexedBits, IndexedFlexibleArray, IndexedList, IndexedOneDimensionalArray, +--R {AssociationList, Bits, CharacterClass, ComplexDoubleFloatVector, DataList, +--R DoubleFloatVector, EqTable, FlexibleArray, GeneralPolynomialSet, +--R GeneralSparseTable, GeneralTriangularSet, HashTable, IndexedBits, +--R IndexedFlexibleArray, IndexedList, IndexedOneDimensionalArray, --R IndexedString, IndexedVector, InnerTable, KeyedAccessFile, Library, List, ---R ListMultiDictionary, Multiset, OneDimensionalArray, Point, PrimitiveArray, ---R RegularChain, RegularTriangularSet, Result, RoutinesTable, Set, ---R SparseTable, SquareFreeRegularTriangularSet, Stream, String, StringTable, ---R Table, Vector, WuWenTsunTriangularSet} +--R ListMultiDictionary, Multiset, NeitherSparseOrDensePowerSeries, +--R OneDimensionalArray, Point, PrimitiveArray, RegularChain, +--R RegularTriangularSet, Result, RoutinesTable, Set, SparseTable, +--R SquareFreeRegularTriangularSet, Stream, String, StringTable, Table, +--R U32Vector, Vector, WuWenTsunTriangularSet} --R Type: Set(Symbol) --E 1 @@ -4231,14 +4233,16 @@ a set of domains which inherit from that category: getDomains 'Collection - {AssociationList, Bits, CharacterClass, DataList, EqTable, FlexibleArray, - GeneralPolynomialSet, GeneralSparseTable, GeneralTriangularSet, HashTable, - IndexedBits, IndexedFlexibleArray, IndexedList, IndexedOneDimensionalArray, + {AssociationList, Bits, CharacterClass, ComplexDoubleFloatVector, DataList, + DoubleFloatVector, EqTable, FlexibleArray, GeneralPolynomialSet, + GeneralSparseTable, GeneralTriangularSet, HashTable, IndexedBits, + IndexedFlexibleArray, IndexedList, IndexedOneDimensionalArray, IndexedString, IndexedVector, InnerTable, KeyedAccessFile, Library, List, - ListMultiDictionary, Multiset, OneDimensionalArray, Point, PrimitiveArray, - RegularChain, RegularTriangularSet, Result, RoutinesTable, Set, - SparseTable, SquareFreeRegularTriangularSet, Stream, String, StringTable, - Table, Vector, WuWenTsunTriangularSet} + ListMultiDictionary, Multiset, NeitherSparseOrDensePowerSeries, + OneDimensionalArray, Point, PrimitiveArray, RegularChain, + RegularTriangularSet, Result, RoutinesTable, Set, SparseTable, + SquareFreeRegularTriangularSet, Stream, String, StringTable, Table, + U32Vector, Vector, WuWenTsunTriangularSet} Type: Set Symbol This can be used to form the set-difference of two categories: @@ -27250,6 +27254,7 @@ makePoly(b) == x + b g := map(makePoly,f) --R --R Compiling function makePoly with type Integer -> Polynomial(Integer) +--R --R --R 4 2 --R (5) (x + 1)(x + 2) (x + 3) (x + 5) @@ -39543,7 +39548,7 @@ mfzn : SQMATRIX(6,DMP([x,y,z],Fraction INT)) := [ [0,1,1,1,1,1], [1,0,1,8/3,x,8/ --R | 8 8 | --R |1 - y - 1 0| --R + 3 3 + ---RType: SquareMatrix(6,DistributedMultivariatePolynomial([x,y,z],Fraction Integer)) +--RType: SquareMatrix(6,DistributedMultivariatePolynomial([x,y,z],Fraction(Integer))) --E 1 --S 2 of 3 @@ -39554,7 +39559,7 @@ eq := determinant mfzn --R 2 2 22 2 25 2 22 2 388 250 25 2 250 14575 --R - x y + -- x y - -- x + -- x y - --- x y - --- x - -- y - --- y + ----- --R 3 9 3 9 27 9 27 81 ---R Type: DistributedMultivariatePolynomial([x,y,z],Fraction Integer) +--R Type: DistributedMultivariatePolynomial([x,y,z],Fraction(Integer)) --E 2 --S 3 of 3 @@ -39585,7 +39590,7 @@ groebnerFactorize [eq,eval(eq, [x,y,z],[y,z,x]), eval(eq,[x,y,z],[z,x,y])] --R 19 5 5 --R [x - --,y + -,z + -]] --R 3 3 3 ---R Type: List List DistributedMultivariatePolynomial([x,y,z],Fraction Integer) +--RType: List(List(DistributedMultivariatePolynomial([x,y,z],Fraction(Integer)))) --E 3 )spool )lisp (bye) @@ -42267,9 +42272,12 @@ GFUN. \cross{GUESS}{shiftHP} && \end{tabular} +The original code would not compile. This is a temporary replacement +with changes marked with my initials. I will pick up the latest +version sometime in the future and hope it compiles. (Tim Daly) \begin{chunk}{package GUESS Guess} )abbrev package GUESS Guess -++ Author: Martin Rubey +++ Author: Martin Rubey, Timothy Daly ++ Description: ++ This package implements guessing of sequences. Packages for the ++ most common cases are provided as \spadtype{GuessInteger}, @@ -42370,17 +42378,9 @@ Guess(F, S, EXPRR, R, retract, coerce): Exports == Implementation where FPOLYS ==> Fraction POLYS FSUPFPOLYS ==> Fraction SUP FPOLYS -\end{chunk} -The following should be commented out when not debugging, see also -Section~\ref{sec:Hermite-Pade}. - -\begin{chunk}{package GUESS Guess} --@<> -- EXT ==> (Integer, V, V) -> FPOLYS -- EXTEXPR ==> (Symbol, F, F) -> EXPRR -\end{chunk} - -\begin{chunk}{package GUESS Guess} Exports == with guess: List F -> GUESSRESULT @@ -42573,2028 +42573,1429 @@ Section~\ref{sec:Hermite-Pade}. ++ options. --@<> -\end{chunk} - -\begin{chunk}{debug exports: Guess} ---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) _ - -> Union("failed", Record(function: EXPRR, order: NNI)) - ---checkResult: (EXPRR, Symbol, Integer, List F, LGOPT) -> NNI --- ---guessBinRatAux: (Symbol, List F, DIFFSPECN, EXT, EXTEXPR, --- List Integer, LGOPT) -> List EXPRR ---guessBinRatAux0: (List F, DIFFSPECN, EXT, EXTEXPR, --- LGOPT) -> GUESSRESULT ---binExt: EXT ---binExtEXPR: EXTEXPR ---defaultD: DIFFSPECN - -\end{chunk} - -\begin{chunk}{package GUESS Guess} Implementation == add -\getchunk{implementation: Guess} -\end{chunk} - -\begin{chunk}{implementation: Guess} -- We have to put this chunk at the beginning, because otherwise it will take -- very long to compile. -\getchunk{implementation: Guess - guessExpRat - Order and Degree} - -\getchunk{implementation: Guess - Utilities} -\getchunk{implementation: Guess - guessExpRat} -\getchunk{implementation: Guess - guessBinRat} -\getchunk{implementation: Guess - Hermite-Pade} -\getchunk{implementation: Guess - guess} -\end{chunk} - -\subsection{general utilities} - -\begin{chunk}{implementation: Guess - 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" - -\end{chunk} -%$ - -\subsection{guessing rational functions with an exponential term} - -\begin{chunk}{implementation: Guess - guessExpRat} -\getchunk{implementation: Guess - guessExpRat - Utilities} -\getchunk{implementation: Guess - guessExpRat - Main} -\end{chunk} - -\subsubsection{Utilities} - -\paragraph{conversion routines} - -\begin{chunk}{implementation: Guess - guessExpRat - Utilities} -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 + 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) - pnum/pden + 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)) + + + 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)) -POLYF2EXPRR(p: POLYF): EXPRR == - map(convert(#1)@Symbol::EXPRR, coerce(#1)@EXPRR, p) - $PolynomialCategoryLifting(IndexedExponents V, V, - F, POLYF, EXPRR) + 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) -- this needs documentation. In particular, why is V appearing here? -GF ==> GeneralizedMultivariateFactorize(SingletonAsOrderedSet, - IndexedExponents V, F, F, - SUP F) + GF ==> GeneralizedMultivariateFactorize(SingletonAsOrderedSet, + IndexedExponents V, F, F, + SUP F) -- does not work: -- 6 -- WARNING (genufact): No known algorithm to factor ? , trying square-free. -- GF ==> GenUFactorize F -\end{chunk} - -\paragraph{mimicking $q$-analoga} - -\begin{chunk}{implementation: Guess - guessExpRat - Utilities} -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) - -\getchunk{implementation: Guess - guessExpRat - evalResultant} -\getchunk{implementation: Guess - guessExpRat - substitute} -\end{chunk} - -\subsubsection{reducing the degree of the polynomials} - -The degree of poly3 is governed by $(a_0+x_m a_1)^{x_m}$. Therefore, we -substitute $A-x_m a1$ for $a_0$, which reduces the degree in $a_1$ by -$x_m-x_{i+1}$. - -\begin{chunk}{implementation: Guess - guessExpRat - substitute} -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)) - -\end{chunk} - -\begin{chunk}{not interesting implementation: Guess - guessExpRat - substitute} -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) - -\end{chunk} - -\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. - -\begin{chunk}{implementation: Guess - guessExpRat - Order and Degree} -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)) - -\end{chunk} - -evalResultant evaluates the resultant of p1 and p2 at d-o+1 -points, so that we can recover it by interpolation. - -\begin{chunk}{implementation: Guess - guessExpRat - evalResultant} -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) - -\end{chunk} - -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. - -\begin{chunk}{implementation: Guess - guessExpRat - evalResultant} - 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) - -\end{chunk} -%$ - -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. - -\begin{chunk}{implementation: Guess - guessExpRat - evalResultant} - reverse res - -\end{chunk} - -\subsubsection{The main routine} - -\begin{chunk}{implementation: Guess - guessExpRat - Main} -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)) - -\end{chunk} - -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}$. - -\begin{chunk}{implementation: Guess - guessExpRat - Main} - 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 - -\end{chunk} - -\begin{verbatim} - 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{verbatim} - -\begin{verbatim} - 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{verbatim} - -\begin{chunk}{implementation: Guess - guessExpRat - Main} - 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 --- for polynomial domains S this seems to be very costly! - 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, []) + 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 -guessExpRat(list: List F, options: LGOPT): GUESSRESULT == - guessExpRatAux0(list, defaultD, options) + 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)) -if F has RetractableTo Symbol and S has RetractableTo Symbol then + 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 + -- tpd: this is undefined since maxDegree is always nonnegative + -- if maxDegree(options)$GOPT0 = -1 + -- then maxDeg := len-1 + -- else maxDeg := min(maxDegree(options)$GOPT0, len-1) + -- maxDeg := min(maxDegree(options)$GOPT0, len-1) + tpd:Integer := (maxDegree(options)$GOPT0)::NNI::Integer + maxDeg:Integer:=min(tpd,len-1) + + for i in 0..maxDeg repeat + if debug(options)$GOPT0 then + output(hconcat("degree ExpRat "::OutputForm, i::OutputForm)) + $OutputPackage - guessExpRat(q: Symbol): GUESSER == + 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 + -- for polynomial domains S this seems to be very costly! + 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) -\end{chunk} - -\subsection{guessing rational functions with a binomial term} - -\begin{verbatim} - It is not clear whether one should take the model - \begin{equation*} - \binom{a+bn}{n}q(n), - \end{equation*} - which includes rational functions, or - \begin{equation*} - (a+bn)(a+bn+1)\dots (a+bn+n)q(n). - \end{equation*} - which includes rational functions times $n!$. We choose the former, since - dividing by $n!$ is a common normalisation. The question remains whether we - should do the same for guessExpRat. -\end{verbatim} - - -\begin{chunk}{implementation: Guess - guessBinRat} - -EXT ==> (Integer, V, V) -> FPOLYS -EXTEXPR ==> (Symbol, F, F) -> EXPRR - -binExt: EXT -binExt(i: Integer, va1: V, vA: V): FPOLYS == - numl: List POLYS := [(vA::POLYS) + i * (va1::POLYS) - (l::POLYS) _ - for l in 0..i-1] - num: POLYS := reduce(_*, numl, 1) - - num/(factorial(i)::POLYS) - -binExtEXPR: EXTEXPR -binExtEXPR(i: Symbol, a1v: F, Av: F): EXPRR == - binomial(coerce Av + coerce a1v * (i::EXPRR), i::EXPRR) - - -guessBinRatAux(xx: Symbol, list: List F, - basis: DIFFSPECN, ext: EXT, extEXPR: EXTEXPR, - 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)) - -\end{chunk} - -We try to fit the data $(s1,s2,\dots)$ to the model $\binom{a+b n}{n} y(n)$, -$r$ being a rational function. To obtain $y$, we compute -$y(n)=s_n*\binom{a+bn}{n}^-1$. - -\begin{chunk}{implementation: Guess - guessBinRat} - y: NNI -> FPOLYS := - F2FPOLYS(list.#1) / _ - ext((xValues.#1)::Integer, a1, A) - - 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 BinRat "::OutputForm, i::OutputForm)) --- $OutputPackage - -\end{chunk} - -\begin{verbatim} - 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{verbatim} - -\begin{verbatim} - 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{verbatim} - -\begin{chunk}{implementation: Guess - guessBinRat} --- if debug(options)$GOPT0 then --- output("interpolating..."::OutputForm)$OutputPackage - - ri: FSUPFPOLYS - := interpolate(xlist, ylist, (len-1-i)::NNI) _ - $FFFG(FPOLYS, SUP FPOLYS) - --- if debug(options)$GOPT0 then --- output(hconcat("ri ", ri::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) - --- if debug(options)$GOPT0 then --- output(hconcat("poly1 ", poly1::OutputForm))$OutputPackage --- output(hconcat("poly2 ", poly2::OutputForm))$OutputPackage --- output(hconcat("poly3 ", poly3::OutputForm))$OutputPackage - - - n:Integer := len - i - res1: SUP S := univariate(resultant(poly1, poly3, a1)) - res2: SUP S := univariate(resultant(poly2, poly3, a1)) - if debug(options)$GOPT0 then --- output(hconcat("res1 ", res1::OutputForm))$OutputPackage --- output(hconcat("res2 ", res2::OutputForm))$OutputPackage - --- 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(minimumDegree res1)) - $OutputPackage - output("res1 deg: " string(degree res1)) - $OutputPackage - output("res2 ord: " string(minimumDegree res2)) - $OutputPackage - output("res2 deg: " string(degree res2)) - $OutputPackage - - if debug(options)$GOPT0 then - output("computing gcd..."::OutputForm)$OutputPackage - --- we want to solve over F - res3: SUP F := SUPS2SUPF(primitivePart(gcd(res1, res2))) - --- if debug(options)$GOPT0 then --- output(hconcat("res3 ", res3::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 - --- if debug(options)$GOPT0 then --- output(hconcat("Av: ", Av::OutputForm))$OutputPackage - --- 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 - --- if debug(options)$GOPT0 then --- output(hconcat("a1v: ", a1v::OutputForm))$OutputPackage - --- 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 debug(options)$GOPT0 then --- output(hconcat("t1: ", t1::OutputForm))$OutputPackage - - if zero? t1 then - t2 := eval(POLYS2POLYF poly2, [a1, A]::List V, - [a1v, Av]::List F) - --- if debug(options)$GOPT0 then --- output(hconcat("t2: ", t2::OutputForm))$OutputPackage - - if zero? t2 then - - ri1: Fraction SUP POLYS - := SUPFPOLYS2FSUPPOLYS(numer ri) - / SUPFPOLYS2FSUPPOLYS(denom ri) - --- if debug(options)$GOPT0 then --- output(hconcat("ri1: ", ri1::OutputForm))$OutputPackage - - numr: SUP F := SUPPOLYS2SUPF(numer ri1, a1v, Av) - denr: SUP F := SUPPOLYS2SUPF(denom ri1, a1v, Av) - --- if debug(options)$GOPT0 then --- output(hconcat("numr: ", numr::OutputForm))$OutputPackage --- output(hconcat("denr: ", denr::OutputForm))$OutputPackage - - if not zero? denr then - res4: EXPRR := eval(FSUPF2EXPRR(xx, numr/denr), - kernel(xx), - basis(xx::EXPRR)) - * extEXPR(xx, a1v, Av) - --- if debug(options)$GOPT0 then --- output(hconcat("res4: ", res4::OutputForm))$OutputPackage - - 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 - -guessBinRatAux0(list: List F, - basis: DIFFSPECN, ext: EXT, extEXPR: EXTEXPR, - options: LGOPT): GUESSRESULT == - - if zero? safety(options)$GOPT0 then - error "Guess: guessBinRat does not support zero safety" --- guesses Functions of the form binomial(a+b*n, 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 guessBinRatAux(xx, newlist, basis, ext, extEXPR, 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) - -guessBinRat(list : List F): GUESSRESULT == - guessBinRatAux0(list, defaultD, binExt, binExtEXPR, []) - -guessBinRat(list: List F, options: LGOPT): GUESSRESULT == - guessBinRatAux0(list, defaultD, binExt, binExtEXPR, options) - - -if F has RetractableTo Symbol and S has RetractableTo Symbol then - - qD: Symbol -> DIFFSPECN - qD q == (q::EXPRR)**#1 - - - qBinExtAux(q: Symbol, i: Integer, va1: V, vA: V): FPOLYS == - fl: List FPOLYS - := [(1$FPOLYS - _ - va1::POLYS::FPOLYS * (vA::POLYS::FPOLYS)**(i-1) * _ - F2FPOLYS(q::F)**l) / (1-F2FPOLYS(q::F)**l) _ - for l in 1..i] - reduce(_*, fl, 1) - - qBinExt: Symbol -> EXT - qBinExt q == qBinExtAux(q, #1, #2, #3) - - qBinExtEXPRaux(q: Symbol, i: Symbol, a1v: F, Av: F): EXPRR == - l: Symbol := 'l - product((1$EXPRR - _ - coerce a1v * (coerce Av) ** (coerce i - 1$EXPRR) * _ - (q::EXPRR) ** coerce(l)) / _ - (1$EXPRR - (q::EXPRR) ** coerce(l)), _ - equation(l, 1$EXPRR..i::EXPRR)) - - qBinExtEXPR: Symbol -> EXTEXPR - qBinExtEXPR q == qBinExtEXPRaux(q, #1, #2, #3) - - guessBinRat(q: Symbol): GUESSER == - guessBinRatAux0(#1, qD q, qBinExt q, qBinExtEXPR q, #2)_ - -\end{chunk} - - -\subsection{Hermite Pad\'e interpolation}\label{sec:Hermite-Pade} - -\begin{chunk}{implementation: Guess - Hermite-Pade} -\getchunk{implementation: Guess - Hermite-Pade - Types for Operators} -\getchunk{implementation: Guess - Hermite-Pade - Streams} -\getchunk{implementation: Guess - Hermite-Pade - Operators} -\getchunk{implementation: Guess - Hermite-Pade - Utilities} -\getchunk{implementation: Guess - guessHPaux} -\getchunk{implementation: Guess - guessHP} -\end{chunk} - -\subsubsection{Types for Operators} -\begin{chunk}{implementation: Guess - Hermite-Pade - 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 - -\end{chunk} - -\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. - -\begin{chunk}{implementation: Guess - Hermite-Pade - Streams} -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) - -\end{chunk} -%$ - -It is not clear whether we should \lq prefer\rq\ shifting and differentiation over -powering. Currently, we produce the stream -\[ - \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} -\] - -Maybe it would be better to produce -\[ - \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} -\] -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 - -\begin{chunk}{implementation: Guess - Hermite-Pade - Streams} -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) - -\end{chunk} -%$ - -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. - -\begin{chunk}{implementation: Guess - Hermite-Pade - Streams} -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) - -\end{chunk} -\subsubsection{Operators} - -We need to define operators that transform series for differentiation and -shifting. We also provide operators for $q$-analogs. 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)$. - -\begin{chunk}{implementation: Guess - Hermite-Pade - Operators} -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) - -\end{chunk} - -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)$. - -\begin{chunk}{implementation: Guess - Hermite-Pade - Operators} -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" - -\end{chunk} - -\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)$. - -\begin{chunk}{implementation: Guess - Hermite-Pade - Operators} -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" - -\end{chunk} - -\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? -\begin{chunk}{implementation: Guess - Hermite-Pade - Operators} -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) - -\end{chunk} -%$ - -As before, next three functions extract the coefficient of $x^k$ in $z^l f(x)$. - -\begin{chunk}{implementation: Guess - Hermite-Pade - Operators} -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 - -\end{chunk} - -\paragraph{$q$-Shifting} The $q$-shift also transforms $u(n)$ into $u(n+1)$, -and we can reuse the corresponding 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? -\begin{chunk}{implementation: Guess - Hermite-Pade - Operators} -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 - -\end{chunk} -%$ -\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. - -\begin{chunk}{implementation: Guess - Hermite-Pade - Operators} -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) - -\end{chunk} -%$ - -\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. - -\begin{chunk}{implementation: Guess - Hermite-Pade - Utilities} -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 - -\end{chunk} - -SUPF2SUPSUPF interprets each coefficient as a univariate polynomial. - -\begin{chunk}{implementation: Guess - Hermite-Pade - Utilities} -SUPF2SUPSUPF(p: SUP F): SUP SUP F == - map(#1::SUP F, p)$SparseUnivariatePolynomialFunctions2(F, SUP F) - -\end{chunk} - -getListSUPF returns the first o elements of the stream, each truncated -after degree deg. - -\begin{chunk}{implementation: Guess - Hermite-Pade - Utilities} -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" - -\getchunk{implementation: Guess - guessInterpolate} -\getchunk{implementation: Guess - guessInterpolate2} -\getchunk{implementation: Guess - testInterpolant} -\end{chunk} -%$ - -guessInterpolate calls the appropriate generalInterpolation from -FFFG, for one vector of degrees, namely eta. - -\begin{chunk}{implementation: Guess - guessInterpolate} -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" -\end{chunk} - -guessInterpolate2 calls the appropriate generalInterpolation from -FFFG, for all degree vectors with given sumEta and maxEta. - -\begin{chunk}{implementation: Guess - guessInterpolate2} -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" -\end{chunk} - -testInterpolant checks whether p is really a solution. -\begin{chunk}{implementation: Guess - testInterpolant} -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) - : Union("failed", Record(function: EXPRR, order: NNI)) == -\end{chunk} - -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. - -\begin{chunk}{implementation: Guess - testInterpolant} - ((maxDegree(options)$GOPT0 = -1) and - (allDegrees(options)$GOPT0 = false) and - zero?(last resi)) - => return "failed" -\end{chunk} - -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. - -\begin{chunk}{implementation: Guess - testInterpolant} - 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 -\end{chunk} - -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{verbatim} - Check that not all coefficients in resi can vanish simultaneously. -\end{verbatim} - -\begin{chunk}{implementation: Guess - testInterpolant} - if not zero? nonZeroCoefficient then - (freeOf?(exprList.nonZeroCoefficient, name op)) => return "failed" - - for e in list repeat - if not zero? e then return "failed" -\end{chunk} - -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. - -\begin{chunk}{implementation: Guess - testInterpolant} - 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" -\end{chunk} - -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. - -If S is a Field, the gcd will always be $1$. Thus, in this case we -make the result monic. - -\begin{chunk}{implementation: Guess - testInterpolant} - g: SUP S - if S has Field - then g := leadingCoefficient(find(not zero? #1, reverse resi)::SUP(S))::SUP(S) - else g := 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) - -\end{chunk} - -\subsubsection{The main routine} - -The following is the main guessing routine, called by all others -- except -guessExpRat. - -\begin{chunk}{implementation: Guess - guessHPaux} -guessHPaux(list: List F, D: HPSPEC, options: LGOPT): GUESSRESULT == - reslist: GUESSRESULT := [] - - listDegree := #list-1-safety(options)$GOPT0 - if listDegree < 0 then return reslist -\end{chunk} -%$ - -\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. - -\begin{chunk}{implementation: Guess - guessHPaux} - a := functionName(options)$GOPT0 - op := operator a - x := variableName(options)$GOPT0 - dummy := new$Symbol - - initials: List EXPRR := [coerce(e)@EXPRR for e in list] - -\end{chunk} -%$ - -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{verbatim} - 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{verbatim} - - -\begin{chunk}{implementation: Guess - guessHPaux} - guessS := (D.guessStream)(list2UFPSF list) - degreeS := D.degreeStream - testS := (D.testStream)(list2UFPSSUPF list) - exprS := (D.exprStream)(op(dummy::EXPRR)::EXPRR, dummy) -\end{chunk} - -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{verbatim} - f &= l_0 + l_1 x + l_2 x^2\\ - f^\prime &= l_1 + 2l_2 x\\ - f^{\prime\prime} &= 2l_2. -\end{verbatim} -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. - -\begin{chunk}{implementation: Guess - guessHPaux} - 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 -\end{chunk} -%$" - -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. - -\begin{chunk}{implementation: Guess - guessHPaux} - if allDegrees(options)$GOPT0 then - (o > guessDegree+2) => return reslist - - maxEta: Integer := 1+maxDegree(options)$GOPT0 - if maxEta = 0 - then maxEta := guessDegree+1 - else -\end{chunk} + EXT ==> (Integer, V, V) -> FPOLYS + EXTEXPR ==> (Symbol, F, F) -> EXPRR + + binExt: EXT + binExt(i: Integer, va1: V, vA: V): FPOLYS == + numl: List POLYS := [(vA::POLYS) + i * (va1::POLYS) - (l::POLYS) _ + for l in 0..i-1] + num: POLYS := reduce(_*, numl, 1) + + num/(factorial(i)::POLYS) + + binExtEXPR: EXTEXPR + binExtEXPR(i: Symbol, a1v: F, Av: F): EXPRR == + binomial(coerce Av + coerce a1v * (i::EXPRR), i::EXPRR) + + + guessBinRatAux(xx: Symbol, list: List F, + basis: DIFFSPECN, ext: EXT, extEXPR: EXTEXPR, + 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) / _ + ext((xValues.#1)::Integer, a1, A) + + 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 + -- tpd: this is undefined since maxDegree is always nonnegative + -- if maxDegree(options)$GOPT0 = -1 + -- then maxDeg := len-1 + -- else maxDeg := min(maxDegree(options)$GOPT0, len-1) + -- maxDeg := min(maxDegree(options)$GOPT0, len-1) + tpd:Integer := (maxDegree(options)$GOPT0)::NNI::Integer + maxDeg:Integer:=min(tpd,len-1) +-- maxDeg:Integer := (maxDegree(options)$GOPT0)::NNI::Integer + + for i in 0..maxDeg repeat + -- if debug(options)$GOPT0 then + -- output(hconcat("degree BinRat "::OutputForm, i::OutputForm)) + -- $OutputPackage -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. + -- if debug(options)$GOPT0 then + -- output("interpolating..."::OutputForm)$OutputPackage + + ri: FSUPFPOLYS + := interpolate(xlist, ylist, (len-1-i)::NNI) _ + $FFFG(FPOLYS, SUP FPOLYS) + + -- if debug(options)$GOPT0 then + -- output(hconcat("ri ", ri::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) + + -- if debug(options)$GOPT0 then + -- output(hconcat("poly1 ", poly1::OutputForm))$OutputPackage + -- output(hconcat("poly2 ", poly2::OutputForm))$OutputPackage + -- output(hconcat("poly3 ", poly3::OutputForm))$OutputPackage + + + n:Integer := len - i + res1: SUP S := univariate(resultant(poly1, poly3, a1)) + res2: SUP S := univariate(resultant(poly2, poly3, a1)) + if debug(options)$GOPT0 then + -- output(hconcat("res1 ", res1::OutputForm))$OutputPackage + -- output(hconcat("res2 ", res2::OutputForm))$OutputPackage + + -- 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(minimumDegree res1)) + $OutputPackage + output("res1 deg: " string(degree res1)) + $OutputPackage + output("res2 ord: " string(minimumDegree res2)) + $OutputPackage + output("res2 deg: " string(degree res2)) + $OutputPackage + + if debug(options)$GOPT0 then + output("computing gcd..."::OutputForm)$OutputPackage + + -- we want to solve over F + res3: SUP F := SUPS2SUPF(primitivePart(gcd(res1, res2))) + + -- if debug(options)$GOPT0 then + -- output(hconcat("res3 ", res3::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 + + -- if debug(options)$GOPT0 then + -- output(hconcat("Av: ", Av::OutputForm))$OutputPackage + + -- 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 + + -- if debug(options)$GOPT0 then + -- output(hconcat("a1v: ", a1v::OutputForm))$OutputPackage + + -- 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 debug(options)$GOPT0 then + -- output(hconcat("t1: ", t1::OutputForm))$OutputPackage + + if zero? t1 then + t2 := eval(POLYS2POLYF poly2, [a1, A]::List V, + [a1v, Av]::List F) + + -- if debug(options)$GOPT0 then + -- output(hconcat("t2: ", t2::OutputForm))$OutputPackage + + if zero? t2 then + + ri1: Fraction SUP POLYS + := SUPFPOLYS2FSUPPOLYS(numer ri) + / SUPFPOLYS2FSUPPOLYS(denom ri) + + -- if debug(options)$GOPT0 then + -- output(hconcat("ri1: ", ri1::OutputForm))$OutputPackage + + numr: SUP F := SUPPOLYS2SUPF(numer ri1, a1v, Av) + denr: SUP F := SUPPOLYS2SUPF(denom ri1, a1v, Av) + + -- if debug(options)$GOPT0 then + -- output(hconcat("numr: ", numr::OutputForm))$OutputPackage + -- output(hconcat("denr: ", denr::OutputForm))$OutputPackage + + if not zero? denr then + res4: EXPRR := eval(FSUPF2EXPRR(xx, numr/denr), + kernel(xx), + basis(xx::EXPRR)) + * extEXPR(xx, a1v, Av) + + -- if debug(options)$GOPT0 then + -- output(hconcat("res4: ", res4::OutputForm))$OutputPackage + + 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 + + guessBinRatAux0(list: List F, + basis: DIFFSPECN, ext: EXT, extEXPR: EXTEXPR, + options: LGOPT): GUESSRESULT == + + if zero? safety(options)$GOPT0 then + error "Guess: guessBinRat does not support zero safety" + -- guesses Functions of the form binomial(a+b*n, 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 guessBinRatAux(xx, newlist, basis, ext, extEXPR, 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) + + guessBinRat(list : List F): GUESSRESULT == + guessBinRatAux0(list, defaultD, binExt, binExtEXPR, []) + + guessBinRat(list: List F, options: LGOPT): GUESSRESULT == + guessBinRatAux0(list, defaultD, binExt, binExtEXPR, options) + + + if F has RetractableTo Symbol and S has RetractableTo Symbol then + + qD: Symbol -> DIFFSPECN + qD q == (q::EXPRR)**#1 + + + qBinExtAux(q: Symbol, i: Integer, va1: V, vA: V): FPOLYS == + fl: List FPOLYS + := [(1$FPOLYS - _ + va1::POLYS::FPOLYS * (vA::POLYS::FPOLYS)**(i-1) * _ + F2FPOLYS(q::F)**l) / (1-F2FPOLYS(q::F)**l) _ + for l in 1..i] + reduce(_*, fl, 1) + + qBinExt: Symbol -> EXT + qBinExt q == qBinExtAux(q, #1, #2, #3) + + qBinExtEXPRaux(q: Symbol, i: Symbol, a1v: F, Av: F): EXPRR == + l: Symbol := 'l + product((1$EXPRR - _ + coerce a1v * (coerce Av) ** (coerce i - 1$EXPRR) * _ + (q::EXPRR) ** coerce(l)) / _ + (1$EXPRR - (q::EXPRR) ** coerce(l)), _ + equation(l, 1$EXPRR..i::EXPRR)) + + qBinExtEXPR: Symbol -> EXTEXPR + qBinExtEXPR q == qBinExtEXPRaux(q, #1, #2, #3) + + guessBinRat(q: Symbol): GUESSER == + guessBinRatAux0(#1, qD q, qBinExt q, qBinExtEXPR q, #2)_ -\begin{chunk}{implementation: Guess - guessHPaux} - maxParams := divide(guessDegree::NNI+1, o) - if debug(options)$GOPT0 - then output(hconcat("maxParams: ", maxParams::OutputForm)) - $OutputPackage -\end{chunk} + -- 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 -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. + 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) -\begin{chunk}{implementation: Guess - guessHPaux} - if maxParams.quotient = 0 and maxParams.remainder < o-1 - then return reslist -\end{chunk} + FilteredPartitionStream(options: LGOPT): Stream List Integer == + --tpd: must force types to NNI and Integer + maxD:Integer := 1+maxDerivative(options)$GOPT0::NNI::Integer + maxP:Integer := maxPower(options)$GOPT0::NNI::Integer + + 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 + --tpd: force the Boolean branch + tpd2:Boolean:=homogeneous(options)$GOPT0::Boolean + if tpd2 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) -If maxDegree is set, we skip the first few partitions, unless we cannot -increase the order anymore. -\begin{verbatim} - I have no satisfactory way to determine whether we can increase the order or - not. -\end{verbatim} + 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" -\begin{chunk}{implementation: Guess - guessHPaux} - 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 if ((maxDegree(options)$GOPT0 ~= -1) and - (maxParams.quotient > maxDegree(options)$GOPT0)) - then - guessDegree := o*(1+maxDegree(options)$GOPT0)-2 - eta: List NNI - := [(if i < o _ - then maxDegree(options)$GOPT0 + 1 _ - else maxDegree(options)$GOPT0)::NNI _ - for i in 1..o] - else eta: List NNI - := [(if i <= maxParams.remainder _ - then maxParams.quotient + 1 _ - else maxParams.quotient)::NNI for i in 1..o] -\end{chunk} - -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. + 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" -\begin{verbatim} - maxDegree should be handled differently, maybe: we should only pass as - many coefficients to FFFG as maxDegree implies! That is, if we cannot - increase the order anymore, we should decrease guessDegree to - $o\cdot maxDegree - 2$ and set eta accordingly. I might have to take - care of allDegrees, too. -\end{verbatim} + 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) -\begin{chunk}{implementation: Guess - guessHPaux} - 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) - - (res case "failed") => "iterate" - - if not member?(res, reslist) - then reslist := cons(res, reslist) - - if one(options)$GOPT0 then return reslist - - MS := rest MS + 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 - 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) - (res case "failed") => "iterate" - - if not member?(res, reslist) - then reslist := cons(res, reslist) - - if one(options)$GOPT0 then return reslist - - reslist - -\end{chunk} - -\begin{verbatim} - The duplicated block at the end should really go into testInterpolant, I - guess. Furthermore, it would be better to remove duplicates already there. -\end{verbatim} - -\subsubsection{Specialisations} - -For convenience we provide some specialisations that cover the most common -situations. - -\paragraph{generating functions} - -\begin{chunk}{implementation: Guess - guessHP} -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, []) -\end{chunk} - -\paragraph{$q$-generating functions} - -\begin{chunk}{implementation: Guess - guessHP} -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) - -\end{chunk} -%$ - -\paragraph{coefficients} - -\begin{chunk}{implementation: Guess - guessHP} -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, []) - -\end{chunk} -%$ - -\paragraph{$q$-coefficients} - -\begin{chunk}{implementation: Guess - guessHP} -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) + 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 -\end{chunk} -%$ + 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 + + 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) -\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} - -\begin{chunk}{implementation: Guess - guess} -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 + 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) + : Union("failed", Record(function: EXPRR, order: NNI)) == + --tpd: maxDegree is defined to be nonnegative + -- ((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 + if S has Field + then g := leadingCoefficient(find(not zero? #1, reverse resi)::SUP(S))::SUP(S) + else g := 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) - 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) + 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 := [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 + --tpd: force NNI and Integer + maxEta: Integer := 1+maxDegree(options)$GOPT0::NNI::Integer + 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 + --tpd: maxDegree is defined to be nonnegative + -- if ((maxDegree(options)$GOPT0 ~= -1) and + if ((maxDegree(options)$GOPT0::NNI::Integer < 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 + --tpd:maxDegree is defined to be nonnegative + -- else if ((maxDegree(options)$GOPT0 ~= -1) and + if (maxParams.quotient > maxDegree(options)$GOPT0::NNI::Integer) + then + --tpd:maxDegree is defined to be nonnegative + guessDegree := o*(1+maxDegree(options)$GOPT0::NNI::Integer)-2 + eta: List NNI + := [(if i < o _ + then maxDegree(options)$GOPT0::NNI + 1 _ + else maxDegree(options)$GOPT0::NNI) _ + for i in 1..o] + 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) + + (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) + (res case "failed") => "iterate" + + if not member?(res, reslist) + then reslist := cons(res, reslist) + + if one(options)$GOPT0 then return reslist + + reslist - res + 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, []) -guess(list: List F): GUESSRESULT == - guess(list, [guessADE$%, guessRec$%], ['guessProduct, 'guessSum], []) + 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) -guess(list: List F, options: LGOPT): GUESSRESULT == - guess(list, [guessADE$%, guessRat$%], ['guessProduct, 'guessSum], - options) + 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, []) -guess(list: List F, guessers: List GUESSER, ops: List Symbol) - : GUESSRESULT == - guess(list, guessers, ops, []) + 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)] + + -- tpd: maxLevel is NNI + 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, options)$%] +-- tpd: this is broken +-- append([(indexName(var)$GuessOption)::Symbol,_ +-- (maxLevel(maxLevel-1)$GuessOption)::NNI],_ +-- 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)] + -- tpd:maxLevel is NNI + 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, options)$%] +--tpd: this is broken +-- for guess in guess(sumList, guessers, ops,_ +-- append([(indexName(var)$GuessOption)::Symbol,_ +-- (maxLevel(maxLevel-1)$GuessOption)::NNI],_ +-- 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{chunk} \begin{chunk}{GUESS.dotabb} diff --git a/books/bookvol10.5.pamphlet b/books/bookvol10.5.pamphlet index d5e32c5..f308ad0 100644 --- a/books/bookvol10.5.pamphlet +++ b/books/bookvol10.5.pamphlet @@ -353,7 +353,7 @@ t1:Complex DoubleFloat := complex(1.0,0) --R --R --R (1) 1. ---R Type: Complex DoubleFloat +--R Type: Complex(DoubleFloat) --E 1 --S 2 of 10 @@ -369,7 +369,7 @@ t2:Complex DoubleFloat := complex(1.0,1.0) --R --R --R (3) 1. + %i ---R Type: Complex DoubleFloat +--R Type: Complex(DoubleFloat) --E 3 --S 4 of 10 @@ -385,7 +385,7 @@ t3:Complex DoubleFloat := complex(1.0,-1.0) --R --R --R (5) 1. - %i ---R Type: Complex DoubleFloat +--R Type: Complex(DoubleFloat) --E 5 --S 6 of 10 @@ -401,7 +401,7 @@ t4:Complex DoubleFloat := complex(-1.0,-1.0) --R --R --R (7) - 1. - %i ---R Type: Complex DoubleFloat +--R Type: Complex(DoubleFloat) --E 7 --S 8 of 10 @@ -417,7 +417,7 @@ t5:Complex DoubleFloat := complex(-2.0,-2.0) --R --R --R (9) - 2. - 2. %i ---R Type: Complex DoubleFloat +--R Type: Complex(DoubleFloat) --E 9 --S 10 of 10 @@ -601,7 +601,7 @@ t1:Complex DoubleFloat := complex(1.0,0) --R --R --R (1) 1. ---R Type: Complex DoubleFloat +--R Type: Complex(DoubleFloat) --E 1 --S 2 of 10 @@ -617,7 +617,7 @@ t2:Complex DoubleFloat := complex(1.0,1.0) --R --R --R (3) 1. + %i ---R Type: Complex DoubleFloat +--R Type: Complex(DoubleFloat) --E 3 --S 4 of 10 @@ -633,7 +633,7 @@ t3:Complex DoubleFloat := complex(1.0,-1.0) --R --R --R (5) 1. - %i ---R Type: Complex DoubleFloat +--R Type: Complex(DoubleFloat) --E 5 --S 6 of 10 @@ -649,7 +649,7 @@ t4:Complex DoubleFloat := complex(-1.0,-1.0) --R --R --R (7) - 1. - %i ---R Type: Complex DoubleFloat +--R Type: Complex(DoubleFloat) --E 7 --S 8 of 10 @@ -665,7 +665,7 @@ t5:Complex DoubleFloat := complex(-2.0,-2.0) --R --R --R (9) - 2. - 2. %i ---R Type: Complex DoubleFloat +--R Type: Complex(DoubleFloat) --E 9 --S 10 of 10 @@ -797,7 +797,7 @@ a:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0] ] --R --R --R (1) [1.,2.,3.,4.,5.,6.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 1 --S 2 of 28 @@ -1242,7 +1242,7 @@ a:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0] ] --R --R --R (1) [1.,2.,3.,4.,5.,6.,7.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 1 --S 2 of 22 @@ -1250,7 +1250,7 @@ b:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0] ] --R --R --R (2) [1.,2.,3.,4.,5.,6.,7.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 2 --S 3 of 22 @@ -1258,7 +1258,7 @@ daxpy(3,2.0,a,1,b,1) --R --R --R (3) [3.,6.,9.,4.,5.,6.,7.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 3 --S 4 of 22 @@ -1266,7 +1266,7 @@ b:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0] ] --R --R --R (4) [1.,2.,3.,4.,5.,6.,7.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 4 --S 5 of 22 @@ -1274,7 +1274,7 @@ daxpy(7,2.0,a,1,b,1) --R --R --R (5) [3.,6.,9.,12.,15.,18.,21.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 5 --S 6 of 22 @@ -1282,7 +1282,7 @@ b:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0] ] --R --R --R (6) [1.,2.,3.,4.,5.,6.,7.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 6 --S 7 of 22 @@ -1290,7 +1290,7 @@ daxpy(8,2.0,a,1,b,1) --R --R --R (7) [1.,2.,3.,4.,5.,6.,7.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 7 --S 8 of 22 @@ -1298,7 +1298,7 @@ b:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0] ] --R --R --R (8) [1.,2.,3.,4.,5.,6.,7.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 8 --S 9 of 22 @@ -1306,7 +1306,7 @@ daxpy(3,2.0,a,3,b,3) --R --R --R (9) [3.,2.,3.,12.,5.,6.,21.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 9 --S 10 of 22 @@ -1314,7 +1314,7 @@ b:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0] ] --R --R --R (10) [1.,2.,3.,4.,5.,6.,7.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 10 --S 11 of 22 @@ -1322,7 +1322,7 @@ daxpy(4,2.0,a,2,b,2) --R --R --R (11) [3.,2.,9.,4.,15.,6.,21.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 11 --S 12 of 22 @@ -1330,7 +1330,7 @@ b:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0] ] --R --R --R (12) [1.,2.,3.,4.,5.,6.,7.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 12 --S 13 of 22 @@ -1338,7 +1338,7 @@ daxpy(5,2.0,a,2,b,2) --R --R --R (13) [1.,2.,3.,4.,5.,6.,7.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 13 --S 14 of 22 @@ -1346,7 +1346,7 @@ b:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0] ] --R --R --R (14) [1.,2.,3.,4.,5.,6.,7.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 14 --S 15 of 22 @@ -1354,7 +1354,7 @@ daxpy(3,2.0,a,2,b,2) --R --R --R (15) [3.,2.,9.,4.,15.,6.,7.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 15 --S 16 of 22 @@ -1362,7 +1362,7 @@ b:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0] ] --R --R --R (16) [1.,2.,3.,4.,5.,6.,7.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 16 --S 17 of 22 @@ -1370,7 +1370,7 @@ daxpy(3,-2.0,a,2,b,2) --R --R --R (17) [- 1.,2.,- 3.,4.,- 5.,6.,7.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 17 --S 18 of 22 @@ -1378,7 +1378,7 @@ a:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0] ] --R --R --R (18) [1.,2.,3.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 18 --S 19 of 22 @@ -1386,7 +1386,7 @@ b:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0] ] --R --R --R (19) [1.,2.,3.,4.,5.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 19 --S 20 of 22 @@ -1394,7 +1394,7 @@ daxpy(3,-2.0,a,1,b,2) --R --R --R (20) [- 1.,2.,- 1.,4.,- 1.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 20 --S 21 of 22 @@ -1402,7 +1402,7 @@ b:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0] ] --R --R --R (21) [1.,2.,3.,4.,5.,6.,7.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 21 --S 22 of 22 @@ -1410,7 +1410,7 @@ daxpy(3,0.0,a,1,b,2) --R --R --R (22) [1.,2.,3.,4.,5.,6.,7.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 22 )spool @@ -1628,7 +1628,7 @@ a:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0] ] --R --R --R (1) [1.,2.,3.,4.,5.,6.,7.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 1 --S 2 of 23 @@ -1636,7 +1636,7 @@ b:PRIMARR(DFLOAT):=[ [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] ] --R --R --R (2) [0.,0.,0.,0.,0.,0.,0.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 2 --S 3 of 23 @@ -1644,7 +1644,7 @@ dcopy(3,a,1,b,1) --R --R --R (3) [1.,2.,3.,0.,0.,0.,0.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 3 --S 4 of 23 @@ -1652,7 +1652,7 @@ b:PRIMARR(DFLOAT):=[ [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] ] --R --R --R (4) [0.,0.,0.,0.,0.,0.,0.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 4 --S 5 of 23 @@ -1660,7 +1660,7 @@ dcopy(7,a,1,b,1) --R --R --R (5) [1.,2.,3.,4.,5.,6.,7.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 5 --S 6 of 23 @@ -1668,7 +1668,7 @@ b:PRIMARR(DFLOAT):=[ [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] ] --R --R --R (6) [0.,0.,0.,0.,0.,0.,0.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 6 --S 7 of 23 @@ -1676,7 +1676,7 @@ dcopy(8,a,1,b,1) --R --R --R (7) [0.,0.,0.,0.,0.,0.,0.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 7 --S 8 of 23 @@ -1684,7 +1684,7 @@ b:PRIMARR(DFLOAT):=[ [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] ] --R --R --R (8) [0.,0.,0.,0.,0.,0.,0.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 8 --S 9 of 23 @@ -1692,7 +1692,7 @@ dcopy(3,a,3,b,3) --R --R --R (9) [1.,0.,0.,4.,0.,0.,7.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 9 --S 10 of 23 @@ -1700,7 +1700,7 @@ b:PRIMARR(DFLOAT):=[ [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] ] --R --R --R (10) [0.,0.,0.,0.,0.,0.,0.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 10 --S 11 of 23 @@ -1708,7 +1708,7 @@ dcopy(4,a,2,b,2) --R --R --R (11) [1.,0.,3.,0.,5.,0.,7.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 11 --S 12 of 23 @@ -1716,7 +1716,7 @@ b:PRIMARR(DFLOAT):=[ [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] ] --R --R --R (12) [0.,0.,0.,0.,0.,0.,0.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 12 --S 13 of 23 @@ -1724,7 +1724,7 @@ dcopy(5,a,2,b,2) --R --R --R (13) [0.,0.,0.,0.,0.,0.,0.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 13 --S 14 of 23 @@ -1732,7 +1732,7 @@ b:PRIMARR(DFLOAT):=[ [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] ] --R --R --R (14) [0.,0.,0.,0.,0.,0.,0.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 14 --S 15 of 23 @@ -1740,7 +1740,7 @@ dcopy(3,a,2,b,2) --R --R --R (15) [1.,0.,3.,0.,5.,0.,0.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 15 --S 16 of 23 @@ -1748,7 +1748,7 @@ a:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0] ] --R --R --R (16) [1.,2.,3.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 16 --S 17 of 23 @@ -1756,7 +1756,7 @@ b:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0] ] --R --R --R (17) [1.,2.,3.,4.,5.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 17 --S 18 of 23 @@ -1764,7 +1764,7 @@ dcopy(3,a,1,b,1) --R --R --R (18) [1.,2.,3.,4.,5.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 18 --S 19 of 23 @@ -1772,7 +1772,7 @@ b:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0] ] --R --R --R (19) [1.,2.,3.,4.,5.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 19 --S 20 of 23 @@ -1780,7 +1780,7 @@ dcopy(3,a,1,b,2) --R --R --R (20) [1.,2.,2.,4.,3.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 20 --S 21 of 23 @@ -1788,7 +1788,7 @@ a:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0, 4.0, 5.0] ] --R --R --R (21) [1.,2.,3.,4.,5.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 21 --S 22 of 23 @@ -1796,7 +1796,7 @@ b:PRIMARR(DFLOAT):=[ [ 1.0, 2.0, 3.0] ] --R --R --R (22) [1.,2.,3.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 22 --S 23 of 23 @@ -1804,7 +1804,7 @@ dcopy(5,a,1,b,1) --R --R --R (23) [1.,2.,3.] ---R Type: PrimitiveArray DoubleFloat +--R Type: PrimitiveArray(DoubleFloat) --E 23 )spool diff --git a/changelog b/changelog index 2c8bb2e..858c1a2 100644 --- a/changelog +++ b/changelog @@ -1,3 +1,10 @@ +20120325 tpd src/axiom-website/patches.html 20120325.01.tpd.patch +20120325 tpd src/input/hyperell.input fix 7217 +20120325 tpd src/input/exampleagcode.input fix bug 7217 +20120325 tpd src/algebra/Makefile fix bug 7217 +20120325 tpd books/bookvol10.5 fix bug 7217 +20120325 tpd books/bookvol10.4 fix bug 7217 +20120325 tpd books/bookvol10.3 fix bug 7217 20120324 tpd src/axiom-website/patches.html 20120324.02.tpd.patch 20120324 tpd books/bookvol10.4 add getAncestors function to API domain 20120324 tpd src/axiom-website/patches.html 20120324.01.tpd.patch diff --git a/src/algebra/Makefile.pamphlet b/src/algebra/Makefile.pamphlet index d2177ad..caff50f 100644 --- a/src/algebra/Makefile.pamphlet +++ b/src/algebra/Makefile.pamphlet @@ -16939,23 +16939,23 @@ LAYER23=\ \subsection{Order} The final order of the layers is determined here. The GUESS package is broken so we remove the layers involved until this can be resolved. -\begin{verbatim} +<>= ORDER=\ ${LAYER0BOOTSTRAP} ${LAYER0} ${LAYER1} ${LAYER2} ${LAYER3} \ ${LAYER4} ${LAYER5} ${LAYER6} ${LAYER7} ${LAYER8} ${LAYER9} \ ${LAYER10} ${LAYER11} ${LAYER12} ${LAYER13} ${LAYER14} ${LAYER15} \ ${LAYER16} ${LAYER17} ${LAYER18} ${LAYER19} ${LAYER20} ${LAYER21} \ ${LAYER22} ${LAYER23} ${LAYER0COPY} -\end{verbatim} -<>= +@ +\begin{verbatim} ORDER=\ ${LAYER0BOOTSTRAP} ${LAYER0} ${LAYER1} ${LAYER2} ${LAYER3} \ ${LAYER4} ${LAYER5} ${LAYER6} ${LAYER7} ${LAYER8} ${LAYER9} \ ${LAYER10} ${LAYER11} ${LAYER12} ${LAYER13} ${LAYER14} ${LAYER15} \ ${LAYER16} ${LAYER17} ${LAYER18} ${LAYER19} ${LAYER0COPY} -@ +\end{verbatim} \section{Cliques} The algebra code sometimes have circular references. The compiler can resolve these references directly if all of the required sources @@ -16965,7 +16965,7 @@ So the idea to remove the BOOTSTRAP code is to cluster the spad sources into "cliqueN.spad" files and feed them all to the compiler at once. <>= -CLIQUE1FILES = ${MID}/MYUP.spad ${MID}/MYEXPR.spad +CLIQUE1FILES = ${IN}/MYUP.spad ${IN}/MYEXPR.spad ${MID}/clique1.spad: ${CLIQUE1FILES} @echo cl1 making ${OUT}/MYUP.o from ${MID}/clique1.spad @@ -16977,6 +16977,8 @@ ${MID}/clique1.spad: ${CLIQUE1FILES} else \ echo ")co clique1.spad" | ${INTERPSYS} >${TMP}/trace ; \ fi ) + @ cp ${MID}/MYUP.nrlib/code.o ${OUT}/MYUP.o + @ cp ${MID}/MYEXPR.nrlib/code.o ${OUT}/MYEXPR.o @ Here we have the general case where two files are co-dependent, that is, @@ -16984,7 +16986,7 @@ PAFF and PAFFFF both have to be compiled together. They also have a set of prerequired files that must be loaded since they are not yet in the new database. <>= -CLIQUE2FILES = ${MID}/PAFF.spad ${MID}/PAFFFF.spad +CLIQUE2FILES = ${IN}/PAFF.spad ${IN}/PAFFFF.spad CLIQUE2DEPS = BLMETCT GPAFF PFORP PACOFF PROJPLPS PLACESPS NSDPS LOCPOWC \ DIV SETCATD PLACESC DIVCAT INFCLSPS INFCLCT DSTREE DSTRCAT \ PRSPCAT UTSZ PACFFC PACPERC PROJPL PLACES INFCLSPT PROJPL ICP diff --git a/src/axiom-website/patches.html b/src/axiom-website/patches.html index de946ea..e5f0fae 100644 --- a/src/axiom-website/patches.html +++ b/src/axiom-website/patches.html @@ -3848,5 +3848,7 @@ src/input/testpackage.input added
books/bookvol10.4 add PadeApproximants test
20120324.02.tpd.patch books/bookvol10.4 add getAncestors function to API domain
+20120325.01.tpd.patch +books/bookvol10.3,4,5 fix 7217
diff --git a/src/input/exampleagcode.input.pamphlet b/src/input/exampleagcode.input.pamphlet index 6239bf4..f765867 100644 --- a/src/input/exampleagcode.input.pamphlet +++ b/src/input/exampleagcode.input.pamphlet @@ -114,8 +114,8 @@ P1:= PAFFFF(K1,[X,Y,Z],BLQT) --R --R --R (3) ---R PackageForAlgebraicFunctionFieldOverFiniteField(PrimeField(2),[X,Y,Z],BlowUpWi ---R thQuadTrans) +--R PackageForAlgebraicFunctionFieldOverFiniteField(PrimeField(2),[X,Y,Z],BlowUpW +--R ithQuadTrans) --R Type: Domain --E 3 @@ -134,7 +134,7 @@ setCurve(C1)$P1 --R --R 5 2 3 4 --R (5) X + Y Z + Y Z ---R Type: DistributedMultivariatePolynomial([X,Y,Z],PrimeField(2)) +--R Type: DistributedMultivariatePolynomial([X,Y,Z],PrimeField(2)) --E 5 --S 6 of 19 @@ -143,7 +143,7 @@ plc3:= placesOfDegree(3)$P1 --R --R 2 3 2 3 --R (6) [[%D5:%D5 :1] ,[%D5:%D5 + 1:1] ] ---R Type: List(PlacesOverPseudoAlgebraicClosureOfFiniteField(PrimeField(2))) +--R Type: List(PlacesOverPseudoAlgebraicClosureOfFiniteField(PrimeField(2))) --E 6 -- %D4 is an elemenent created by the domain PACOFF: it is a root of the @@ -160,7 +160,7 @@ a:= elt( first % , 1 ) --R --R --R (7) %D5 ---R Type: PseudoAlgebraicClosureOfFiniteField(PrimeField(2)) +--R Type: PseudoAlgebraicClosureOfFiniteField(PrimeField(2)) --E 7 --S 8 of 19 @@ -177,7 +177,7 @@ a**3 + a**2 + 1 --R --R --R (9) 0 ---R Type: PseudoAlgebraicClosureOfFiniteField(PrimeField(2)) +--R Type: PseudoAlgebraicClosureOfFiniteField(PrimeField(2)) --E 9 -- As you can see, %D4 is the root of an irreducible poynomial of degree 3. @@ -192,7 +192,7 @@ D:= 2 * reduce(+,(plc3 :: List DIV PLACESPS PF 2)) --R --R 2 3 2 3 --R (10) 2 [%D5:%D5 :1] + 2 [%D5:%D5 + 1:1] ---R Type: Divisor(PlacesOverPseudoAlgebraicClosureOfFiniteField(PrimeField(2))) +--R Type: Divisor(PlacesOverPseudoAlgebraicClosureOfFiniteField(PrimeField(2))) --E 10 -- Now we compute a basis of L(D) @@ -280,7 +280,7 @@ plc1 := placesOfDegree(1)$P4 --R [%A :%A :1] , [%A :%A :1] , [%A :%A :1] , [%A :%A :1] , [%A :%A :1] , --R 1 1 1 --R [0:0:1] , [0:1:1] , %I3 ] ---IType: List(PlacesOverPseudoAlgebraicClosureOfFiniteField ... +--IType: List(PlacesOverPseudoAlgebraicClosureOfFiniteField(... --E 17 -- Now, we can construct the matrix of the AG-code, which code-words consist @@ -380,7 +380,7 @@ mG:= matrix [ [ eval( f, lB1.den, pl )$P4 for pl in plc1 ] for f in lB1.num ] --R %A , %A , %A , %A , 1, 1, %A , %A , %A , %A , %A , %A , %A , %A , --R 0, 0, 0] --R ] ---R Type: Matrix(FiniteFieldCyclicGroup(2,4)) +--R Type: Matrix(FiniteFieldCyclicGroup(2,4)) --E 18 -- The preceding matrix is the generator matrix of a [n,k,d]-code over GF(2^4) diff --git a/src/input/hyperell.input.pamphlet b/src/input/hyperell.input.pamphlet index 170a44c..ffe91ca 100644 --- a/src/input/hyperell.input.pamphlet +++ b/src/input/hyperell.input.pamphlet @@ -46,17 +46,17 @@ P:=PAFFFF( K, [x,y,z], BLQT) --R --R --R (4) ---R PackageForAlgebraicFunctionFieldOverFiniteField(PrimeField(1048583),[x,y,z],Bl ---R owUpWithQuadTrans) ---R--R--R Type: DistributedMultivariatePolynomial([x,y,z],PrimeField 1048583) +--R PackageForAlgebraicFunctionFieldOverFiniteField(PrimeField(1048583),[x,y,z],B +--R lowUpWithQuadTrans) +--R Type: Domain --E 4 --S 5 of 26 ProjPl := PROJPLPS PrimeField p --R --R - (4) - ProjectivePlaneOverPseudoAlgebraicClosureOfFiniteField(PrimeField(1048583)) +--R (5) +--R ProjectivePlaneOverPseudoAlgebraicClosureOfFiniteField(PrimeField(1048583)) --R Type: Domain --E 5 @@ -75,7 +75,7 @@ fh:R:= homogenize( f , 3 )$P --R --R 5 4 3 2 2 3 4 2 3 5 --R (7) 1048582x + 15x z + 1048498x z + 225x z + 1048309x z + y z + 120z ---R--R--RType: Divisor PlacesOverPseudoAlgebraicClosureOfFiniteField PrimeField 1048583 +--R Type: DistributedMultivariatePolynomial([x,y,z],PrimeField(1048583)) --E 7 --S 8 of 26 @@ -84,7 +84,7 @@ setCurve(fh)$P --R --R 5 4 3 2 2 3 4 2 3 5 --R (8) 1048582x + 15x z + 1048498x z + 225x z + 1048309x z + y z + 120z ---RType: ProjectivePlaneOverPseudoAlgebraicClosureOfFiniteField(PrimeField(1048583)) +--R Type: DistributedMultivariatePolynomial([x,y,z],PrimeField(1048583)) --E 8 --S 9 of 26 @@ -92,7 +92,7 @@ g:=genus()$P --R --R --R (9) 2 ---RType: ProjectivePlaneOverPseudoAlgebraicClosureOfFiniteField(PrimeField(1048583)) +--R Type: NonNegativeInteger --E 9 --S 10 of 26 @@ -101,7 +101,7 @@ divZ := intersectionDivisor(z)$P --R --R 1 --I (10) 5 %I7 ---RType: ProjectivePlaneOverPseudoAlgebraicClosureOfFiniteField(PrimeField(1048583)) +--RType: Divisor(PlacesOverPseudoAlgebraicClosureOfFiniteField(PrimeField(1048583))) --E 10 --S 11 of 26 @@ -110,7 +110,7 @@ pInf:= first supp divZ --R --R 1 --I (11) %I7 ---RType: ProjectivePlaneOverPseudoAlgebraicClosureOfFiniteField(PrimeField(1048583)) +--R Type: PlacesOverPseudoAlgebraicClosureOfFiniteField(PrimeField(1048583)) --E 11 --S 12 of 26 @@ -128,7 +128,7 @@ pl1:= first placesAbove( p1 )$P --R --R 1 --R (13) [1:0:1] ---R Type: DistributedMultivariatePolynomial([x,y,z],PrimeField(1048583)) +--R Type: PlacesOverPseudoAlgebraicClosureOfFiniteField(PrimeField(1048583)) --E 13 --S 14 of 26 @@ -137,7 +137,7 @@ p2:= projectivePoint( [2,0,1] :: List K )$ProjPl --R --R 1 --R (14) (2:0:1) ---RType: Divisor PlacesOverPseudoAlgebraicClosureOfFiniteField(PrimeField(1048583)) +--RType: ProjectivePlaneOverPseudoAlgebraicClosureOfFiniteField(PrimeField(1048583)) --E 14 --S 15 of 26 @@ -146,7 +146,7 @@ pl2:= first placesAbove( p2 )$P --R --R 1 --R (15) [2:0:1] ---RType: ProjectivePlaneOverPseudoAlgebraicClosureOfFiniteField(PrimeField(1048583)) +--R Type: PlacesOverPseudoAlgebraicClosureOfFiniteField(PrimeField(1048583)) --E 15 --S 16 of 26 @@ -155,7 +155,7 @@ p3:= projectivePoint( [3,0,1] :: List K )$ProjPl --R --R 1 --R (16) (3:0:1) ---R Type: PlacesOverPseudoAlgebraicClosureOfFiniteField(PrimeField(1048583)) +--RType: ProjectivePlaneOverPseudoAlgebraicClosureOfFiniteField(PrimeField(1048583)) --E 16 --S 17 of 26 @@ -164,7 +164,7 @@ pl3:= first placesAbove( p3 )$P --R --R 1 --R (17) [3:0:1] ---RType: Divisor PlacesOverPseudoAlgebraicClosureOfFiniteField(PrimeField(1048583)) +--R Type: PlacesOverPseudoAlgebraicClosureOfFiniteField(PrimeField(1048583)) --E 17 --S 18 of 26 @@ -173,7 +173,7 @@ p4:= projectivePoint( [4,0,1] :: List K )$ProjPl --R --R 1 --R (18) (4:0:1) ---R Type: DistributedMultivariatePolynomial([x,y,z],PrimeField(1048583)) +--RType: ProjectivePlaneOverPseudoAlgebraicClosureOfFiniteField(PrimeField(1048583)) --E 18 --S 19 of 26 @@ -182,7 +182,7 @@ pl4:= first placesAbove( p4 )$P --R --R 1 --R (19) [4:0:1] ---RType: Divisor PlacesOverPseudoAlgebraicClosureOfFiniteField(PrimeField 1048583)) +--R Type: PlacesOverPseudoAlgebraicClosureOfFiniteField(PrimeField(1048583)) --E 19 --S 20 of 26 @@ -200,7 +200,7 @@ pl5:= first placesAbove( p5 )$P --R --R 1 --R (21) [5:0:1] ---R Type: PlacesOverPseudoAlgebraicClosureOfFiniteField(PrimeField(1048583)) +--R Type: PlacesOverPseudoAlgebraicClosureOfFiniteField(PrimeField(1048583)) --E 21 --S 22 of 26 @@ -209,7 +209,7 @@ D:= pl1+pl2+ 3*pl3 - 5* pInf --R --R 1 1 1 1 --I (22) [1:0:1] + [2:0:1] + 3 [3:0:1] - 5 %I7 ---RType: Divisor PlacesOverPseudoAlgebraicClosureOfFiniteField(PrimeField(1048583)) +--RType: Divisor(PlacesOverPseudoAlgebraicClosureOfFiniteField(PrimeField(1048583))) --E 22 --S 23 of 26 @@ -225,7 +225,7 @@ lb:= lBasis( D + g*pInf )$P --R (23) --R 2 3 2 2 2 3 --R [num= [873819x y z + y z ],den= 873819x + x z + 174762x z + 87382y z + z ] ---IType: Record(num: List DistributedMultivariatePolynomial(... +--IType: Record(num: List(DistributedMultivariatePolynomial(... --E 23 --S 24 of 26 @@ -234,7 +234,7 @@ g1:= first lb.num --R --R 2 --R (24) 873819x y z + y z ---R Type: DistributedMultivariatePolynomial([x,y,z],PrimeField(1048583)) +--R Type: DistributedMultivariatePolynomial([x,y,z],PrimeField(1048583)) --E 24 --S 25 of 26 @@ -243,7 +243,7 @@ g0:= lb.den --R --R 3 2 2 2 3 --R (25) 873819x + x z + 174762x z + 87382y z + z ---R Type: DistributedMultivariatePolynomial([x,y,z],PrimeField(1048583)) +--R Type: DistributedMultivariatePolynomial([x,y,z],PrimeField(1048583)) --E 25 -- Voici le diviseur equivalent a D ayant un diviseur des zeros (partie effective)