|
Here is the source code rec.spad
axiom )lib UFPS EXPRSOL SUPEXPR
UnivariateFormalPowerSeries is now explicitly exposed in frame
initial
UnivariateFormalPowerSeries will be automatically loaded when needed
from /var/zope2/var/LatexWiki/UFPS.NRLIB/code
ExpressionSolve is now explicitly exposed in frame initial
ExpressionSolve will be automatically loaded when needed from
/var/zope2/var/LatexWiki/EXPRSOL.NRLIB/code
SparseUnivariatePolynomialExpressions is now explicitly exposed in
frame initial
SparseUnivariatePolynomialExpressions will be automatically loaded
when needed from /var/zope2/var/LatexWiki/SUPEXPR.NRLIB/code
spad -- Copyright (C) 2006 Martin Rubey <Martin.Rubey@univie.ac.at> - -- This program is free software; you can redistribute it and/or
-- modify it under the terms of the GNU General Public License as
-- published by the Free Software Foundation; either version 2 of
-- the License, or (at your option) any later version.
--
-- This program is distributed in the hope that it will be
-- useful, but WITHOUT ANY WARRANTY; without even the implied
-- warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
-- PURPOSE. See the GNU General Public License for more details.
)abbrev package RECOP RecurrenceOperator
++ Author: Martin Rubey
++ Description:
++ This package provides an operator for the n-th term of a recurrence and an
++ operator for the coefficient of x^n in a function specified by a functional
++ equation.
RecurrenceOperator(R, F): Exports == Implementation where
R: Join(OrderedSet, IntegralDomain, ConvertibleTo InputForm)
F: Join(FunctionSpace R, AbelianMonoid, RetractableTo Integer, _
RetractableTo Symbol, PartialDifferentialRing Symbol)
--RecurrenceOperator(F): Exports == Implementation where
-- F: Join(ExpressionSpace, AbelianMonoid, RetractableTo Integer,
-- RetractableTo Symbol, PartialDifferentialRing Symbol)
Exports == with
evalRec: (BasicOperator, Symbol, F, F, F, List F) -> F
++ \spad{evalRec(u, dummy, n, n0, eq, values)} creates an expression that
++ stands for u(n0), where u(n) is given by the equation eq. However, for
++ technical reasons the variable n has to be replaced by a dummy
++ variable dummy in eq. The argument values specifies the initial values
++ of the recurrence u(0), u(1),...
++ For the moment we don't allow recursions that contain u inside of
++ another operator.
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),...
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
-- should be local
if R has Ring then
getShiftRec: (BasicOperator, Kernel F, Symbol) -> Union(Integer, "failed")
shiftInfoRec: (BasicOperator, Symbol, F) ->
Record(max: Union(Integer, "failed"),
ord: Union(Integer, "failed"),
ker: Kernel F)
Implementation == add
oprecur := operator("rootOfRec"::Symbol)$BasicOperator
opADE := operator("rootOfADE"::Symbol)$BasicOperator
setProperty(oprecur, "%dummyVar", 2 pretend None)
setProperty(opADE, "%dummyVar", 2 pretend None)
eqAsF: List F -> F
eqAsF l == l.1
dummy: List F -> Symbol
dummy l == retract(l.2)@Symbol
dummyAsF: List F -> F
dummyAsF l == l.2
displayVariable: List F -> F
displayVariable l == l.3
operatorName: List F -> BasicOperator
operatorName l == operator(kernels(l.4).1)
operatorNameAsF: List F -> F
operatorNameAsF l == l.4
operatorArgument: List F -> F
operatorArgument l == argument(kernels(l.4).1).1
initialValues: List F -> List F
initialValues l == rest(l, 4)
if R has Ring then
getShiftRec(op: BasicOperator, f: Kernel F, n: Symbol)
: Union(Integer, "failed") ==
a := argument f
if every?(freeOf?(#1, n::F), a) then return 0
if #a ~= 1 then error "RECOP: operator should have only one argument"
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
and every?(freeOf?(#1, n::F), coefficients num)
then return retractIfCan(coefficient(num, 0))
else return "failed"
-- if the recurrence is of the form
-- $p(n, f(n+m-o), f(n+m-o+1), \dots, f(n+m)) = 0$
-- in which case shiftInfoRec returns [m, o, f(n+m)].
shiftInfoRec(op: BasicOperator, argsym: Symbol, eq: F):
Record(max: Union(Integer, "failed"),
ord: Union(Integer, "failed"),
ker: Kernel F) ==
-- ord and ker are valid only if all shifts are Integers
-- ker is the kernel of the maximal shift.
maxShift: Integer
minShift: Integer
nextKernel: Kernel F
-- We consider only those kernels that have op as operator. If there is none,
-- we raise an error. For the moment we don't allow recursions that contain op
-- inside of another operator.
error? := true
for f in kernels eq repeat
if is?(f, op) then
shift := getShiftRec(op, f, argsym)
if error? then
error? := false
nextKernel := f
if shift case Integer then
maxShift := shift
minShift := shift
else return ["failed", "failed", nextKernel]
else
if shift case Integer then
if maxShift < shift then
maxShift := shift
nextKernel := f
if minShift > shift then
minShift := shift
else return ["failed", "failed", nextKernel]
if error? then error "evalRec: equation does not contain operator"
[maxShift, maxShift - minShift, nextKernel]
evalRec(op, argsym, argdisp, arg, eq, values) ==
if ((n := retractIfCan(arg)@Union(Integer, "failed")) case "failed")
or (n < 0)
then
shiftInfo := shiftInfoRec(op, argsym, eq)
if (shiftInfo.ord case "failed") or ((shiftInfo.ord)::Integer > 0)
then
kernel(oprecur,
append([eq, argsym::F, argdisp, op(arg)], values))
else
p := univariate(eq, shiftInfo.ker)
num := numer p
-- If the degree is 1, we can return the function explicitly.
if degree num = 1 then
eval(-coefficient(num, 0)/coefficient(num, 1), argsym::F,
arg::F-(shiftInfo.max)::Integer::F)
else
kernel(oprecur,
append([eq, argsym::F, argdisp, op(arg)], values))
else
len: Integer := #values
if n < len
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)
newval := eval(nextval, op,
evalRec(op, argsym, argdisp, #1, eq, values))
evalRec(op, argsym, argdisp, arg, eq, cons(newval, values))
else
kernel(oprecur,
append([eq, argsym::F, argdisp, op(arg)], values))
else
kernel(oprecur,
append([eq, argsym::F, argdisp, op(arg)], values))
numberOfValuesNeeded(numberOfValues: Integer,
op: BasicOperator, argsym: Symbol, eq: F): Integer ==
order := shiftInfoRec(op, argsym, eq).ord
if order case Integer
then min(numberOfValues, retract(order)@Integer)
else numberOfValues
else
evalRec(op, argsym, argdisp, arg, eq, values) ==
kernel(oprecur,
append([eq, argsym::F, argdisp, op(arg)], values))
numberOfValuesNeeded(numberOfValues: Integer,
op: BasicOperator, argsym: Symbol, eq: F): Integer ==
numberOfValues
irecur: List F -> F
irecur l ==
evalRec(operatorName l,
dummy l, displayVariable l,
operatorArgument l, eqAsF l, initialValues l)
evaluate(oprecur, irecur)$BasicOperatorFunctions1(F)
ddrec: List F -> OutputForm
ddrec l ==
op := operatorName l
values := reverse l
eq := eqAsF l
numberOfValues := numberOfValuesNeeded(#values-4, op, dummy l, eq)
vals: List OutputForm
:= cons(eval(eq, dummyAsF l, displayVariable l)::OutputForm = _
0::OutputForm,
[elt(op::OutputForm, [(i-1)::OutputForm]) = _
(values.i)::OutputForm
for i in 1..numberOfValues])
bracket(hconcat([(operatorNameAsF l)::OutputForm,
": ",
commaSeparate vals]))
setProperty(oprecur, "%specialDisp",
ddrec@(List F -> OutputForm) pretend None)
getOrder(op: BasicOperator, f: Kernel F): NonNegativeInteger ==
res: NonNegativeInteger := 0
g := f
while is?(g, %diff) repeat
g := kernels(argument(g).1).1
res := res+1
if is?(g, op) then res else 0
evalADE(op, argsym, argdisp, arg, eq, values) ==
if not freeOf?(eq, retract(argdisp)@Symbol)
then error "RECOP: The argument should not be used in the equation of the_
ADE"
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))
else
if n < 0
then 0
else
keq := kernels eq
order := getOrder(op, keq.1)
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
if degree numer p > 1
then
kernel(opADE,
append([eq, argsym::F, argdisp, op(arg)], values))
else
s := seriesSolve(eq, op, argsym, first(reverse values, order))
$ExpressionSolve(R, F,
UnivariateFormalPowerSeries F,
UnivariateFormalPowerSeries
SparseUnivariatePolynomialExpressions F)
elt(s, n::Integer::NonNegativeInteger)
iADE: List F -> F
-- This is just a wrapper that allows us to write a recurrence relation as an
-- operator.
iADE l ==
evalADE(operatorName l,
dummy l, displayVariable l,
operatorArgument l, eqAsF l, initialValues l)
evaluate(opADE, iADE)$BasicOperatorFunctions1(F)
getADE(f: F): F ==
ker := kernels f
if one?(#ker) and is?(operator(ker.1), "rootOfADE"::Symbol) then
l := argument(ker.1)
eval(eqAsF l, dummyAsF l, displayVariable l)
else
error "getADE: argument should be a single rootOfADE object"
ddADE: List F -> OutputForm
ddADE l ==
op := operatorName l
values := reverse l
vals: List OutputForm
:= 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
for i in 0..min(4,#values-5)])
bracket(hconcat([bracket((displayVariable l)::OutputForm ** _
(operatorArgument l)::OutputForm),
(op(displayVariable l))::OutputForm, ": ",
commaSeparate vals]))
setProperty(opADE, "%specialDisp",
ddADE@(List F -> OutputForm) pretend None)
spad Compiling FriCAS source code from file
/var/zope2/var/LatexWiki/496568161422582843-25px002.spad using
old system compiler.
RECOP abbreviates package RecurrenceOperator
------------------------------------------------------------------------
initializing NRLIB RECOP for RecurrenceOperator
compiling into NRLIB RECOP
****** Domain: R already in scope
****** Domain: R already in scope
****** Domain: F already in scope
****** Domain: F already in scope
****** Domain: F already in scope
****** Domain: F already in scope
compiling local eqAsF : List F -> F
Time: 0.21 SEC.
compiling local dummy : List F -> Symbol
Time: 0.02 SEC.
compiling local dummyAsF : List F -> F
Time: 0.01 SEC.
compiling local displayVariable : List F -> F
Time: 0.02 SEC.
compiling local operatorName : List F -> BasicOperator
Time: 0.10 SEC.
compiling local operatorNameAsF : List F -> F
Time: 0.02 SEC.
compiling local operatorArgument : List F -> F
Time: 0.10 SEC.
compiling local initialValues : List F -> List F
Time: 0.02 SEC.
compiling exported getShiftRec : (BasicOperator,Kernel F,Symbol) -> Union(Integer,failed)
Time: 0.46 SEC.
compiling exported shiftInfoRec : (BasicOperator,Symbol,F) -> Record(max: Union(Integer,failed),ord: Union(Integer,failed),ker: Kernel F)
Time: 0.05 SEC.
compiling exported evalRec : (BasicOperator,Symbol,F,F,F,List F) -> F
Time: 0.63 SEC.
compiling exported numberOfValuesNeeded : (Integer,BasicOperator,Symbol,F) -> Integer
Time: 0.01 SEC.
compiling exported evalRec : (BasicOperator,Symbol,F,F,F,List F) -> F
Time: 0.02 SEC.
compiling exported numberOfValuesNeeded : (Integer,BasicOperator,Symbol,F) -> Integer
RECOP;numberOfValuesNeeded;IBoSFI;14 is replaced by numberOfValues
Time: 0 SEC.
compiling local irecur : List F -> F
Time: 0.02 SEC.
compiling local ddrec : List F -> OutputForm
Time: 0.21 SEC.
compiling local getOrder : (BasicOperator,Kernel F) -> NonNegativeInteger
Time: 0.18 SEC.
compiling exported evalADE : (BasicOperator,Symbol,F,F,F,List F) -> F
Time: 0.43 SEC.
compiling local iADE : List F -> F
Time: 0.02 SEC.
compiling exported getADE : F -> F
Time: 0.12 SEC.
compiling local ddADE : List F -> OutputForm
Time: 0.63 SEC.
(time taken in buildFunctor: 0)
;;; *** |RecurrenceOperator| REDEFINED
;;; *** |RecurrenceOperator| REDEFINED
Time: 0.06 SEC.
Warnings:
[1] operatorName: not known that (OrderedSet) is of mode (CATEGORY package (SIGNATURE evalRec (F (BasicOperator) (Symbol) F F F (List F))) (SIGNATURE evalADE (F (BasicOperator) (Symbol) F F F (List F))) (SIGNATURE getADE (F F)) (SIGNATURE numberOfValuesNeeded ((Integer) (Integer) (BasicOperator) (Symbol) F)) (IF (has R (Ring)) (PROGN (SIGNATURE getShiftRec ((Union (Integer) failed) (BasicOperator) (Kernel F) (Symbol))) (SIGNATURE shiftInfoRec ((Record (: max (Union (Integer) failed)) (: ord (Union (Integer) failed)) (: ker (Kernel F))) (BasicOperator) (Symbol) F))) noBranch))
[2] getShiftRec: not known that (OrderedSet) is of mode (CATEGORY package (SIGNATURE evalRec (F (BasicOperator) (Symbol) F F F (List F))) (SIGNATURE evalADE (F (BasicOperator) (Symbol) F F F (List F))) (SIGNATURE getADE (F F)) (SIGNATURE numberOfValuesNeeded ((Integer) (Integer) (BasicOperator) (Symbol) F)) (IF (has R (Ring)) (PROGN (SIGNATURE getShiftRec ((Union (Integer) failed) (BasicOperator) (Kernel F) (Symbol))) (SIGNATURE shiftInfoRec ((Record (: max (Union (Integer) failed)) (: ord (Union (Integer) failed)) (: ker (Kernel F))) (BasicOperator) (Symbol) F))) noBranch))
[3] getShiftRec: not known that (Ring) is of mode (CATEGORY package (SIGNATURE evalRec (F (BasicOperator) (Symbol) F F F (List F))) (SIGNATURE evalADE (F (BasicOperator) (Symbol) F F F (List F))) (SIGNATURE getADE (F F)) (SIGNATURE numberOfValuesNeeded ((Integer) (Integer) (BasicOperator) (Symbol) F)) (IF (has R (Ring)) (PROGN (SIGNATURE getShiftRec ((Union (Integer) failed) (BasicOperator) (Kernel F) (Symbol))) (SIGNATURE shiftInfoRec ((Record (: max (Union (Integer) failed)) (: ord (Union (Integer) failed)) (: ker (Kernel F))) (BasicOperator) (Symbol) F))) noBranch))
[4] shiftInfoRec: maxShift has no value
[5] shiftInfoRec: minShift has no value
[6] shiftInfoRec: nextKernel has no value
[7] evalRec: ord has no value
[8] evalRec: ker has no value
[9] getOrder: %diff has no value
Cumulative Statistics for Constructor RecurrenceOperator
Time: 3.34 seconds
finalizing NRLIB RECOP
Processing RecurrenceOperator for Browser database:
--------(evalRec (F (BasicOperator) (Symbol) F F F (List F)))---------
--------(evalADE (F (BasicOperator) (Symbol) F F F (List F)))---------
--------(getADE (F F))---------
--->/usr/local/lib/axiom/target/x86_64-unknown-linux/../../src/algebra/RECOP.spad-->RecurrenceOperator((numberOfValuesNeeded ((Integer) (Integer) (BasicOperator) (Symbol) F))): Not documented!!!!
--->/usr/local/lib/axiom/target/x86_64-unknown-linux/../../src/algebra/RECOP.spad-->RecurrenceOperator((getShiftRec ((Union (Integer) failed) (BasicOperator) (Kernel F) (Symbol)))): Not documented!!!!
--->/usr/local/lib/axiom/target/x86_64-unknown-linux/../../src/algebra/RECOP.spad-->RecurrenceOperator((shiftInfoRec ((Record (: max (Union (Integer) failed)) (: ord (Union (Integer) failed)) (: ker (Kernel F))) (BasicOperator) (Symbol) F))): Not documented!!!!
--------constructor---------
------------------------------------------------------------------------
RecurrenceOperator is now explicitly exposed in frame initial
RecurrenceOperator will be automatically loaded when needed from
/var/zope2/var/LatexWiki/RECOP.NRLIB/code
|