diff --git a/books/bookvol10.3.pamphlet b/books/bookvol10.3.pamphlet index 4a12546..7d78657 100644 --- a/books/bookvol10.3.pamphlet +++ b/books/bookvol10.3.pamphlet @@ -32337,6 +32337,151 @@ GraphImage (): Exports == Implementation where @ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +\section{domain GOPT GuessOption} +\pagehead{GuessOption}{GOPT} +\pagepic{ps/v103guessoption.ps}{GOPT}{1.00} +<>= +)abbrev domain GOPT GuessOption +++ Author: Martin Rubey +++ Description: GuessOption is a domain whose elements are various options used +++ by \spadtype{Guess}. +GuessOption(): Exports == Implementation where + + Exports == SetCategory with + + maxDerivative: Integer -> % + ++ maxDerivative(d) specifies the maximum derivative in an algebraic + ++ differential equation. maxDerivative(-1) specifies that the maximum + ++ derivative can be arbitrary. This option is expressed in the form + ++ \spad{maxDerivative == d}. + + maxShift: Integer -> % + ++ maxShift(d) specifies the maximum shift in a recurrence + ++ equation. maxShift(-1) specifies that the maximum shift can be + ++ arbitrary. This option is expressed in the form \spad{maxShift == d}. + + maxPower: Integer -> % + ++ maxPower(d) specifies the maximum degree in an algebraic differential + ++ equation. For example, the degree of (f'')^3 f' is 4. maxPower(-1) + ++ specifies that the maximum exponent can be arbitrary. This option is + ++ expressed in the form \spad{maxPower == d}. + + homogeneous: Boolean -> % + ++ homogeneous(d) specifies whether we allow only homogeneous algebraic + ++ differential equations. This option is expressed in the form + ++ \spad{homogeneous == d}. + + maxLevel: Integer -> % + ++ maxLevel(d) specifies the maximum number of recursion levels operators + ++ guessProduct and guessSum will be applied. maxLevel(-1) specifies + ++ that all levels are tried. This option is expressed in the form + ++ \spad{maxLevel == d}. + + maxDegree: Integer -> % + ++ maxDegree(d) specifies the maximum degree of the coefficient + ++ polynomials in an algebraic differential equation or a recursion with + ++ polynomial coefficients. For rational functions with an exponential + ++ term, \spad{maxDegree} bounds the degree of the denominator + ++ polynomial. + ++ maxDegree(-1) specifies that the maximum + ++ degree can be arbitrary. This option is expressed in the form + ++ \spad{maxDegree == d}. + + allDegrees: Boolean -> % + ++ allDegrees(d) specifies whether all possibilities of the degree vector + ++ - taking into account maxDegree - should be tried. This is mainly + ++ interesting for rational interpolation. This option is expressed in + ++ the form \spad{allDegrees == d}. + + safety: NonNegativeInteger -> % + ++ safety(d) specifies the number of values reserved for testing any + ++ solutions found. This option is expressed in the form \spad{safety == + ++ d}. + + one: Boolean -> % + ++ one(d) specifies whether we are happy with one solution. This option + ++ is expressed in the form \spad{one == d}. + + debug: Boolean -> % + ++ debug(d) specifies whether we want additional output on the + ++ progress. This option is expressed in the form \spad{debug == d}. + + functionName: Symbol -> % + ++ functionName(d) specifies the name of the function given by the + ++ algebraic differential equation or recurrence. This option is + ++ expressed in the form \spad{functionName == d}. + + variableName: Symbol -> % + ++ variableName(d) specifies the variable used in by the algebraic + ++ differential equation. This option is expressed in the form + ++ \spad{variableName == d}. + + indexName: Symbol -> % + ++ indexName(d) specifies the index variable used for the formulas. This + ++ option is expressed in the form \spad{indexName == d}. + + displayAsGF: Boolean -> % + ++ displayAsGF(d) specifies whether the result is a generating function + ++ or a recurrence. This option should not be set by the user, but rather + ++ by the HP-specification. + + option : (List %, Symbol) -> Union(Any, "failed") + ++ option() is not to be used at the top level; + ++ option determines internally which drawing options are indicated in + ++ a draw command. + + option?: (List %, Symbol) -> Boolean + ++ option?() is not to be used at the top level; + ++ option? internally returns true for drawing options which are + ++ indicated in a draw command, or false for those which are not. + + checkOptions: List % -> Void + ++ checkOptions checks whether an option is given twice + + Implementation ==> add + import AnyFunctions1(Boolean) + import AnyFunctions1(Symbol) + import AnyFunctions1(Integer) + import AnyFunctions1(NonNegativeInteger) + + Rep := Record(keyword: Symbol, value: Any) + + maxLevel d == ["maxLevel"::Symbol, d::Any] + maxDerivative d == ["maxDerivative"::Symbol, d::Any] + maxShift d == maxDerivative d + maxDegree d == ["maxDegree"::Symbol, d::Any] + allDegrees d == ["allDegrees"::Symbol, d::Any] + maxPower d == ["maxPower"::Symbol, d::Any] + safety d == ["safety"::Symbol, d::Any] + homogeneous d == ["homogeneous"::Symbol, d::Any] + debug d == ["debug"::Symbol, d::Any] + one d == ["one"::Symbol, d::Any] + functionName d == ["functionName"::Symbol, d::Any] + variableName d == ["variableName"::Symbol, d::Any] + indexName d == ["indexName"::Symbol, d::Any] + displayAsGF d == ["displayAsGF"::Symbol, d::Any] + + coerce(x:%):OutputForm == x.keyword::OutputForm = x.value::OutputForm + x:% = y:% == x.keyword = y.keyword and x.value = y.value + + option?(l, s) == + for x in l repeat + x.keyword = s => return true + false + + option(l, s) == + for x in l repeat + x.keyword = s => return(x.value) + "failed" + + checkOptions l == + if not empty? l then + if find((first l).keyword = #1.keyword, rest l) case "failed" + then checkOptions rest l + else error "GuessOption: Option specified twice" + +@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Chapter H} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{domain HASHTBL HashTable} @@ -69807,11 +69952,13 @@ SparseUnivariatePolynomial(R:Ring): UnivariatePolynomialCategory(R) with @ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{domain SUPEXPR SparseUnivariatePolynomialExpressions} + This domain is a hack, in some sense. What I'd really like to do - automatically - is to provide all operations supported by the coefficient domain, as long as the polynomials can be retracted to that domain, i.e., as long as they are just constants. I don't see another way to do this, unfortunately. + \pagehead{SparseUnivariatePolynomialExpressions}{SUPEXPR} \pagepic{ps/v103sparseunivariatepolynomialexpressions.ps}{SUPEXPR}{1.00} <>= @@ -69826,6 +69973,11 @@ SparseUnivariatePolynomialExpressions(R: Ring): Exports == Implementation where Implementation == SparseUnivariatePolynomial R add if R has TranscendentalFunctionCategory then + log(p: %): % == + ground? p => coerce log ground p + output(hconcat("log p for p= ", p::OutputForm))$OutputPackage + error "SUPTRAFUN: log only defined for elements of the coefficient ring" + exp(p: %): % == ground? p => coerce exp ground p output(hconcat("exp p for p= ", p::OutputForm))$OutputPackage @@ -80133,6 +80285,16 @@ TwoDimensionalViewport ():Exports == Implementation where %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Chapter U} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +\section{domain UFPS UnivariateFormalPowerSeries} +\pagehead{UnivariateFormalPowerSeries}{UFPS} +\pagepic{ps/v103univariateformalpowerseries.ps}{UFPS}{1.00} +<>= +)abbrev domain UFPS UnivariateFormalPowerSeries +UnivariateFormalPowerSeries(Coef: Ring) == + UnivariateTaylorSeries(Coef, 'x, 0$Coef) + +@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{domain ULS UnivariateLaurentSeries} \pagehead{UnivariateLaurentSeries}{ULS} \pagepic{ps/v103univariatelaurentseries.ps}{ULS}{1.00} @@ -91111,6 +91273,7 @@ Note that this code is not included in the generated catdef.spad file. <> <> <> +<> <> <> @@ -91321,6 +91484,7 @@ Note that this code is not included in the generated catdef.spad file. <> <> +<> <> <> <> diff --git a/changelog b/changelog index 7bc32ab..1e0a9cf 100644 --- a/changelog +++ b/changelog @@ -1,3 +1,14 @@ +20081217 tpd src/axiom-website/patches.html 20081217.02.tpd.patch +20081217 tpd src/input/ndftip.input fix regression on map operation count +20081217 tpd src/algebra/Makefile remove wtpol.spad +20081217 tpd src/algebra/mantepse.spad move domains to bookvol10.3 +20081217 tpd src/algebra/ssolve.spad move domain to bookvol10.3 +20081217 tpd src/algebra/ssolve.spad add guess package code +20081217 tpd src/algebra/rec.spad add guess package code +20081217 tpd src/algebra/mantepse.spad add guess package code +20081217 tpd src/algebra/Makefile add guess package code +20081217 tpd src/algebra/fffg.spad add guess package code +20081217 tpd src/algebra/exposed.lsp add guess package code 20081217 tpd src/axiom-website/patches.html 20081217.01.tpd.patch 20081217 tpd src/algebra/padic.spad removed, move domains to bookvol10.3 20081216 tpd src/axiom-website/patches.html 20081216.04.tpd.patch diff --git a/src/algebra/Makefile.pamphlet b/src/algebra/Makefile.pamphlet index 22de514..ce9ddba 100644 --- a/src/algebra/Makefile.pamphlet +++ b/src/algebra/Makefile.pamphlet @@ -422,7 +422,9 @@ xlpoly.spad.pamphlet (MAGMA LWORD LIECAT FLALG XEXPPKG LPOLY PBWLB XPBWPOLY LAYER11=\ ${OUT}/APPLYORE.o ${OUT}/ARRAY1.o ${OUT}/ARRAY12.o ${OUT}/ARRAY2.o \ ${OUT}/ASTACK.o ${OUT}/BTAGG.o ${OUT}/BTAGG-.o ${OUT}/COMBINAT.o \ - ${OUT}/CSTTOOLS.o ${OUT}/D01FCFA.o ${OUT}/E04MBFA.o ${OUT}/FARRAY.o \ + ${OUT}/CSTTOOLS.o ${OUT}/D01FCFA.o ${OUT}/E04MBFA.o \ + ${OUT}/FAMR2.o \ + ${OUT}/FARRAY.o \ ${OUT}/FLALG.o ${OUT}/GALUTIL.o ${OUT}/HEAP.o ${OUT}/IARRAY1.o \ ${OUT}/IARRAY2.o ${OUT}/IFARRAY.o ${OUT}/INTCAT.o ${OUT}/INTHEORY.o \ ${OUT}/IRREDFFX.o ${OUT}/LFCAT.o ${OUT}/LODOCAT.o ${OUT}/LODOCAT-.o \ @@ -579,7 +581,8 @@ LAYER14=\ ${OUT}/D01AQFA.o ${OUT}/EMR.o ${OUT}/EQ.o ${OUT}/ERROR.o \ ${OUT}/EVALCYC.o ${OUT}/E04DGFA.o ${OUT}/E04FDFA.o ${OUT}/E04GCFA.o \ ${OUT}/E04JAFA.o ${OUT}/FACUTIL.o ${OUT}/FF.o ${OUT}/FFCG.o \ - ${OUT}/FFCGX.o ${OUT}/FFHOM.o ${OUT}/FFNB.o ${OUT}/FFNBX.o \ + ${OUT}/FFCGX.o ${OUT}/FFFG.o \ + ${OUT}/FFHOM.o ${OUT}/FFNB.o ${OUT}/FFNBX.o \ ${OUT}/FFPOLY.o ${OUT}/FFX.o ${OUT}/FFSLPE.o ${OUT}/FGLMICPK.o \ ${OUT}/FILE.o ${OUT}/FINAALG.o ${OUT}/FINAALG-.o ${OUT}/FINRALG.o \ ${OUT}/FINRALG-.o ${OUT}/FFF.o ${OUT}/FLOATRP.o ${OUT}/FNAME.o \ @@ -616,7 +619,8 @@ LAYER14=\ ${OUT}/SCPKG.o ${OUT}/SHDP.o ${OUT}/SHP.o ${OUT}/SIGNRF.o \ ${OUT}/SMITH.o ${OUT}/SMP.o ${OUT}/SMTS.o ${OUT}/SOLVEFOR.o \ ${OUT}/SPLTREE.o ${OUT}/STINPROD.o ${OUT}/STTFNC.o ${OUT}/SUBRESP.o \ - ${OUT}/SUMRF.o ${OUT}/SUP.o ${OUT}/SUPFRACF.o ${OUT}/TANEXP.o \ + ${OUT}/SUMRF.o ${OUT}/SUP.o ${OUT}/SUPEXPR.o \ + ${OUT}/SUPFRACF.o ${OUT}/TANEXP.o \ ${OUT}/TEMUTL.o ${OUT}/TEX.o ${OUT}/TEXTFILE.o ${OUT}/TREE.o \ ${OUT}/TWOFACT.o ${OUT}/UNIFACT.o ${OUT}/UP.o ${OUT}/UPCDEN.o \ ${OUT}/UPDECOMP.o ${OUT}/UPDIVP.o ${OUT}/UPMP.o ${OUT}/UPOLYC2.o \ @@ -635,6 +639,7 @@ plot.spad.pamphlet (PLOT PLOT1) LAYER15=\ ${OUT}/DIAGG.o ${OUT}/DIAGG-.o ${OUT}/DSMP.o ${OUT}/EXPUPXS.o \ + ${OUT}/FFFGF.o \ ${OUT}/FRAMALG.o ${OUT}/FRAMALG-.o ${OUT}/MDAGG.o ${OUT}/ODPOL.o \ ${OUT}/PLOT.o ${OUT}/RMCAT2.o ${OUT}/ROIRC.o ${OUT}/SDPOL.o \ ${OUT}/SMATCAT.o ${OUT}/SMATCAT-.o ${OUT}/TUBETOOL.o ${OUT}/UPXSCCA.o \ @@ -762,7 +767,8 @@ LAYER19=\ ${OUT}/FFCAT.o ${OUT}/FFCAT-.o ${OUT}/FFCGP.o ${OUT}/FFNBP.o \ ${OUT}/FFP.o ${OUT}/FLOAT.o ${OUT}/FPARFRAC.o ${OUT}/FR.o \ ${OUT}/FRNAALG.o ${OUT}/FRNAALG-.o ${OUT}/FS.o ${OUT}/FS-.o \ - ${OUT}/FST.o ${OUT}/FUNCTION.o ${OUT}/GDMP.o ${OUT}/HACKPI.o \ + ${OUT}/FST.o ${OUT}/FUNCTION.o ${OUT}/GDMP.o ${OUT}/GOPT.o \ + ${OUT}/HACKPI.o \ ${OUT}/IDEAL.o ${OUT}/INFORM.o ${OUT}/INFORM1.o ${OUT}/IPRNTPK.o \ ${OUT}/IR.o ${OUT}/ISUPS.o ${OUT}/KERNEL.o ${OUT}/LIB.o \ ${OUT}/LMDICT.o ${OUT}/LODOOPS.o ${OUT}/MATRIX.o ${OUT}/MKFLCFN.o \ @@ -865,6 +871,7 @@ LAYER20=\ ${OUT}/FORDER.o ${OUT}/FORTRAN.o ${OUT}/FSRED.o ${OUT}/FSUPFACT.o \ ${OUT}/FRNAAF2.o ${OUT}/FSPECF.o ${OUT}/FS2.o ${OUT}/FS2UPS.o \ ${OUT}/GAUSSFAC.o ${OUT}/GCNAALG.o ${OUT}/GENUFACT.o ${OUT}/GENUPS.o \ + ${OUT}/GOPT0.o \ ${OUT}/GTSET.o ${OUT}/GPOLSET.o ${OUT}/IAN.o ${OUT}/INEP.o \ ${OUT}/INFPROD0.o ${OUT}/INFSP.o ${OUT}/INPRODFF.o ${OUT}/INPRODPF.o \ ${OUT}/INTAF.o ${OUT}/INTALG.o ${OUT}/INTEF.o ${OUT}/INTG0.o \ @@ -881,9 +888,10 @@ LAYER20=\ ${OUT}/RDEEF.o ${OUT}/RDEEFS.o ${OUT}/RDIV.o ${OUT}/RSETCAT.o \ ${OUT}/RSETCAT-.o ${OUT}/RULE.o ${OUT}/RULESET.o ${OUT}/SIMPAN.o \ ${OUT}/SFORT.o ${OUT}/SOLVESER.o ${OUT}/SUMFS.o ${OUT}/SUTS.o \ - ${OUT}/TOOLSIGN.o ${OUT}/TRIGMNIP.o ${OUT}/TRMANIP.o ${OUT}/ULSCCAT.o \ + ${OUT}/TOOLSIGN.o ${OUT}/TRIGMNIP.o ${OUT}/TRMANIP.o ${OUT}/UFPS.o \ + ${OUT}/ULSCCAT.o \ ${OUT}/ULSCCAT-.o ${OUT}/UPXSSING.o ${OUT}/UTSODE.o ${OUT}/UTSODETL.o \ - ${OUT}/UTS2.o ${OUT}/WUTSET.o \ + ${OUT}/UTSSOL.o ${OUT}/UTS2.o ${OUT}/WUTSET.o \ layer20done @ @@ -918,12 +926,15 @@ taylor.spad.pamphlet (ITAYLOR UTS UTS2) LAYER21=\ ${OUT}/DEFINTEF.o ${OUT}/DFINTTLS.o ${OUT}/DEFINTRF.o ${OUT}/D01TRNS.o \ - ${OUT}/EFULS.o ${OUT}/ESCONT.o ${OUT}/EXPR.o ${OUT}/EXPR2UPS.o \ + ${OUT}/EFULS.o ${OUT}/ESCONT.o ${OUT}/EXPR.o ${OUT}/EXPRSOL.o \ + ${OUT}/EXPR2UPS.o \ ${OUT}/FDIV.o ${OUT}/FSCINT.o ${OUT}/FSINT.o ${OUT}/FS2EXPXP.o \ ${OUT}/GSERIES.o ${OUT}/HELLFDIV.o ${OUT}/INVLAPLA.o ${OUT}/IR2F.o \ ${OUT}/IRRF2F.o ${OUT}/LAPLACE.o ${OUT}/LIMITPS.o ${OUT}/LODEEF.o \ - ${OUT}/NODE1.o ${OUT}/ODECONST.o ${OUT}/ODEINT.o ${OUT}/REP.o \ - ${OUT}/SOLVERAD.o ${OUT}/SULS.o ${OUT}/SUPXS.o ${OUT}/ULS.o \ + ${OUT}/NODE1.o ${OUT}/ODECONST.o ${OUT}/ODEINT.o ${OUT}/RECOP.o \ + ${OUT}/REP.o \ + ${OUT}/SOLVERAD.o ${OUT}/SULS.o ${OUT}/SUPXS.o ${OUT}/UFPS1.o \ + ${OUT}/ULS.o \ ${OUT}/ULSCONS.o ${OUT}/UPXS.o ${OUT}/UPXSCONS.o ${OUT}/UTS.o \ layer21done @@ -949,6 +960,7 @@ zerodim.spad.pamphlet (RGCHAIN ZDSOLVE) LAYER22=\ ${OUT}/ASP29.o ${OUT}/COMBF.o ${OUT}/D01AGNT.o ${OUT}/FSPRMELT.o \ + ${OUT}/GUESS.o \ ${OUT}/INBFF.o ${OUT}/LODO.o ${OUT}/LODO1.o ${OUT}/LODO2.o \ ${OUT}/NTSCAT.o ${OUT}/REGSET.o ${OUT}/RGCHAIN.o ${OUT}/RSETGCD.o \ ${OUT}/RSDCMPK.o ${OUT}/SFRTCAT.o ${OUT}/SIGNEF.o ${OUT}/SNTSCAT.o \ @@ -975,12 +987,21 @@ zerodim.spad.pamphlet (LEXTRIPK IRURPK RURPK) <>= LAYER23=\ - ${OUT}/CPIMA.o ${OUT}/IRURPK.o ${OUT}/LAZM3PK.o ${OUT}/LEXTRIPK.o \ + ${OUT}/CPIMA.o ${OUT}/GUESSAN.o ${OUT}/GUESSINT.o \ + ${OUT}/GUESSF1.o ${OUT}/GUESSP.o ${OUT}/GUESSUP.o \ + ${OUT}/IRURPK.o ${OUT}/LAZM3PK.o ${OUT}/LEXTRIPK.o \ ${OUT}/NORMPK.o ${OUT}/QCMPACK.o ${OUT}/RURPK.o ${OUT}/SFRGCD.o \ ${OUT}/SFQCMPK.o ${OUT}/INTRVL.o ${OUT}/ODEEF.o \ layer23done @ +<>= + +LAYER24=\ + ${OUT}/GUESSF.o \ + layer24done + +@ \subsection{User Layer for newly added algebra} Rather than classify newly created algebra into the existing type lattice we add it here. @@ -998,9 +1019,184 @@ ORDER=\ ${LAYER4} ${LAYER5} ${LAYER6} ${LAYER7} ${LAYER8} ${LAYER9} \ ${LAYER10} ${LAYER11} ${LAYER12} ${LAYER13} ${LAYER14} ${LAYER15} \ ${LAYER16} ${LAYER17} ${LAYER18} ${LAYER19} ${LAYER20} ${LAYER21} \ - ${LAYER22} ${LAYER23} ${USERLAYER} ${LAYER0COPY} + ${LAYER22} ${LAYER23} ${LAYER24} ${USERLAYER} ${LAYER0COPY} @ +\section{New Algebra dependencies} + +New algebra files can depend on new algebra files. Since these files +are not part of the default database we need to be explicit about the +new files being loaded with the library command. + +<>= + +FFFGFDEPS = FAMR2 FFFG + +${MID}/FFFGF.nrlib/code.o: ${MID}/FFFGF.spad + @echo S1 making ${MID}/FFFGF.nrlib/code.o from ${MID}/FFFGF.spad + @ (cd ${MID} ; \ + if [ -z "${NOISE}" ] ; then \ + echo -e ")lib ${FFFGFDEPS} \n )co FFFGF.spad" | ${INTERPSYS} ; \ + else \ + echo -e ")lib ${FFFGFDEPS} \n )co FFFGF.spad" \ + | ${INTERPSYS} >${TMP}/trace ; \ + fi ) + +UTSSOLDEPS = SUPEXPR + +${MID}/UTSSOL.nrlib/code.o: ${MID}/UTSSOL.spad + @echo S1 making ${MID}/UTSSOL.nrlib/code.o from ${MID}/UTSSOL.spad + @ (cd ${MID} ; \ + if [ -z "${NOISE}" ] ; then \ + echo -e ")lib ${UTSSOLDEPS} \n )co UTSSOL.spad" | ${INTERPSYS} ; \ + else \ + echo -e ")lib ${UTSSOLDEPS} \n )co UTSSOL.spad" \ + | ${INTERPSYS} >${TMP}/trace ; \ + fi ) + +EXPRSOLDEPS = SUPEXPR UTSSOL + +${MID}/EXPRSOL.nrlib/code.o: ${MID}/EXPRSOL.spad + @echo S1 making ${MID}/EXPRSOL.nrlib/code.o from ${MID}/EXPRSOL.spad + @ (cd ${MID} ; \ + if [ -z "${NOISE}" ] ; then \ + echo -e ")lib ${EXPRSOLDEPS} \n )co EXPRSOL.spad" | ${INTERPSYS} ; \ + else \ + echo -e ")lib ${EXPRSOLDEPS} \n )co EXPRSOL.spad" \ + | ${INTERPSYS} >${TMP}/trace ; \ + fi ) + +GOPT0DEPS = GOPT + +${MID}/GOPT0.nrlib/code.o: ${MID}/GOPT0.spad + @echo S1 making ${MID}/GOPT0.nrlib/code.o from ${MID}/GOPT0.spad + @ (cd ${MID} ; \ + if [ -z "${NOISE}" ] ; then \ + echo -e ")lib ${GOPT0DEPS} \n )co GOPT0.spad" | ${INTERPSYS} ; \ + else \ + echo -e ")lib ${GOPT0DEPS} \n )co GOPT0.spad" \ + | ${INTERPSYS} >${TMP}/trace ; \ + fi ) + +RECOPDEPS = SUPEXPR UTSSOL EXPRSOL UFPS + +${MID}/RECOP.nrlib/code.o: ${MID}/RECOP.spad + @echo S1 making ${MID}/RECOP.nrlib/code.o from ${MID}/RECOP.spad + @ (cd ${MID} ; \ + if [ -z "${NOISE}" ] ; then \ + echo -e ")lib ${RECOPDEPS} \n )co RECOP.spad" | ${INTERPSYS} ; \ + else \ + echo -e ")lib ${RECOPDEPS} \n )co RECOP.spad" \ + | ${INTERPSYS} >${TMP}/trace ; \ + fi ) + +UFPS1DEPS = UFPS + +${MID}/UFPS1.nrlib/code.o: ${MID}/UFPS1.spad + @echo S1 making ${MID}/UFPS1.nrlib/code.o from ${MID}/UFPS1.spad + @ (cd ${MID} ; \ + if [ -z "${NOISE}" ] ; then \ + echo -e ")lib ${UFPS1DEPS} \n )co UFPS1.spad" | ${INTERPSYS} ; \ + else \ + echo -e ")lib ${UFPS1DEPS} \n )co UFPS1.spad" \ + | ${INTERPSYS} >${TMP}/trace ; \ + fi ) + +GUESSDEPS = NEWTON FAMR2 FFFG FFFGF SUPEXPR UTSSOL EXPRSOL GOPT GOPT0 \ + UFPS RECOP UFPS1 + +${MID}/GUESS.nrlib/code.o: ${MID}/GUESS.spad + @echo S1 making ${MID}/GUESS.nrlib/code.o from ${MID}/GUESS.spad + @ (cd ${MID} ; \ + if [ -z "${NOISE}" ] ; then \ + echo -e ")lib ${GUESSDEPS} \n )co GUESS.spad" | ${INTERPSYS} ; \ + else \ + echo -e ")lib ${GUESSDEPS} \n )co GUESS.spad" \ + | ${INTERPSYS} >${TMP}/trace ; \ + fi ) + +GUESSINTDEPS = NEWTON FAMR2 FFFG FFFGF SUPEXPR UTSSOL EXPRSOL GOPT GOPT0 \ + UFPS RECOP UFPS1 GUESS + +${MID}/GUESSINT.nrlib/code.o: ${MID}/GUESSINT.spad + @echo S1 making ${MID}/GUESSINT.nrlib/code.o from ${MID}/GUESSINT.spad + @ (cd ${MID} ; \ + if [ -z "${NOISE}" ] ; then \ + echo -e ")lib ${GUESSINTDEPS} \n )co GUESSINT.spad" | ${INTERPSYS} ; \ + else \ + echo -e ")lib ${GUESSINTDEPS} \n )co GUESSINT.spad" \ + | ${INTERPSYS} >${TMP}/trace ; \ + fi ) + +GUESSF1DEPS = NEWTON FAMR2 FFFG FFFGF SUPEXPR UTSSOL EXPRSOL GOPT GOPT0 \ + UFPS RECOP UFPS1 GUESS + +${MID}/GUESSF1.nrlib/code.o: ${MID}/GUESSF1.spad + @echo S1 making ${MID}/GUESSF1.nrlib/code.o from ${MID}/GUESSF1.spad + @ (cd ${MID} ; \ + if [ -z "${NOISE}" ] ; then \ + echo -e ")lib ${GUESSF1DEPS} \n )co GUESSF1.spad" | ${INTERPSYS} ; \ + else \ + echo -e ")lib ${GUESSF1DEPS} \n )co GUESSF1.spad" \ + | ${INTERPSYS} >${TMP}/trace ; \ + fi ) + +GUESSFDEPS = NEWTON FAMR2 FFFG FFFGF SUPEXPR UTSSOL EXPRSOL GOPT GOPT0 \ + UFPS RECOP UFPS1 GUESS GUESSF1 + +${MID}/GUESSF.nrlib/code.o: ${MID}/GUESSF.spad + @echo S1 making ${MID}/GUESSF.nrlib/code.o from ${MID}/GUESSF.spad + @ (cd ${MID} ; \ + if [ -z "${NOISE}" ] ; then \ + echo -e ")lib ${GUESSFDEPS} \n )co GUESSF.spad" | ${INTERPSYS} ; \ + else \ + echo -e ")lib ${GUESSFDEPS} \n )co GUESSF.spad" \ + | ${INTERPSYS} >${TMP}/trace ; \ + fi ) + +GUESSPDEPS = NEWTON FAMR2 FFFG FFFGF SUPEXPR UTSSOL EXPRSOL GOPT GOPT0 \ + UFPS RECOP UFPS1 GUESS + +${MID}/GUESSP.nrlib/code.o: ${MID}/GUESSP.spad + @echo S1 making ${MID}/GUESSP.nrlib/code.o from ${MID}/GUESSP.spad + @ (cd ${MID} ; \ + if [ -z "${NOISE}" ] ; then \ + echo -e ")lib ${GUESSPDEPS} \n )co GUESSP.spad" | ${INTERPSYS} ; \ + else \ + echo -e ")lib ${GUESSPDEPS} \n )co GUESSP.spad" \ + | ${INTERPSYS} >${TMP}/trace ; \ + fi ) + +GUESSANDEPS = NEWTON FAMR2 FFFG FFFGF SUPEXPR UTSSOL EXPRSOL GOPT GOPT0 \ + UFPS RECOP UFPS1 GUESS + +${MID}/GUESSAN.nrlib/code.o: ${MID}/GUESSAN.spad + @echo S1 making ${MID}/GUESSAN.nrlib/code.o from ${MID}/GUESSAN.spad + @ (cd ${MID} ; \ + if [ -z "${NOISE}" ] ; then \ + echo -e ")lib ${GUESSANDEPS} \n )co GUESSAN.spad" | ${INTERPSYS} ; \ + else \ + echo -e ")lib ${GUESSANDEPS} \n )co GUESSAN.spad" \ + | ${INTERPSYS} >${TMP}/trace ; \ + fi ) + +GUESSUPDEPS = NEWTON FAMR2 FFFG FFFGF SUPEXPR UTSSOL EXPRSOL GOPT GOPT0 \ + UFPS RECOP UFPS1 GUESS + +${MID}/GUESSUP.nrlib/code.o: ${MID}/GUESSUP.spad + @echo S1 making ${MID}/GUESSUP.nrlib/code.o from ${MID}/GUESSUP.spad + @ (cd ${MID} ; \ + if [ -z "${NOISE}" ] ; then \ + echo -e ")lib ${GUESSUPDEPS} \n )co GUESSUP.spad" | ${INTERPSYS} ; \ + else \ + echo -e ")lib ${GUESSUPDEPS} \n )co GUESSUP.spad" \ + | ${INTERPSYS} >${TMP}/trace ; \ + fi ) + + + +@ + \section{Broken Files} These files are Aldor files \begin{verbatim} @@ -1073,6 +1269,16 @@ files. INTERPSYS=${OBJ}/${SYS}/bin/interpsys @ +\subsection{The shell variable} +We use the ``-e'' flag to echo which is not supported by the ``sh'' +shell but is supported by ``bash''. The ``-e'' flag to echo causes +it to interpret special characters, in our case newlines. + +<>= + +SHELL=bash + +@ \subsection{The SPADFILES list} Note that we have excluded {\bf mlift.spad.jhd} from this list. We need to figure out which mlift.spad to keep. @@ -1204,7 +1410,7 @@ SPADFILES= \ ${OUTSRC}/viewdef.spad ${OUTSRC}/viewpack.spad \ ${OUTSRC}/void.spad \ ${OUTSRC}/weier.spad \ - ${OUTSRC}/xlpoly.spad \ + ${OUTSRC}/xlpoly.spad \ ${OUTSRC}/ystream.spad \ ${OUTSRC}/zerodim.spad @@ -1356,7 +1562,7 @@ DOCFILES= \ ${DOC}/viewdef.spad.dvi ${DOC}/viewpack.spad.dvi \ ${DOC}/void.spad.dvi \ ${DOC}/weier.spad.dvi \ - ${DOC}/xlpoly.spad.dvi \ + ${DOC}/xlpoly.spad.dvi \ ${DOC}/ystream.spad.dvi \ ${DOC}/zerodim.spad.dvi @@ -1729,122 +1935,127 @@ layer0copy: layer0done: @ echo ================================== - @ echo === layer 0 of 23 complete ====== + @ echo === layer 0 of 24 complete ====== @ echo ================================== layer1done: @ echo ================================== - @ echo === layer 1 of 23 complete ====== + @ echo === layer 1 of 24 complete ====== @ echo ================================== layer2done: @ echo ================================== - @ echo === layer 2 of 23 complete ====== + @ echo === layer 2 of 24 complete ====== @ echo ================================== layer3done: @ echo ================================== - @ echo === layer 3 of 23 complete ====== + @ echo === layer 3 of 24 complete ====== @ echo ================================== layer4done: @ echo ================================== - @ echo === layer 4 of 23 complete ====== + @ echo === layer 4 of 24 complete ====== @ echo ================================== layer5done: @ echo ================================== - @ echo === layer 5 of 23 complete ====== + @ echo === layer 5 of 24 complete ====== @ echo ================================== layer6done: @ echo ================================== - @ echo === layer 6 of 23 complete ====== + @ echo === layer 6 of 24 complete ====== @ echo ================================== layer7done: @ echo ================================== - @ echo === layer 7 of 23 complete ====== + @ echo === layer 7 of 24 complete ====== @ echo ================================== layer8done: @ echo ================================== - @ echo === layer 8 of 23 complete ====== + @ echo === layer 8 of 24 complete ====== @ echo ================================== layer9done: @ echo ================================== - @ echo === layer 9 of 23 complete ====== + @ echo === layer 9 of 24 complete ====== @ echo ================================== layer10done: @ echo ================================== - @ echo === layer 10 of 23 complete ====== + @ echo === layer 10 of 24 complete ====== @ echo ================================== layer11done: @ echo ================================== - @ echo === layer 11 of 23 complete ====== + @ echo === layer 11 of 24 complete ====== @ echo ================================== layer12done: @ echo ================================== - @ echo === layer 12 of 23 complete ====== + @ echo === layer 12 of 24 complete ====== @ echo ================================== layer13done: @ echo ================================== - @ echo === layer 13 of 23 complete ====== + @ echo === layer 13 of 24 complete ====== @ echo ================================== layer14done: @ echo ================================== - @ echo === layer 14 of 23 complete ====== + @ echo === layer 14 of 24 complete ====== @ echo ================================== layer15done: @ echo ================================== - @ echo === layer 15 of 23 complete ====== + @ echo === layer 15 of 24 complete ====== @ echo ================================== layer16done: @ echo ================================== - @ echo === layer 16 of 23 complete ====== + @ echo === layer 16 of 24 complete ====== @ echo ================================== layer17done: @ echo ================================== - @ echo === layer 17 of 23 complete ====== + @ echo === layer 17 of 24 complete ====== @ echo ================================== layer18done: @ echo ================================== - @ echo === layer 18 of 23 complete ====== + @ echo === layer 18 of 24 complete ====== @ echo ================================== layer19done: @ echo ================================== - @ echo === layer 19 of 23 complete ====== + @ echo === layer 19 of 24 complete ====== @ echo ================================== layer20done: @ echo ================================== - @ echo === layer 20 of 23 complete ====== + @ echo === layer 20 of 24 complete ====== @ echo ================================== layer21done: @ echo ================================== - @ echo === layer 21 of 23 complete ====== + @ echo === layer 21 of 24 complete ====== @ echo ================================== layer22done: @ echo ================================== - @ echo === layer 22 of 23 complete ====== + @ echo === layer 22 of 24 complete ====== @ echo ================================== layer23done: @ echo ================================== - @ echo === layer 23 of 23 complete ====== + @ echo === layer 23 of 24 complete ====== + @ echo ================================== + +layer24done: + @ echo ================================== + @ echo === layer 24 of 24 complete ====== @ echo ================================== @ @@ -2476,18 +2687,23 @@ ${HELP}/Library.help: ${BOOKS}/bookvol10.3.pamphlet ${HELP}/LieExponentials.help: ${BOOKS}/bookvol10.3.pamphlet @echo 7040 create LieExponentials.help from \ ${BOOKS}/bookvol10.3.pamphlet - @${TANGLE} -R"LieExponentials.help" ${BOOKS}/bookvol10.3.pamphlet \ + @${TANGLE} -R"LieExponentials.help" \ + ${BOOKS}/bookvol10.3.pamphlet \ >${HELP}/LieExponentials.help @cp ${HELP}/LieExponentials.help ${HELP}/LEXP.help - @${TANGLE} -R"LieExponentials.input" ${BOOKS}/bookvol10.3.pamphlet \ + @${TANGLE} -R"LieExponentials.input" \ + ${BOOKS}/bookvol10.3.pamphlet \ >${INPUT}/LieExponentials.input ${HELP}/LiePolynomial.help: ${BOOKS}/bookvol10.3.pamphlet - @echo 7041 create LiePolynomial.help from ${BOOKS}/bookvol10.3.pamphlet - @${TANGLE} -R"LiePolynomial.help" ${BOOKS}/bookvol10.3.pamphlet \ + @echo 7041 create LiePolynomial.help from \ + ${BOOKS}/bookvol10.3.pamphlet + @${TANGLE} -R"LiePolynomial.help" \ + ${BOOKS}/bookvol10.3.pamphlet \ >${HELP}/LiePolynomial.help @cp ${HELP}/LiePolynomial.help ${HELP}/LPOLY.help - @${TANGLE} -R"LiePolynomial.input" ${BOOKS}/bookvol10.3.pamphlet \ + @${TANGLE} -R"LiePolynomial.input" \ + ${BOOKS}/bookvol10.3.pamphlet \ >${INPUT}/LiePolynomial.input ${HELP}/LinearOrdinaryDifferentialOperator.help: ${BOOKS}/bookvol10.3.pamphlet @@ -2534,19 +2750,25 @@ ${HELP}/List.help: ${BOOKS}/bookvol10.3.pamphlet >${INPUT}/List.input ${HELP}/LyndonWord.help: ${BOOKS}/bookvol10.3.pamphlet - @echo 7046 create LyndonWord.help from ${BOOKS}/bookvol10.3.pamphlet - @${TANGLE} -R"LyndonWord.help" ${BOOKS}/bookvol10.3.pamphlet \ + @echo 7046 create LyndonWord.help from \ + ${BOOKS}/bookvol10.3.pamphlet + @${TANGLE} -R"LyndonWord.help" \ + ${BOOKS}/bookvol10.3.pamphlet \ >${HELP}/LyndonWord.help @cp ${HELP}/LyndonWord.help ${HELP}/LWORD.help - @${TANGLE} -R"LyndonWord.input" ${BOOKS}/bookvol10.3.pamphlet \ + @${TANGLE} -R"LyndonWord.input" \ + ${BOOKS}/bookvol10.3.pamphlet \ >${INPUT}/LyndonWord.input -${HELP}/Magma.help: ${BOOKS}/bookvol10.3.pamphlet - @echo 7047 create Magma.help from ${BOOKS}/bookvol10.3.pamphlet - @${TANGLE} -R"Magma.help" ${BOOKS}/bookvol10.3.pamphlet \ +${HELP}/Magma.help: ${BOOKS}/bookvol10.3.pamphlet + @echo 7047 create Magma.help from \ + ${BOOKS}/bookvol10.3.pamphlet + @${TANGLE} -R"Magma.help" \ + ${BOOKS}/bookvol10.3.pamphlet \ >${HELP}/Magma.help @-cp ${HELP}/Magma.help ${HELP}/MAGMA.help - @${TANGLE} -R"Magma.input" ${BOOKS}/bookvol10.3.pamphlet \ + @${TANGLE} -R"Magma.input" \ + ${BOOKS}/bookvol10.3.pamphlet \ >${INPUT}/Magma.input ${HELP}/MakeFunction.help: ${IN}/mkfunc.spad.pamphlet @@ -2954,20 +3176,25 @@ ${HELP}/XPBWPolynomial.help: ${BOOKS}/bookvol10.3.pamphlet >${INPUT}/XPBWPolynomial.input ${HELP}/XPolynomial.help: ${BOOKS}/bookvol10.3.pamphlet - @echo 7092 create XPolynomial.help from ${BOOKS}/bookvol10.3.pamphlet - @${TANGLE} -R"XPolynomial.help" ${BOOKS}/bookvol10.3.pamphlet \ + @echo 7092 create XPolynomial.help from \ + ${BOOKS}/bookvol10.3.pamphlet + @${TANGLE} -R"XPolynomial.help" \ + ${BOOKS}/bookvol10.3.pamphlet \ >${HELP}/XPolynomial.help @cp ${HELP}/XPolynomial.help ${HELP}/XPOLY.help - @${TANGLE} -R"XPolynomial.input" ${BOOKS}/bookvol10.3.pamphlet \ + @${TANGLE} -R"XPolynomial.input" \ + ${BOOKS}/bookvol10.3.pamphlet \ >${INPUT}/XPolynomial.input ${HELP}/XPolynomialRing.help: ${BOOKS}/bookvol10.3.pamphlet @echo 7093 create XPolynomialRing.help from \ ${BOOKS}/bookvol10.3.pamphlet - @${TANGLE} -R"XPolynomialRing.help" ${BOOKS}/bookvol10.3.pamphlet \ + @${TANGLE} -R"XPolynomialRing.help" \ + ${BOOKS}/bookvol10.3.pamphlet \ >${HELP}/XPolynomialRing.help @cp ${HELP}/XPolynomialRing.help ${HELP}/XPR.help - @${TANGLE} -R"XPolynomialRing.input" ${BOOKS}/bookvol10.3.pamphlet \ + @${TANGLE} -R"XPolynomialRing.input" \ + ${BOOKS}/bookvol10.3.pamphlet \ >${INPUT}/XPolynomialRing.input ${HELP}/ZeroDimensionalSolvePackage.help: ${IN}/zerodim.spad.pamphlet @@ -3022,6 +3249,8 @@ all: src ${OUT}/libdb.text ${DOCFILES} ${SPADBIN}/index.html \ ${SPADHELP} @ echo 4302 finished ${IN} +<> + ${SPADBIN}/index.html: @ echo 5000 making ${SPADBIN}/index.html @ echo "Axiom Algebra" >${SPADBIN}/index.html diff --git a/src/algebra/exposed.lsp.pamphlet b/src/algebra/exposed.lsp.pamphlet index 9701e7f..a577f11 100644 --- a/src/algebra/exposed.lsp.pamphlet +++ b/src/algebra/exposed.lsp.pamphlet @@ -160,6 +160,12 @@ (|GraphicsDefaults| . GRDEF) (|GroebnerPackage| . GB) (|GroebnerFactorizationPackage| . GBF) + (|Guess| . GUESS) + (|GuessOption| . GOPT) + (|GuessAlgebraicNumber| . GUESSAN) + (|GuessInteger| . GUESSINT) + (|GuessPolynomial| . GUESSP) + (|GuessFinite| . GUESSF) (|HallBasis| . HB) (|Heap| . HEAP) (|HexadecimalExpansion| . HEXADEC) diff --git a/src/algebra/fffg.spad.pamphlet b/src/algebra/fffg.spad.pamphlet index c49b7f4..8402fb0 100644 --- a/src/algebra/fffg.spad.pamphlet +++ b/src/algebra/fffg.spad.pamphlet @@ -1,5 +1,5 @@ \documentclass{article} -\usepackage{axiom,amsthm} +\usepackage{axiom,amsthm,amsmath,amssymb} \newtheorem{ToDo}{ToDo}[section] \begin{document} \title{fffg.spad} @@ -48,7 +48,7 @@ FiniteAbelianMonoidRingFunctions2(E: OrderedAbelianMonoid, ++ Applications 22. FractionFreeFastGaussian(D, V): Exports == Implementation where D: Join(IntegralDomain, GcdDomain) - V: AbelianMonoidRing(D, NonNegativeInteger) + V: AbelianMonoidRing(D, NonNegativeInteger) -- for example, SUP D SUP ==> SparseUnivariatePolynomial @@ -62,8 +62,18 @@ FractionFreeFastGaussian(D, V): Exports == Implementation where ++ \spad{fffg} is the general algorithm as proposed by Beckermann and ++ Labahn. ++ - ++ The second argument, c(k, M) computes c_k(M) as in Equation (2). Note - ++ that the information about f is therefore encoded in c. + + ++ The first argument is the list of c_{i,i}. These are the only values + ++ of C explicitely needed in \spad{fffg}. + ++ + ++ The second argument c, computes c_k(M), i.e., c_k(.) is the dual basis + ++ of the vector space V, but also knows about the special multiplication + ++ rule as descibed in Equation (2). Note that the information about f + ++ is therefore encoded in c. + ++ + ++ The third argument is the vector of degree bounds n, as introduced in + ++ Definition 2.1. In particular, the sum of the entries is the order of + ++ the Mahler system computed. interpolate: (List D, List D, NonNegativeInteger) -> Fraction SUP D ++ \spad{interpolate(xlist, ylist, deg} returns the rational function with @@ -77,11 +87,11 @@ FractionFreeFastGaussian(D, V): Exports == Implementation where \begin{ToDo} The following function could be moved to [[FFFGF]], parallel to - [[generalInterpolation]]. However, the reason for moving + [[generalInterpolation]]. However, the reason for moving [[generalInterpolation]] for fractions to a separate package was the need of - a generic signature, hence the extra argument [[VF]] to [[FFFGF]]. In the + a generic signature, hence the extra argument [[VF]] to [[FFFGF]]. In the special case of rational interpolation, this extra argument is not necessary, - since we are always returning a fraction of [[SUP]]s, and ignore [[V]]. In + since we are always returning a fraction of [[SUP]]s, and ignore [[V]]. In fact, [[V]] is not needed for [[fffg]] itself, only if we want to specify a [[CoeffAction]]. @@ -229,7 +239,7 @@ FractionFreeFastGaussian(D, V): Exports == Implementation where generalInterpolation(C: List D, coeffAction: CoeffAction, f: Vector V, eta: List NonNegativeInteger): Matrix SUP D == - + c: cFunction := generalCoefficient(coeffAction, f, (#1-1)::NonNegativeInteger, #2) fffg(C, c, eta) @@ -270,7 +280,7 @@ non-negative components smaller than [[p]] with the same sum as [[v]]. else v.pos := v.pos + 1 v.(pos-1) := (v.(pos-1) - 1)::NonNegativeInteger - + v @ @@ -401,83 +411,107 @@ interpolation problem, if [[M.1.(2,1)]] vanishes. ------------------------------------------------------------------------------- -- fffg ------------------------------------------------------------------------------- +@ + +[[recurrence]] computes the new matrix $M$, according to the following formulas +(cf. Table~2 in Beckermann and Labahn): +\begin{align*} + &\text{Increase order}\\ + &\quad\quad\text{for $\ell=1\dots m$, $\ell\neq\pi$}\\ + &\quad\quad\quad\quad\mathbf M_{\sigma+1}^{(.,\ell)} := + \left(\mathbf M_{\sigma}^{(.,\ell)}r^{(\pi)} + - \mathbf M_{\sigma}^{(.,\pi)}r^{(\ell)}\right)/d_\sigma\\ + &\text{Increase order in column $\pi$}\\ + &\quad\quad\mathbf M_{\sigma+1}^{(.,\pi)} := + \left(z-c_{\sigma,\sigma}\right)\mathbf M_{\sigma}^{(.,\pi)}\\ + &\text{Adjust degree constraints:}\\ + &\quad\quad\mathbf M_{\sigma+1}^{(.,\pi)} := + \left(\mathbf M_{\sigma+1}^{(.,\pi)}r^{(\pi)} + - \sum_{\ell\neq\pi}\mathbf M_{\sigma+1}^{(.,\ell)}p^{(\ell)} + \right)/d_\sigma +\end{align*} + +Since we do not need the matrix $\mathbf M_{\sigma}$ anymore, we drop the index +and update the matrix destructively. In the following, we write [[Ck]] for +$c_{\sigma,\sigma}$. +<>= -- a major part of the time is spent here - recurrence(M: Matrix SUP D, lambda: NonNegativeInteger, m: NonNegativeInteger, + recurrence(M: Matrix SUP D, pi: NonNegativeInteger, m: NonNegativeInteger, r: Vector D, d: D, z: SUP D, Ck: D, p: Vector D): Matrix SUP D == - rLambda := qelt(r, lambda) - polyf := rLambda * (z - Ck::SUP D) - - for i in 1..m repeat - Milambda := qelt(M, i, lambda) - newMilambda := polyf * Milambda - + rPi: D := qelt(r, pi) + polyf: SUP D := rPi * (z - Ck::SUP D) --- update columns ~= lambda and calculate their sum - for l in 1..m | l ~= lambda repeat - rl := qelt(r, l) - Mil := ((qelt(M, i, l) * rLambda - Milambda * rl) exquo d)::SUP D - qsetelt!(M, i, l, Mil) + for i in 1..m repeat + MiPi: SUP D := qelt(M, i, pi) + newMiPi: SUP D := polyf * MiPi - pl := qelt(p, l) - newMilambda := newMilambda - pl * Mil +-- update columns ~= pi and calculate their sum + for l in 1..m | l ~= pi repeat + rl: D := qelt(r, l) +-- I need the coercion to SUP D, since exquo returns an element of +-- Union("failed", SUP D)... + Mil: SUP D := ((qelt(M, i, l) * rPi - MiPi * rl) exquo d)::SUP D + qsetelt!(M, i, l, Mil) --- update column lambda + pl: D := qelt(p, l) + newMiPi := newMiPi - pl * Mil - qsetelt!(M, i, lambda, (newMilambda exquo d)::SUP D) +-- update column pi + qsetelt!(M, i, pi, (newMiPi exquo d)::SUP D) - M + M fffg(C: List D, c: cFunction, eta: List NonNegativeInteger): Matrix SUP D == - -- eta is the vector of degrees. We compute M with degrees eta+e_i-1, i=1..m - z: SUP D := monomial(1, 1) - m: NonNegativeInteger := #eta - M: Matrix SUP D := scalarMatrix(m, 1) - d: D := 1 - K: NonNegativeInteger := reduce(_+, eta) - etak: Vector NonNegativeInteger := zero(m) - r: Vector D := zero(m) - p: Vector D := zero(m) - Lambda: List Integer - lambdaMax: Integer - lambda: NonNegativeInteger - - for k in 1..K repeat + z: SUP D := monomial(1, 1) + m: NonNegativeInteger := #eta + M: Matrix SUP D := scalarMatrix(m, 1) + d: D := 1 + K: NonNegativeInteger := reduce(_+, eta) + etak: Vector NonNegativeInteger := zero(m) + r: Vector D := zero(m) + p: Vector D := zero(m) + Lambda: List Integer + lambdaMax: Integer + lambda: NonNegativeInteger + + for k in 1..K repeat -- k = sigma+1 - for l in 1..m repeat r.l := c(k, column(M, l)) + for l in 1..m repeat r.l := c(k, column(M, l)) - Lambda := [eta.l-etak.l for l in 1..m | r.l ~= 0] + Lambda := [eta.l-etak.l for l in 1..m | r.l ~= 0] -- if Lambda is empty, then M, d and etak remain unchanged. Otherwise, we look -- for the next closest para-normal point. - (empty? Lambda) => "iterate" + (empty? Lambda) => "iterate" - lambdaMax := reduce(max, Lambda) - lambda := 1 - while eta.lambda-etak.lambda < lambdaMax or r.lambda = 0 repeat - lambda := lambda + 1 + lambdaMax := reduce(max, Lambda) + lambda := 1 + while eta.lambda-etak.lambda < lambdaMax or r.lambda = 0 repeat + lambda := lambda + 1 -- Calculate leading coefficients - for l in 1..m | l ~= lambda repeat - if etak.l > 0 then - p.l := coefficient(M.(l, lambda), (etak.l-1)::NonNegativeInteger) - else - p.l := 0 + for l in 1..m | l ~= lambda repeat + if etak.l > 0 then + p.l := coefficient(M.(l, lambda), + (etak.l-1)::NonNegativeInteger) + else + p.l := 0 -- increase order and adjust degree constraints - M := recurrence(M, lambda, m, r, d, z, C.k, p) + M := recurrence(M, lambda, m, r, d, z, C.k, p) - d := r.lambda - etak.lambda := etak.lambda + 1 + d := r.lambda + etak.lambda := etak.lambda + 1 - M + M @ %$ @@ -540,7 +574,7 @@ FractionFreeFastGaussianFractions(D, V, VF): Exports == Implementation where n := #f g: Vector V := new(n, 0) den: Vector D := new(n, 0) - + for i in 1..n repeat c := coefficients(f.i) den.i := commonDenominator(c)$CommonDenominator(D, F, List F) @@ -548,7 +582,7 @@ FractionFreeFastGaussianFractions(D, V, VF): Exports == Implementation where $FAMR2(NonNegativeInteger, Fraction D, VF, D, V) M := generalInterpolation(C, coeffAction, g, eta)$FFFG(D, V) - + -- The following is necessary since I'm multiplying each row with a factor, not -- each column. Possibly I could factor out gcd den, but I'm not sure whether -- this is efficient. @@ -560,24 +594,24 @@ FractionFreeFastGaussianFractions(D, V, VF): Exports == Implementation where sumEta: NonNegativeInteger, maxEta: NonNegativeInteger) : Stream Matrix SUP D == - + n := #f g: Vector V := new(n, 0) den: Vector D := new(n, 0) - + for i in 1..n repeat c := coefficients(f.i) den.i := commonDenominator(c)$CommonDenominator(D, F, List F) g.i := map(retract(#1*den.i)@D, f.i) $FAMR2(NonNegativeInteger, Fraction D, VF, D, V) - + c: cFunction := generalCoefficient(coeffAction, g, (#1-1)::NonNegativeInteger, #2)$FFFG(D, V) MS: Stream Matrix SUP D := generalInterpolation(C, coeffAction, g, sumEta, maxEta)$FFFG(D, V) - + -- The following is necessary since I'm multiplying each row with a factor, not -- each column. Possibly I could factor out gcd den, but I'm not sure whether -- this is efficient. diff --git a/src/algebra/mantepse.spad.pamphlet b/src/algebra/mantepse.spad.pamphlet index f21b766..5addd74 100644 --- a/src/algebra/mantepse.spad.pamphlet +++ b/src/algebra/mantepse.spad.pamphlet @@ -17,15 +17,6 @@ \GFUN. \end{abstract} \tableofcontents -\section{domain UFPS UnivariateFormalPowerSeries} -<>= -)abbrev domain UFPS UnivariateFormalPowerSeries -UnivariateFormalPowerSeries(Coef: Ring) == - UnivariateTaylorSeries(Coef, 'x, 0$Coef) - -@ -%$ - \section{package UFPS1 UnivariateFormalPowerSeriesFunctions} <>= )abbrev package UFPS1 UnivariateFormalPowerSeriesFunctions @@ -47,148 +38,6 @@ UnivariateFormalPowerSeriesFunctions(Coef: Ring): Exports == Implementation @ %$ -\section{domain GOPT GuessOption} -<>= -)abbrev domain GOPT GuessOption -++ Author: Martin Rubey -++ Description: GuessOption is a domain whose elements are various options used -++ by \spadtype{Guess}. -GuessOption(): Exports == Implementation where - - Exports == SetCategory with - - maxDerivative: Integer -> % - ++ maxDerivative(d) specifies the maximum derivative in an algebraic - ++ differential equation. maxDerivative(-1) specifies that the maximum - ++ derivative can be arbitrary. This option is expressed in the form - ++ \spad{maxDerivative == d}. - - maxShift: Integer -> % - ++ maxShift(d) specifies the maximum shift in a recurrence - ++ equation. maxShift(-1) specifies that the maximum shift can be - ++ arbitrary. This option is expressed in the form \spad{maxShift == d}. - - maxPower: Integer -> % - ++ maxPower(d) specifies the maximum degree in an algebraic differential - ++ equation. For example, the degree of (f'')^3 f' is 4. maxPower(-1) - ++ specifies that the maximum exponent can be arbitrary. This option is - ++ expressed in the form \spad{maxPower == d}. - - homogeneous: Boolean -> % - ++ homogeneous(d) specifies whether we allow only homogeneous algebraic - ++ differential equations. This option is expressed in the form - ++ \spad{homogeneous == d}. - - maxLevel: Integer -> % - ++ maxLevel(d) specifies the maximum number of recursion levels operators - ++ guessProduct and guessSum will be applied. maxLevel(-1) specifies - ++ that all levels are tried. This option is expressed in the form - ++ \spad{maxLevel == d}. - - maxDegree: Integer -> % - ++ maxDegree(d) specifies the maximum degree of the coefficient - ++ polynomials in an algebraic differential equation or a recursion with - ++ polynomial coefficients. For rational functions with an exponential - ++ term, \spad{maxDegree} bounds the degree of the denominator - ++ polynomial. - ++ maxDegree(-1) specifies that the maximum - ++ degree can be arbitrary. This option is expressed in the form - ++ \spad{maxDegree == d}. - - allDegrees: Boolean -> % - ++ allDegrees(d) specifies whether all possibilities of the degree vector - ++ - taking into account maxDegree - should be tried. This is mainly - ++ interesting for rational interpolation. This option is expressed in - ++ the form \spad{allDegrees == d}. - - safety: NonNegativeInteger -> % - ++ safety(d) specifies the number of values reserved for testing any - ++ solutions found. This option is expressed in the form \spad{safety == - ++ d}. - - one: Boolean -> % - ++ one(d) specifies whether we are happy with one solution. This option - ++ is expressed in the form \spad{one == d}. - - debug: Boolean -> % - ++ debug(d) specifies whether we want additional output on the - ++ progress. This option is expressed in the form \spad{debug == d}. - - functionName: Symbol -> % - ++ functionName(d) specifies the name of the function given by the - ++ algebraic differential equation or recurrence. This option is - ++ expressed in the form \spad{functionName == d}. - - variableName: Symbol -> % - ++ variableName(d) specifies the variable used in by the algebraic - ++ differential equation. This option is expressed in the form - ++ \spad{variableName == d}. - - indexName: Symbol -> % - ++ indexName(d) specifies the index variable used for the formulas. This - ++ option is expressed in the form \spad{indexName == d}. - - displayAsGF: Boolean -> % - ++ displayAsGF(d) specifies whether the result is a generating function - ++ or a recurrence. This option should not be set by the user, but rather - ++ by the HP-specification. - - option : (List %, Symbol) -> Union(Any, "failed") - ++ option() is not to be used at the top level; - ++ option determines internally which drawing options are indicated in - ++ a draw command. - - option?: (List %, Symbol) -> Boolean - ++ option?() is not to be used at the top level; - ++ option? internally returns true for drawing options which are - ++ indicated in a draw command, or false for those which are not. - - checkOptions: List % -> Void - ++ checkOptions checks whether an option is given twice - - Implementation ==> add - import AnyFunctions1(Boolean) - import AnyFunctions1(Symbol) - import AnyFunctions1(Integer) - import AnyFunctions1(NonNegativeInteger) - - Rep := Record(keyword: Symbol, value: Any) - - maxLevel d == ["maxLevel"::Symbol, d::Any] - maxDerivative d == ["maxDerivative"::Symbol, d::Any] - maxShift d == maxDerivative d - maxDegree d == ["maxDegree"::Symbol, d::Any] - allDegrees d == ["allDegrees"::Symbol, d::Any] - maxPower d == ["maxPower"::Symbol, d::Any] - safety d == ["safety"::Symbol, d::Any] - homogeneous d == ["homogeneous"::Symbol, d::Any] - debug d == ["debug"::Symbol, d::Any] - one d == ["one"::Symbol, d::Any] - functionName d == ["functionName"::Symbol, d::Any] - variableName d == ["variableName"::Symbol, d::Any] - indexName d == ["indexName"::Symbol, d::Any] - displayAsGF d == ["displayAsGF"::Symbol, d::Any] - - coerce(x:%):OutputForm == x.keyword::OutputForm = x.value::OutputForm - x:% = y:% == x.keyword = y.keyword and x.value = y.value - - option?(l, s) == - for x in l repeat - x.keyword = s => return true - false - - option(l, s) == - for x in l repeat - x.keyword = s => return(x.value) - "failed" - - checkOptions l == - if not empty? l then - if find((first l).keyword = #1.keyword, rest l) case "failed" - then checkOptions rest l - else error "GuessOption: Option specified twice" - -@ \section{package GOPT0 GuessOptionFunctions0} <>= @@ -213,7 +62,7 @@ GuessOptionFunctions0(): Exports == Implementation where maxShift: LGOPT -> Integer ++ maxShift returns the specified maxShift or -1 as default. - + maxDegree: LGOPT -> Integer ++ maxDegree returns the specified maxDegree or -1 as default. @@ -267,7 +116,7 @@ GuessOptionFunctions0(): Exports == Implementation where retract(opt :: Any)$AnyFunctions1(Integer) maxShift l == maxDerivative l - + maxDegree l == if (opt := option(l, "maxDegree" :: Symbol)) case "failed" then -1 @@ -279,13 +128,13 @@ GuessOptionFunctions0(): Exports == Implementation where false else retract(opt :: Any)$AnyFunctions1(Boolean) - + maxPower l == if (opt := option(l, "maxPower" :: Symbol)) case "failed" then -1 else retract(opt :: Any)$AnyFunctions1(Integer) - + safety l == if (opt := option(l, "safety" :: Symbol)) case "failed" then 1$NonNegativeInteger @@ -345,17 +194,17 @@ GuessOptionFunctions0(): Exports == Implementation where ++ \spadtype{GuessPolynomial}, etc. Guess(F, S, EXPRR, R, retract, coerce): Exports == Implementation where F: Field -- zB.: FRAC POLY PF 5 --- in F wird interpoliert und gerechnet +-- in F we interpolate und check S: GcdDomain --- in guessExpRat möchte ich die Wurzeln von Polynomen über F bestimmen. Wenn F --- ein Quotientenkörper ist, dann kann ich den Nenner loswerden. --- F wäre dann also QFCAT S +-- in guessExpRat I would like to determine the roots of polynomials in F. When +-- F is a quotientfield, I can get rid of the denominator. In this case F is +-- roughly QFCAT S R: Join(OrderedSet, IntegralDomain) -- zB.: FRAC POLY INT --- die Ergebnisse werden aber in EXPRR dargestellt +-- results are given as elements of EXPRR -- EXPRR: Join(ExpressionSpace, IntegralDomain, EXPRR: Join(FunctionSpace Integer, IntegralDomain, RetractableTo R, RetractableTo Symbol, @@ -369,16 +218,16 @@ Guess(F, S, EXPRR, R, retract, coerce): Exports == Implementation where ground? : % -> Boolean -- zB.: EXPR INT --- EXPR FRAC POLY INT ist ja verboten. Deshalb kann ich nicht einfach EXPR R --- verwenden --- EXPRR existiert, falls irgendwann einmal EXPR PF 5 unterstützt wird. +-- EXPR FRAC POLY INT is forbidden. Thus i cannot just use EXPR R + +-- EXPRR exists, in case at some point there is support for EXPR PF 5. --- das folgende möchte ich eigentlich loswerden. - +-- the following I really would like to get rid of + retract: R -> F -- zB.: i+->i coerce: F -> EXPRR -- zB.: i+->i --- Achtung EXPRR ~= EXPR R +-- attention: EXPRR ~= EXPR R LGOPT ==> List GuessOption GOPT0 ==> GuessOptionFunctions0 @@ -445,6 +294,8 @@ Section~\ref{sec:Hermite-Pade}. <>= --@<> + -- EXT ==> (Integer, V, V) -> FPOLYS + -- EXTEXPR ==> (Symbol, F, F) -> EXPRR @ <>= @@ -455,168 +306,184 @@ Section~\ref{sec:Hermite-Pade}. ++ \spadfun{guessADE} to the successive differences and quotients of ++ the list. Default options as described in ++ \spadtype{GuessOptionFunctions0} are used. - + guess: (List F, LGOPT) -> GUESSRESULT ++ \spad{guess(l, options)} applies recursively \spadfun{guessRec} ++ and \spadfun{guessADE} to the successive differences and quotients ++ of the list. The given options are used. - + guess: (List F, List GUESSER, List Symbol) -> GUESSRESULT ++ \spad{guess(l, guessers, ops)} applies recursively the given ++ guessers to the successive differences if ops contains the symbol ++ guessSum and quotients if ops contains the symbol guessProduct to ++ the list. Default options as described in ++ \spadtype{GuessOptionFunctions0} are used. - + guess: (List F, List GUESSER, List Symbol, LGOPT) -> GUESSRESULT ++ \spad{guess(l, guessers, ops)} applies recursively the given ++ guessers to the successive differences if ops contains the symbol ++ \spad{guessSum} and quotients if ops contains the symbol ++ \spad{guessProduct} to the list. The given options are used. - + guessExpRat: List F -> GUESSRESULT ++ \spad{guessExpRat l} tries to find a function of the form ++ n+->(a+b n)^n r(n), where r(n) is a rational function, that fits ++ l. - + guessExpRat: (List F, LGOPT) -> GUESSRESULT ++ \spad{guessExpRat(l, options)} tries to find a function of the ++ form n+->(a+b n)^n r(n), where r(n) is a rational function, that ++ fits l. - + + guessBinRat: List F -> GUESSRESULT + ++ \spad{guessBinRat(l, options)} tries to find a function of the + ++ form n+->binomial(a+b n, n) r(n), where r(n) is a rational + ++ function, that fits l. + + guessBinRat: (List F, LGOPT) -> GUESSRESULT + ++ \spad{guessBinRat(l, options)} tries to find a function of the + ++ form n+->binomial(a+b n, n) r(n), where r(n) is a rational + ++ function, that fits l. + if F has RetractableTo Symbol and S has RetractableTo Symbol then guessExpRat: Symbol -> GUESSER ++ \spad{guessExpRat q} returns a guesser that tries to find a - ++ function of the form n+->(a+b q^n)^n r(q^n), where r(n) is a - ++ rational function, that fits l. + ++ function of the form n+->(a+b q^n)^n r(q^n), where r(q^n) is a + ++ q-rational function, that fits l. + + guessBinRat: Symbol -> GUESSER + ++ \spad{guessBinRat q} returns a guesser that tries to find a + ++ function of the form n+->qbinomial(a+b n, n) r(n), where r(q^n) is a + ++ q-rational function, that fits l. guessHP: (LGOPT -> HPSPEC) -> GUESSER ++ \spad{guessHP f} constructs an operation that applies Hermite-Pade ++ approximation to the series generated by the given function f. - + guessADE: List F -> GUESSRESULT ++ \spad{guessADE l} tries to find an algebraic differential equation ++ for a generating function whose first Taylor coefficients are ++ given by l, using the default options described in ++ \spadtype{GuessOptionFunctions0}. - + guessADE: (List F, LGOPT) -> GUESSRESULT ++ \spad{guessADE(l, options)} tries to find an algebraic ++ differential equation for a generating function whose first Taylor ++ coefficients are given by l, using the given options. - + guessAlg: List F -> GUESSRESULT ++ \spad{guessAlg l} tries to find an algebraic equation for a ++ generating function whose first Taylor coefficients are given by ++ l, using the default options described in ++ \spadtype{GuessOptionFunctions0}. It is equivalent to ++ \spadfun{guessADE}(l, maxDerivative == 0). - + guessAlg: (List F, LGOPT) -> GUESSRESULT ++ \spad{guessAlg(l, options)} tries to find an algebraic equation ++ for a generating function whose first Taylor coefficients are ++ given by l, using the given options. It is equivalent to ++ \spadfun{guessADE}(l, options) with \spad{maxDerivative == 0}. - + guessHolo: List F -> GUESSRESULT ++ \spad{guessHolo l} tries to find an ordinary linear differential ++ equation for a generating function whose first Taylor coefficients ++ are given by l, using the default options described in ++ \spadtype{GuessOptionFunctions0}. It is equivalent to ++ \spadfun{guessADE}\spad{(l, maxPower == 1)}. - + guessHolo: (List F, LGOPT) -> GUESSRESULT ++ \spad{guessHolo(l, options)} tries to find an ordinary linear ++ differential equation for a generating function whose first Taylor ++ coefficients are given by l, using the given options. It is ++ equivalent to \spadfun{guessADE}\spad{(l, options)} with ++ \spad{maxPower == 1}. - + guessPade: (List F, LGOPT) -> GUESSRESULT ++ \spad{guessPade(l, options)} tries to find a rational function ++ whose first Taylor coefficients are given by l, using the given ++ options. It is equivalent to \spadfun{guessADE}\spad{(l, ++ maxDerivative == 0, maxPower == 1, allDegrees == true)}. - + guessPade: List F -> GUESSRESULT ++ \spad{guessPade(l, options)} tries to find a rational function ++ whose first Taylor coefficients are given by l, using the default ++ options described in \spadtype{GuessOptionFunctions0}. It is ++ equivalent to \spadfun{guessADE}\spad{(l, options)} with ++ \spad{maxDerivative == 0, maxPower == 1, allDegrees == true}. - + guessRec: List F -> GUESSRESULT ++ \spad{guessRec l} tries to find an ordinary difference equation ++ whose first values are given by l, using the default options ++ described in \spadtype{GuessOptionFunctions0}. - + guessRec: (List F, LGOPT) -> GUESSRESULT ++ \spad{guessRec(l, options)} tries to find an ordinary difference ++ equation whose first values are given by l, using the given ++ options. - + guessPRec: (List F, LGOPT) -> GUESSRESULT ++ \spad{guessPRec(l, options)} tries to find a linear recurrence ++ with polynomial coefficients whose first values are given by l, ++ using the given options. It is equivalent to ++ \spadfun{guessRec}\spad{(l, options)} with \spad{maxPower == 1}. - + guessPRec: List F -> GUESSRESULT ++ \spad{guessPRec l} tries to find a linear recurrence with ++ polynomial coefficients whose first values are given by l, using ++ the default options described in ++ \spadtype{GuessOptionFunctions0}. It is equivalent to ++ \spadfun{guessRec}\spad{(l, maxPower == 1)}. - + guessRat: (List F, LGOPT) -> GUESSRESULT ++ \spad{guessRat(l, options)} tries to find a rational function ++ whose first values are given by l, using the given options. It is ++ equivalent to \spadfun{guessRec}\spad{(l, maxShift == 0, maxPower ++ == 1, allDegrees == true)}. - + guessRat: List F -> GUESSRESULT ++ \spad{guessRat l} tries to find a rational function whose first ++ values are given by l, using the default options described in ++ \spadtype{GuessOptionFunctions0}. It is equivalent to ++ \spadfun{guessRec}\spad{(l, maxShift == 0, maxPower == 1, ++ allDegrees == true)}. - + diffHP: LGOPT -> HPSPEC ++ \spad{diffHP options} returns a specification for Hermite-Pade ++ approximation with the differential operator - + shiftHP: LGOPT -> HPSPEC ++ \spad{shiftHP options} returns a specification for Hermite-Pade ++ approximation with the shift operator - + if F has RetractableTo Symbol and S has RetractableTo Symbol then + shiftHP: Symbol -> (LGOPT -> HPSPEC) ++ \spad{shiftHP options} returns a specification for ++ Hermite-Pade approximation with the $q$-shift operator - + diffHP: Symbol -> (LGOPT -> HPSPEC) ++ \spad{diffHP options} returns a specification for Hermite-Pade ++ approximation with the $q$-dilation operator - + guessRec: Symbol -> GUESSER ++ \spad{guessRec q} returns a guesser that finds an ordinary ++ q-difference equation whose first values are given by l, using ++ the given options. - + guessPRec: Symbol -> GUESSER ++ \spad{guessPRec q} returns a guesser that tries to find ++ a linear q-recurrence with polynomial coefficients whose first ++ values are given by l, using the given options. It is ++ equivalent to \spadfun{guessRec}\spad{(q)} with ++ \spad{maxPower == 1}. - + guessRat: Symbol -> GUESSER ++ \spad{guessRat q} returns a guesser that tries to find a ++ q-rational function whose first values are given by l, using ++ the given options. It is equivalent to \spadfun{guessRec} with ++ \spad{(l, maxShift == 0, maxPower == 1, allDegrees == true)}. - + guessADE: Symbol -> GUESSER ++ \spad{guessADE q} returns a guesser that tries to find an ++ algebraic differential equation for a generating function whose @@ -627,22 +494,32 @@ Section~\ref{sec:Hermite-Pade}. @ <>= -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 +--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 + @ <>= @@ -659,6 +536,7 @@ testInterpolant: (List SUP S, List F, List UFPSSUPF, List EXPRR, List EXPRR, _ <> <> +<> <> <> @ @@ -686,6 +564,7 @@ SUPS2SUPF(p: SUP S): SUP F == to Fraction S" @ +%$ \subsection{guessing rational functions with an exponential term} @@ -740,7 +619,7 @@ SUPPOLYS2SUPF(p: SUP POLYS, a1v: F, Av: F): SUP F == 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) @@ -754,9 +633,19 @@ POLYF2EXPRR(p: POLYF): EXPRR == $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) + +-- does not work: +-- 6 +-- WARNING (genufact): No known algorithm to factor ? , trying square-free. + +-- GF ==> GenUFactorize F @ -\paragraph{mimicking $q$-anaologa} +\paragraph{mimicking $q$-analoga} <>= defaultD: DIFFSPECN @@ -772,8 +661,8 @@ DN2DL(DN, i) == retract(retract(DN(i::EXPRR))@R) \subsubsection{reducing the degree of the polynomials} -The degree of [[poly3]] is governed by $(a0+x_m a_1)^{x_m}$. Therefore, we -substitute $A-x_m a1$ for $a$, which reduces the degree in $a_1$ by +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}$. <>= @@ -896,10 +785,6 @@ still obtain a polynomial. \subsubsection{The main routine} <>= -GF ==> GeneralizedMultivariateFactorize(SingletonAsOrderedSet, - IndexedExponents V, F, F, - SUP F) - guessExpRatAux(xx: Symbol, list: List F, basis: DIFFSPECN, xValues: List Integer, options: LGOPT): List EXPRR == @@ -909,7 +794,7 @@ guessExpRatAux(xx: Symbol, list: List F, basis: DIFFSPECN, 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)) @@ -1034,7 +919,8 @@ a rational function. To obtain $y$, we compute $y(n)=s_n*(a+b n)^{-n}$. systemCommand("sys date +%s")$MoreSystemCommands output("computing gcd..."::OutputForm)$OutputPackage --- we want to solve over F +-- 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 @@ -1153,7 +1039,7 @@ guessExpRatAux0(list: List F, basis: DIFFSPECN, options: LGOPT): GUESSRESULT == 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)) @@ -1172,6 +1058,334 @@ if F has RetractableTo Symbol and S has RetractableTo Symbol then @ +\subsection{guessing rational functions with a binomial term} + +\begin{ToDo} + 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{ToDo} + + +<>= + +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)) + +@ + +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$. + +<>= + 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 + +@ + +\begin{ToDo} + Shouldn't we use the algorithm over [[POLYS]] here? Strange enough, for + polynomial interpolation, it is faster, but for rational interpolation + \emph{much} slower. This should be investigated. +\end{ToDo} + +\begin{ToDo} + It seems that [[maxDeg]] bounds the degree of the denominator, rather than + the numerator? This is now added to the documentation of [[maxDegree]], it + should make sense. +\end{ToDo} + +<>= +-- if debug(options)$GOPT0 then +-- 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)_ + +@ + + \subsection{Hermite Pad\'e interpolation}\label{sec:Hermite-Pade} <>= @@ -1193,12 +1407,12 @@ DIFFSPECX ==> (EXPRR, Symbol, NonNegativeInteger) -> EXPRR -- 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 @@ -1209,8 +1423,8 @@ DIFFSPEC1X ==> Symbol -> EXPRR \subsubsection{Streams}\label{sec:streams} -In this section we define some functions that provide streams for [[HermitePade]]. - +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 @@ -1317,7 +1531,7 @@ does not. Note that we conjugate all partitions at the very end of the following procedure\dots - + <>= FilteredPartitionStream(options: LGOPT): Stream List Integer == maxD := 1+maxDerivative(options)$GOPT0 @@ -1340,7 +1554,7 @@ FilteredPartitionStream(options: LGOPT): Stream List Integer == 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 @@ -1364,7 +1578,7 @@ applying one of the power and shift or differentiation operators. More precisely, the $n$\textsuperscript{th} entry of the stream takes into account all partitions up to index $n$. Thus, the entries of the stream are weakly increasing. - + <>= ADEdegreeStream(partitions: Stream List Integer): Stream NNI == scan(0, max((if empty? #1 then 0 else (first #1 - 1)::NNI), #2), @@ -1389,7 +1603,7 @@ ADEEXPRRStream(f: EXPRR, xx: Symbol, partitions: Stream List Integer, \subsubsection{Operators} We need to define operators that transform series for differentiation and -shifting. We also provide operators for $q$-analoga. The functionality +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}. @@ -1571,7 +1785,7 @@ shiftHP options == @ \paragraph{$q$-Shifting} The $q$-shift also transforms $u(n)$ into $u(n+1)$, -and we can reuse the corrseponding functions of the previous paragraph. +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. @@ -1716,7 +1930,7 @@ guessInterpolate2(guessList: List SUP F, 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" @ @@ -1809,8 +2023,14 @@ 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. + <>= - g: SUP S := gcd resi + 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) @@ -1877,12 +2097,10 @@ important if we apply, for example, [[guessRat]] recursively. In this case, x := variableName(options)$GOPT0 dummy := new$Symbol - initials: List EXPRR - if displayAsGF(options)$GOPT0 - then initials := [coerce(e)@EXPRR*factorial(k::EXPRR) _ - for e in list for k in 0..] - else initials := [coerce(e)@EXPRR for e in list] + initials: List EXPRR := [coerce(e)@EXPRR for e in list] + @ +%$ We need to create several streams. Let $P$ be the univariate power series whose first few elements are given by [[list]]. As an example, consider the @@ -1942,7 +2160,7 @@ which all the generated series are correct, taking into account [[safety]]. output(hconcat("guessDegree is ", guessDegree::OutputForm)) $OutputPackage @ -%$ +%$" We now have to distinguish between the case where we try all combination of degrees and the case where we try only an (nearly) evenly distributed vector of @@ -2004,10 +2222,19 @@ increase the order anymore. .quotient = 0) and (newMaxParams.remainder < o))) then iterate? := true - else eta: List NNI - := [(if i <= maxParams.remainder _ - then maxParams.quotient + 1 _ - else maxParams.quotient)::NNI for i in 1..o] + 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] @ We distribute the parameters as evenly as possible. Is it better to have @@ -2019,7 +2246,10 @@ HermitePade. \begin{ToDo} [[maxDegree]] should be handled differently, maybe: we should only pass as - many coefficients to [[FFFG]] as [[maxDegree]] implies! + 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{ToDo} <>= @@ -2035,7 +2265,7 @@ HermitePade. 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, @@ -2043,7 +2273,7 @@ HermitePade. repeat (empty? MS) => break M := first MS - + for i in 1..o repeat res := testInterpolant(entries column(M, i), list, @@ -2057,13 +2287,13 @@ HermitePade. 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, @@ -2073,10 +2303,10 @@ HermitePade. 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 @@ -2168,7 +2398,7 @@ 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) @@ -2219,7 +2449,7 @@ guess(list: List F, guessers: List GUESSER, ops: List Symbol, 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 @@ -2270,7 +2500,7 @@ guess(list: List F, guessers: List GUESSER, ops: List Symbol, for guess in sumGuess | not any?(guess.function = #1.function, res) repeat res := cons(guess, res) - + res guess(list: List F): GUESSRESULT == @@ -2308,15 +2538,15 @@ GuessInteger() == Guess(Fraction Integer, Integer, Expression Integer, ++ This package exports guessing of sequences of numbers in a finite field GuessFiniteFunctions(F:Join(FiniteFieldCategory, ConvertibleTo Integer)): Exports == Implementation where - + EXPRR ==> Expression Integer - + Exports == with F2EXPRR: F -> EXPRR - + Implementation == add - + F2EXPRR(p: F): EXPRR == convert(p)@Integer::EXPRR @ @@ -2345,6 +2575,320 @@ GuessPolynomial() == Guess(Fraction Polynomial Integer, Polynomial Integer, @ +\section{package GUESSAN GuessAlgebraicNumber} +<>= +)abbrev package GUESSAN GuessAlgebraicNumber +++ Description: +++ This package exports guessing of sequences of rational functions +GuessAlgebraicNumber() == Guess(AlgebraicNumber, AlgebraicNumber, + Expression Integer, + AlgebraicNumber, + id$MappingPackage1(AlgebraicNumber), + coerce$Expression(Integer)) + +@ + + +\section{package GUESSUP GuessUnivariatePolynomial} +<>= +)abbrev domain MYUP MyUnivariatePolynomial +MyUnivariatePolynomial(x:Symbol, R:Ring): + UnivariatePolynomialCategory(R) with + RetractableTo Symbol; + coerce: Variable(x) -> % + ++ coerce(x) converts the variable x to a univariate polynomial. + fmecg: (%,NonNegativeInteger,R,%) -> % + ++ fmecg(p1,e,r,p2) finds X : p1 - r * X**e * p2 + if R has univariate: (R, Symbol) -> SparseUnivariatePolynomial R + then coerce: R -> % + coerce: Polynomial R -> % + == SparseUnivariatePolynomial(R) add + Rep := SparseUnivariatePolynomial(R) + coerce(p: %):OutputForm == outputForm(p, outputForm x) + coerce(x: Symbol): % == monomial(1, 1) + coerce(v: Variable(x)):% == monomial(1, 1) + retract(p: %): Symbol == + retract(p)@SingletonAsOrderedSet + x + if R has univariate: (R, Symbol) -> SparseUnivariatePolynomial R + then coerce(p: R): % == univariate(p, x)$R + + coerce(p: Polynomial R): % == + poly: SparseUnivariatePolynomial(Polynomial R) + := univariate(p, x)$(Polynomial R) + map(retract(#1), poly)$UnivariatePolynomialCategoryFunctions2(Polynomial R, + SparseUnivariatePolynomial Polynomial R, + R, %) + + +)abbrev domain MYEXPR MyExpression +MyExpression(q: Symbol, R): Exports == Implementation where + + R: Join(Ring, OrderedSet, IntegralDomain) + UP ==> MyUnivariatePolynomial(q, R) + + Exports == Join(FunctionSpace R, IntegralDomain, + RetractableTo UP, RetractableTo Symbol, + RetractableTo Integer, CombinatorialOpsCategory, + PartialDifferentialRing Symbol) with + _* : (%,%) -> % + _/ : (%,%) -> % + _*_* : (%,%) -> % + numerator : % -> % + denominator : % -> % + ground? : % -> Boolean + + coerce: Fraction UP -> % + retract: % -> Fraction UP + + Implementation == Expression R add + Rep := Expression R + + iunivariate(p: Polynomial R): UP == + poly: SparseUnivariatePolynomial(Polynomial R) + := univariate(p, q)$(Polynomial R) + map(retract(#1), poly)$UnivariatePolynomialCategoryFunctions2(Polynomial R, + SparseUnivariatePolynomial Polynomial R, + R, UP) + + retract(p: %): Fraction UP == + poly: Fraction Polynomial R := retract p + upoly: UP := iunivariate numer poly + vpoly: UP := iunivariate denom poly + + upoly / vpoly + + retract(p: %): UP == iunivariate retract p + + coerce(r: Fraction UP): % == + num: SparseUnivariatePolynomial R := makeSUP numer r + den: SparseUnivariatePolynomial R := makeSUP denom r + u: Polynomial R := multivariate(num, q) + v: Polynomial R := multivariate(den, q) + + quot: Fraction Polynomial R := u/v + + quot::(Expression R) + + + +)abbrev package GUESSUP GuessUnivariatePolynomial +++ Description: +++ This package exports guessing of sequences of univariate rational functions +GuessUnivariatePolynomial(q: Symbol): Exports == Implementation where + + UP ==> MyUnivariatePolynomial(q, Integer) + EXPRR ==> MyExpression(q, Integer) + F ==> Fraction UP + S ==> UP + NNI ==> NonNegativeInteger + LGOPT ==> List GuessOption + GUESSRESULT ==> List Record(function: EXPRR, order: NNI) + GUESSER ==> (List F, LGOPT) -> GUESSRESULT + + Exports == with + + guess: List F -> GUESSRESULT + ++ \spad{guess l} applies recursively \spadfun{guessRec} and + ++ \spadfun{guessADE} to the successive differences and quotients of + ++ the list. Default options as described in + ++ \spadtype{GuessOptionFunctions0} are used. + + guess: (List F, LGOPT) -> GUESSRESULT + ++ \spad{guess(l, options)} applies recursively \spadfun{guessRec} + ++ and \spadfun{guessADE} to the successive differences and quotients + ++ of the list. The given options are used. + + guess: (List F, List GUESSER, List Symbol) -> GUESSRESULT + ++ \spad{guess(l, guessers, ops)} applies recursively the given + ++ guessers to the successive differences if ops contains the symbol + ++ guessSum and quotients if ops contains the symbol guessProduct to + ++ the list. Default options as described in + ++ \spadtype{GuessOptionFunctions0} are used. + + guess: (List F, List GUESSER, List Symbol, LGOPT) -> GUESSRESULT + ++ \spad{guess(l, guessers, ops)} applies recursively the given + ++ guessers to the successive differences if ops contains the symbol + ++ \spad{guessSum} and quotients if ops contains the symbol + ++ \spad{guessProduct} to the list. The given options are used. + + guessExpRat: List F -> GUESSRESULT + ++ \spad{guessExpRat l} tries to find a function of the form + ++ n+->(a+b n)^n r(n), where r(n) is a rational function, that fits + ++ l. + + guessExpRat: (List F, LGOPT) -> GUESSRESULT + ++ \spad{guessExpRat(l, options)} tries to find a function of the + ++ form n+->(a+b n)^n r(n), where r(n) is a rational function, that + ++ fits l. + + guessBinRat: List F -> GUESSRESULT + ++ \spad{guessBinRat(l, options)} tries to find a function of the + ++ form n+->binomial(a+b n, n) r(n), where r(n) is a rational + ++ function, that fits l. + + guessBinRat: (List F, LGOPT) -> GUESSRESULT + ++ \spad{guessBinRat(l, options)} tries to find a function of the + ++ form n+->binomial(a+b n, n) r(n), where r(n) is a rational + ++ function, that fits l. + +-- if F has RetractableTo Symbol and S has RetractableTo Symbol then + + guessExpRat: Symbol -> GUESSER + ++ \spad{guessExpRat q} returns a guesser that tries to find a + ++ function of the form n+->(a+b q^n)^n r(q^n), where r(q^n) is a + ++ q-rational function, that fits l. + + guessBinRat: Symbol -> GUESSER + ++ \spad{guessBinRat q} returns a guesser that tries to find a + ++ function of the form n+->qbinomial(a+b n, n) r(n), where r(q^n) is a + ++ q-rational function, that fits l. + +------------------------------------------------------------------------------- + + guessHP: (LGOPT -> HPSPEC) -> GUESSER + ++ \spad{guessHP f} constructs an operation that applies Hermite-Pade + ++ approximation to the series generated by the given function f. + + guessADE: List F -> GUESSRESULT + ++ \spad{guessADE l} tries to find an algebraic differential equation + ++ for a generating function whose first Taylor coefficients are + ++ given by l, using the default options described in + ++ \spadtype{GuessOptionFunctions0}. + + guessADE: (List F, LGOPT) -> GUESSRESULT + ++ \spad{guessADE(l, options)} tries to find an algebraic + ++ differential equation for a generating function whose first Taylor + ++ coefficients are given by l, using the given options. + + guessAlg: List F -> GUESSRESULT + ++ \spad{guessAlg l} tries to find an algebraic equation for a + ++ generating function whose first Taylor coefficients are given by + ++ l, using the default options described in + ++ \spadtype{GuessOptionFunctions0}. It is equivalent to + ++ \spadfun{guessADE}(l, maxDerivative == 0). + + guessAlg: (List F, LGOPT) -> GUESSRESULT + ++ \spad{guessAlg(l, options)} tries to find an algebraic equation + ++ for a generating function whose first Taylor coefficients are + ++ given by l, using the given options. It is equivalent to + ++ \spadfun{guessADE}(l, options) with \spad{maxDerivative == 0}. + + guessHolo: List F -> GUESSRESULT + ++ \spad{guessHolo l} tries to find an ordinary linear differential + ++ equation for a generating function whose first Taylor coefficients + ++ are given by l, using the default options described in + ++ \spadtype{GuessOptionFunctions0}. It is equivalent to + ++ \spadfun{guessADE}\spad{(l, maxPower == 1)}. + + guessHolo: (List F, LGOPT) -> GUESSRESULT + ++ \spad{guessHolo(l, options)} tries to find an ordinary linear + ++ differential equation for a generating function whose first Taylor + ++ coefficients are given by l, using the given options. It is + ++ equivalent to \spadfun{guessADE}\spad{(l, options)} with + ++ \spad{maxPower == 1}. + + guessPade: (List F, LGOPT) -> GUESSRESULT + ++ \spad{guessPade(l, options)} tries to find a rational function + ++ whose first Taylor coefficients are given by l, using the given + ++ options. It is equivalent to \spadfun{guessADE}\spad{(l, + ++ maxDerivative == 0, maxPower == 1, allDegrees == true)}. + + guessPade: List F -> GUESSRESULT + ++ \spad{guessPade(l, options)} tries to find a rational function + ++ whose first Taylor coefficients are given by l, using the default + ++ options described in \spadtype{GuessOptionFunctions0}. It is + ++ equivalent to \spadfun{guessADE}\spad{(l, options)} with + ++ \spad{maxDerivative == 0, maxPower == 1, allDegrees == true}. + + guessRec: List F -> GUESSRESULT + ++ \spad{guessRec l} tries to find an ordinary difference equation + ++ whose first values are given by l, using the default options + ++ described in \spadtype{GuessOptionFunctions0}. + + guessRec: (List F, LGOPT) -> GUESSRESULT + ++ \spad{guessRec(l, options)} tries to find an ordinary difference + ++ equation whose first values are given by l, using the given + ++ options. + + guessPRec: (List F, LGOPT) -> GUESSRESULT + ++ \spad{guessPRec(l, options)} tries to find a linear recurrence + ++ with polynomial coefficients whose first values are given by l, + ++ using the given options. It is equivalent to + ++ \spadfun{guessRec}\spad{(l, options)} with \spad{maxPower == 1}. + + guessPRec: List F -> GUESSRESULT + ++ \spad{guessPRec l} tries to find a linear recurrence with + ++ polynomial coefficients whose first values are given by l, using + ++ the default options described in + ++ \spadtype{GuessOptionFunctions0}. It is equivalent to + ++ \spadfun{guessRec}\spad{(l, maxPower == 1)}. + + guessRat: (List F, LGOPT) -> GUESSRESULT + ++ \spad{guessRat(l, options)} tries to find a rational function + ++ whose first values are given by l, using the given options. It is + ++ equivalent to \spadfun{guessRec}\spad{(l, maxShift == 0, maxPower + ++ == 1, allDegrees == true)}. + + guessRat: List F -> GUESSRESULT + ++ \spad{guessRat l} tries to find a rational function whose first + ++ values are given by l, using the default options described in + ++ \spadtype{GuessOptionFunctions0}. It is equivalent to + ++ \spadfun{guessRec}\spad{(l, maxShift == 0, maxPower == 1, + ++ allDegrees == true)}. + + diffHP: LGOPT -> HPSPEC + ++ \spad{diffHP options} returns a specification for Hermite-Pade + ++ approximation with the differential operator + + shiftHP: LGOPT -> HPSPEC + ++ \spad{shiftHP options} returns a specification for Hermite-Pade + ++ approximation with the shift operator + +-- if F has RetractableTo Symbol and S has RetractableTo Symbol then + + shiftHP: Symbol -> (LGOPT -> HPSPEC) + ++ \spad{shiftHP options} returns a specification for + ++ Hermite-Pade approximation with the $q$-shift operator + + diffHP: Symbol -> (LGOPT -> HPSPEC) + ++ \spad{diffHP options} returns a specification for Hermite-Pade + ++ approximation with the $q$-dilation operator + + guessRec: Symbol -> GUESSER + ++ \spad{guessRec q} returns a guesser that finds an ordinary + ++ q-difference equation whose first values are given by l, using + ++ the given options. + + guessPRec: Symbol -> GUESSER + ++ \spad{guessPRec q} returns a guesser that tries to find + ++ a linear q-recurrence with polynomial coefficients whose first + ++ values are given by l, using the given options. It is + ++ equivalent to \spadfun{guessRec}\spad{(q)} with + ++ \spad{maxPower == 1}. + + guessRat: Symbol -> GUESSER + ++ \spad{guessRat q} returns a guesser that tries to find a + ++ q-rational function whose first values are given by l, using + ++ the given options. It is equivalent to \spadfun{guessRec} with + ++ \spad{(l, maxShift == 0, maxPower == 1, allDegrees == true)}. + + guessADE: Symbol -> GUESSER + ++ \spad{guessADE q} returns a guesser that tries to find an + ++ algebraic differential equation for a generating function whose + ++ first Taylor coefficients are given by l, using the given + ++ options. + +------------------------------------------------------------------------------- + + Implementation == Guess(Fraction UP, UP, + MyExpression(q, Integer), + Fraction UP, + id$MappingPackage1(Fraction UP), + coerce$MyExpression(q, Integer)) +@ + + \section{License} <>= --Copyright (c) 2006-2007, Martin Rubey @@ -2377,14 +2921,14 @@ GuessPolynomial() == Guess(Fraction Polynomial Integer, Polynomial Integer, <<*>>= <> -<> <> -<> <> <> <> +<> <> <> <> +<> @ \end{document} diff --git a/src/algebra/rec.spad.pamphlet b/src/algebra/rec.spad.pamphlet index 5be79c2..8cedd33 100644 --- a/src/algebra/rec.spad.pamphlet +++ b/src/algebra/rec.spad.pamphlet @@ -21,7 +21,8 @@ RecurrenceOperator(R, F): Exports == Implementation where R: Join(OrderedSet, IntegralDomain, ConvertibleTo InputForm) F: Join(FunctionSpace R, AbelianMonoid, RetractableTo Integer, _ - RetractableTo Symbol, PartialDifferentialRing Symbol) + RetractableTo Symbol, PartialDifferentialRing Symbol, _ + CombinatorialOpsCategory) --RecurrenceOperator(F): Exports == Implementation where -- F: Join(ExpressionSpace, AbelianMonoid, RetractableTo Integer, -- RetractableTo Symbol, PartialDifferentialRing Symbol) @@ -39,14 +40,20 @@ RecurrenceOperator(R, F): Exports == Implementation where evalADE: (BasicOperator, Symbol, F, F, F, List F) -> F ++ \spad{evalADE(f, dummy, x, n, eq, values)} creates an expression that - ++ stands for the coefficient of x^n in f(x), where f(x) is given by the - ++ functional equation eq. However, for technical reasons the variable x - ++ has to be replaced by a dummy variable dummy in eq. The argument - ++ values specifies f(0), f'(0), f''(0),... + ++ stands for the coefficient of x^n in the Taylor expansion of f(x), + ++ where f(x) is given by the functional equation eq. However, for + ++ technical reasons the variable x has to be replaced by a dummy + ++ variable dummy in eq. The argument values specifies the first few + ++ Taylor coefficients. + + getEq: F -> F + ++ \spad{getEq f} returns the defining equation, if f represents the + ++ coefficient of an ADE or a recurrence. + + getOp: F -> BasicOperator + ++ \spad{getOp f}, if f represents the coefficient of a recurrence or + ++ ADE, returns the operator representing the solution - getADE: F -> F - ++ \spad{getADE f} returns the defining equation, if f stands for the - ++ coefficient of an ADE. -- should be local numberOfValuesNeeded: (Integer, BasicOperator, Symbol, F) -> Integer @@ -165,7 +172,7 @@ how many initial values are necessary. p := univariate(a.1, retract(n::F)@Kernel(F)) if denominator p ~= 1 then return "failed" - + num := numer p if degree num = 1 and coefficient(num, 1) = 1 @@ -249,14 +256,14 @@ how many initial values are necessary. then values.(len-n) else shiftInfo := shiftInfoRec(op, argsym, eq) - + if shiftInfo.max case Integer then p := univariate(eq, shiftInfo.ker) num := numer p - + if degree num = 1 then - + next := -coefficient(num, 0)/coefficient(num, 1) nextval := eval(next, argsym::F, (len-(shiftInfo.max)::Integer)::F) @@ -266,7 +273,7 @@ how many initial values are necessary. else kernel(oprecur, append([eq, argsym::F, argdisp, op(arg)], values)) - + else kernel(oprecur, append([eq, argsym::F, argdisp, op(arg)], values)) @@ -334,6 +341,24 @@ an operator. \subsection{Functional Equations} +\subsubsection{Determining the number of initial values for ADE's} + +We use Joris van der Hoeven's instructions for ADE's. Given +$Q=p(f,f',\dots,f^{(r)})$ we first need to differentiate $Q$ with respect to +$f^{(i)}$ for $i\in\{0,1,\dots,r\}$, plug in the given truncated power series +solution and determine the valuation. + +<>= + getValuation(op, argsym, eq, maxorder, values): Integer == + max: Integer := -1; + ker: Kernel F + for i in 0..maxorder repeat + ker := D(op(argsym), argsym, i)::Kernel F + pol := univariate(eq, ker) + dif := D pol + ground numer(dif.D(op(argsym), argsym, i)) +@ + \subsubsection{Extracting some information from the functional equation} [[getOrder]] returns the maximum derivative of [[op]] occurring in [[f]]. @@ -360,51 +385,55 @@ an operator. if ((n := retractIfCan(arg)@Union(Integer, "failed")) case "failed") then --- --- The following would need yet another operator, namely "coefficient of". --- --- keq := kernels eq --- order := getOrder(op, keq.1) --- for k in rest keq repeat order := max(order, getOrder(op, k)) --- --- if zero? order then --- -- in this case, we have an algebraic equation --- --- p: Fraction SparseUnivariatePolynomial F --- := univariate(eq, kernels(D(op(argsym::F), argsym, order)).1)$F --- --- num := numer p --- --- if one? degree num then --- -- the equation is in fact linear --- return eval(-coefficient(num, 0)/coefficient(num, 1), argsym::F, arg::F) --- - kernel(opADE, - append([eq, argsym::F, argdisp, op(arg)], values)) +-- try to determine the necessary number of initial values + keq := kernels eq + order := getOrder(op, keq.1) + for k in rest keq repeat order := max(order, getOrder(op, k)) + + p: Fraction SparseUnivariatePolynomial F + := univariate(eq, kernels(D(op(argsym::F), argsym, order)).1)$F + + if one? degree numer p +-- the equation is holonomic + then kernel(opADE, + append([eq, argsym::F, argdisp, op(arg)], + reverse first(reverse values, order))) + else kernel(opADE, + append([eq, argsym::F, argdisp, op(arg)], values)) else if n < 0 then 0 else keq := kernels eq order := getOrder(op, keq.1) - output(hconcat("The order is ", order::OutputForm))$OutputPackage +-- output(hconcat("The order is ", order::OutputForm))$OutputPackage for k in rest keq repeat order := max(order, getOrder(op, k)) p: Fraction SparseUnivariatePolynomial F := univariate(eq, kernels(D(op(argsym::F), argsym, order)).1)$F - output(hconcat("p: ", p::OutputForm))$OutputPackage +-- output(hconcat("p: ", p::OutputForm))$OutputPackage + if degree numer p > 1 then - kernel(opADE, - append([eq, argsym::F, argdisp, op(arg)], values)) +-- kernel(opADE, +-- append([eq, argsym::F, argdisp, op(arg)], values)) + + s := seriesSolve(eq, op, argsym, reverse values) + $ExpressionSolve(R, F, + UnivariateFormalPowerSeries F, + UnivariateFormalPowerSeries + SparseUnivariatePolynomialExpressions F) + + elt(s, n::Integer::NonNegativeInteger) + else s := seriesSolve(eq, op, argsym, first(reverse values, order)) $ExpressionSolve(R, F, UnivariateFormalPowerSeries F, UnivariateFormalPowerSeries SparseUnivariatePolynomialExpressions F) - + elt(s, n::Integer::NonNegativeInteger) @@ -418,13 +447,25 @@ an operator. evaluate(opADE, iADE)$BasicOperatorFunctions1(F) - getADE(f: F): F == + getEq(f: F): F == ker := kernels f - if one?(#ker) and is?(operator(ker.1), "rootOfADE"::Symbol) then + if one?(#ker) and _ + (is?(operator(ker.1), "rootOfADE"::Symbol) or _ + is?(operator(ker.1), "rootOfRec"::Symbol)) then l := argument(ker.1) eval(eqAsF l, dummyAsF l, displayVariable l) else - error "getADE: argument should be a single rootOfADE object" + error "getEq: argument should be a single rootOfADE or rootOfRec object" + + getOp(f: F): BasicOperator == + ker := kernels f + if one?(#ker) and _ + (is?(operator(ker.1), "rootOfADE"::Symbol) or _ + is?(operator(ker.1), "rootOfRec"::Symbol)) then + operatorName argument(ker.1) + else + error "getOp: argument should be a single rootOfADE or rootOfRec object" + @ %$ @@ -440,7 +481,9 @@ an operator. := cons(eval(eqAsF l, dummyAsF l, displayVariable l)::OutputForm = _ 0::OutputForm, [eval(D(op(dummyAsF l), dummy l, i), _ - dummyAsF l=0)::OutputForm = (values.(i+1))::OutputForm + dummyAsF l=0)::OutputForm = _ + (values.(i+1))::OutputForm * _ + factorial(box(i::R::F)$F)::OutputForm _ for i in 0..min(4,#values-5)]) bracket(hconcat([bracket((displayVariable l)::OutputForm ** _ diff --git a/src/algebra/ssolve.spad.pamphlet b/src/algebra/ssolve.spad.pamphlet index cce6ca7..f761c4b 100644 --- a/src/algebra/ssolve.spad.pamphlet +++ b/src/algebra/ssolve.spad.pamphlet @@ -9,6 +9,156 @@ \begin{abstract} \end{abstract} \tableofcontents +\section{package UTSSOL TaylorSolve} + +[[UTSSOL]] is a facility to compute the first few coefficients of a Taylor +series given only implicitely by a function [[f]] that vanishes when applied to +the Taylor series. + +It uses the method of undetermined coefficients. + +\begin{ToDo} + Could I either + \begin{itemize} + \item take a function [[UTSCAT F -> UTSCAT F]] and still be able to compute + with undetermined coefficients, or + \item take a function [[F -> F]], and do likewise? + \end{itemize} + + Let's see. + + Try to compute the equation without resorting to power series. I.e., % + [[c: SUP SUP F]], and [[f: SUP SUP F -> SUP SUP F]]. Won't this make the + computation of coefficients terribly slow? + + I could also try to replace transcendental kernels with variables\dots + + Unfortunately, [[SUP F]] does not have [[TRANFUN]] -- well, it can't, of + course. However, I'd like to be able to compute % + [[sin(1+monomial(1,1)$UFPS SUP EXPR INT)]]. +\end{ToDo} + +<>= +)abb package UTSSOL TaylorSolve +TaylorSolve(F, UTSF, UTSSUPF): Exports == Implementation where + F: Field + SUP ==> SparseUnivariatePolynomialExpressions + UTSF: UnivariateTaylorSeriesCategory F + UTSSUPF: UnivariateTaylorSeriesCategory SUP F + NNI ==> NonNegativeInteger + + Exports == with + seriesSolve: (UTSSUPF -> UTSSUPF, List F) -> UTSF + + Implementation == add +<> +@ + +<>= + seriesSolve(f, l) == + c1 := map(#1::(SUP F), l)$ListFunctions2(F, SUP F)::(Stream SUP F) + coeffs: Stream SUP F := concat(c1, generate(monomial(1$F,1$NNI))) +-- coeffs: Stream SUP F := concat(c1, monomial(1$F,1$NNI)) +@ + +[[coeffs]] is the stream of the already computed coefficients of the solution, +plus one which is so far undetermined. We store in [[st.2]] the complete stream +and in [[st.1]] the stream starting with the first coefficient that has +possibly not yet been computed. + +\begin{ToDo} + The mathematics is not quite worked out. If [[coeffs]] is initialized as + stream with all coefficients set to the \emph{same} transcendental value, + and not enough initial values are given, then the missing ones are + implicitely assumed to be all identical. It may well happen that a solution + is produced, although it is not uniquely determined\dots +\end{ToDo} + +<>= + st: List Stream SUP F := [coeffs, coeffs] +@ + +Consider an arbitrary equation $f\big(x, y(x)\big)=0$. When setting $x=0$, we +obtain $f\big(0, y(0)\big)=0$. It is not necessarily the case that this +determines $y(0)$ uniquely, so we need one initial value that satisfies this +equation. +\begin{ToDo} + [[seriesSolve]] should check that the given initial values satisfy $f\big(0, y(0), + y'(0),...\big) = 0$. +\end{ToDo} +Now consider the derivatives of $f$, where we write $y$ instead of $y(x)$ for +better readability: +\begin{equation*} + \frac{d}{dx}f(x, y)=f_1(x, y) + f_2(x, y)y^\prime +\end{equation*} +and +\begin{align*} + \frac{d^2}{dx^2}f(x, y)&=f_{1,1}(x, y)\\ + &+f_{1,2}(x, y)y^\prime\\ + &+f_{2,1}(x, y)y^\prime\\ + &+f_{2,2}(x, y)(y^\prime)^2\\ + &+f_2(x, y)y^{\prime\prime}. +\end{align*} +In general, $\frac{d^2}{dx^2}f(x, y)$ depends only linearly on +$y^{\prime\prime}$. + +\begin{ToDo} + This points to another possibility: Since we know that we only need to solve + linear equations, we could compute two values and then use interpolation. + This might be a bit slower, but more importantly: can we still check that we + have enough initial values? Furthermore, we then really need that $f$ is + analytic, i.e., operators are not necessarily allowed anymore. However, it + seems that composition is allowed. +\end{ToDo} + +<>= + next: () -> F := + nr := st.1 + res: F + + if ground?(coeff: SUP F := nr.1)$(SUP F) +@ +%$ + +If the next element was already calculated, we can simply return it: + +<>= + then + res := ground coeff + st.1 := rest nr + else +@ + +Otherwise, we have to find the first non-satisfied relation and solve it. It +should be linear, or a single non-constant monomial. That is, the solution +should be unique. + +<>= + ns := st.2 + eqs: Stream SUP F := coefficients f series ns + while zero? first eqs repeat eqs := rest eqs + eq: SUP F := first eqs + if degree eq > 1 then + if monomial? eq then res := 0 + else + output(hconcat("The equation is: ", eq::OutputForm)) + $OutputPackage + error "seriesSolve: equation for coefficient not linear" + else res := (-coefficient(eq, 0$NNI)$(SUP F) + /coefficient(eq, 1$NNI)$(SUP F)) + + nr.1 := res::SUP F +-- concat!(st.2, monomial(1$F,1$NNI)) + st.1 := rest nr + + res + + series generate next + +@ +%$ + + \section{package EXPRSOL ExpressionSolve} \begin{ToDo} @@ -34,6 +184,7 @@ ExpressionSolve(R, F, UTSF, UTSSUPF): Exports == Implementation where Exports == with seriesSolve: (F, OP, SY, List F) -> UTSF + replaceDiffs: (F, OP, Symbol) -> F Implementation == add <> @@ -51,18 +202,25 @@ It turns out that the compiler doesn't find the right definition of that's even cleaner. Also, we need to tell the compiler that kernels that are independent of the main variable should be coerced to elements of the coefficient ring, since it will complain otherwise. +\begin{ToDo} + I cannot find an example for this behaviour right now. However, if I do use + the coerce, the following fails: + \begin{verbatim} + seriesSolve(h x -1-x*h x *h(q*x), h, x, [1]) + \end{verbatim} +\end{ToDo} <>= opelt := operator("elt"::Symbol)$OP opdiff := operator("D"::Symbol)$OP opcoerce := operator("coerce"::Symbol)$OP - replaceDiffs: (F, OP, Symbol) -> F +-- replaceDiffs: (F, OP, Symbol) -> F replaceDiffs (expr, op, sy) == lk := kernels expr for k in lk repeat - if freeOf?(coerce k, sy) then - expr := subst(expr, [k], [opcoerce [coerce k]]) +-- if freeOf?(coerce k, sy) then +-- expr := subst(expr, [k], [opcoerce [coerce k]]) if is?(k, op) then arg := first argument k @@ -139,10 +297,8 @@ works. This is probably due to missing [[/]] in [[UFPS]]. <<*>>= <> - -<> - -<> +<> +<> @ \end{document} diff --git a/src/axiom-website/patches.html b/src/axiom-website/patches.html index 3189281..b6d0a87 100644 --- a/src/axiom-website/patches.html +++ b/src/axiom-website/patches.html @@ -817,6 +817,8 @@ bookvol10.3 add domains
bookvol10.3 add domains
20081217.01.tpd.patch padic.spad removed
+20081217.02.tpd.patch +add guess package
- \ No newline at end of file + diff --git a/src/input/ndftip.input.pamphlet b/src/input/ndftip.input.pamphlet index 84578e4..3e4c0ab 100644 --- a/src/input/ndftip.input.pamphlet +++ b/src/input/ndftip.input.pamphlet @@ -586,7 +586,7 @@ hdftsD := nagHermitianDFT seqsD; --S 34 of 45 used to work? map(expand,hdftsD) :: List Vector Complex Float --R ---R There are 68 exposed and 8 unexposed library operations named map +--R There are 68 exposed and 9 unexposed library operations named map --R having 2 argument(s) but none was determined to be applicable. --R Use HyperDoc Browse, or issue --R )display op map