|
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/UFPS
ExpressionSolve is now explicitly exposed in frame initial
ExpressionSolve will be automatically loaded when needed from
/var/zope2/var/LatexWiki/EXPRSOL.NRLIB/EXPRSOL
SparseUnivariatePolynomialExpressions is now explicitly exposed in
frame initial
SparseUnivariatePolynomialExpressions will be automatically loaded
when needed from /var/zope2/var/LatexWiki/SUPEXPR.NRLIB/SUPEXPR
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
****** comp fails at level 2 with expression: ******
(|setProperty| |oprecur| | << %dummyVar >> | (|pretend| 2 (|None|)))
****** level 2 ******
$x:= %dummyVar
$m:= (Symbol)
$f:=
((((* #) (+ #) (< #) (<= #) ...)))
>> Apparent user error:
Cannot coerce %dummyVar
of mode %dummyVar
to mode (Union %dummyVar (Symbol))
|
