Changes http://wiki.axiom-developer.org/RecSpad/diff
--
??changed:
-This page was renamed to rec.spad.pamphlet. You can delete this one if no
longer needed.
-
+Here is the source code 'rec.spad'
+
+\begin{spad}
+)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),...
+
+-- 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)
+
+-- this implies that the second and third arguments of oprecur are dummy
+-- variables and affects tower$ES:
+-- the second argument will not appear in tower$ES, if it does not appear in
+-- any argument but the first and second.
+-- the third argument will not appear in tower$ES, if it does not appear in any
+-- other argument. (%defsum is a good example)
+
+-- The arguments of these operators are as follows:
+
+-- 1. eq, i.e. the vanishing expression
+
+ eqAsF: List F -> F
+ eqAsF l == l.1
+
+-- 2. dummy
+
+ dummy: List F -> Symbol
+ dummy l == retract(l.2)@Symbol
+
+ dummyAsF: List F -> F
+ dummyAsF l == l.2
+
+-- 3. variable, i.e., for display
+
+ displayVariable: List F -> F
+ displayVariable l == l.3
+
+-- 4. operatorname(argument)
+
+ 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
+
+-- rootOfADE: note that although we have arg as argument of the operator,
+-- it is intended to indicate the coefficient, not the argument
+-- of the power series
+
+-- 5.- values in reversed order
+-- rootOfRec:
+-- maybe "values" should be preceded by the index of the first given
+-- value. Currently, the last value is interpreted as f(0)
+
+-- rootOfADE:
+-- values are the first few coefficients of the power series expansion in order
+
+ 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"
+
+
+ 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 non, 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]
+
+-- try to evaluate:
+ 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 degree num = 1 then
+ eval(-coefficient(num, 0)/coefficient(num, 1), argsym::F, arg::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
+
+[143 more lines...]
--
forwarded from http://wiki.axiom-developer.org/[EMAIL PROTECTED]
