open-axiom repository from github
\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{\$SPAD/src/algebra combfunc.spad}
\author{Manuel Bronstein, Martin Rubey}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{category COMBOPC CombinatorialOpsCategory}
<<category COMBOPC CombinatorialOpsCategory>>=
)abbrev category COMBOPC CombinatorialOpsCategory
++ Category for summations and products
++ Author: Manuel Bronstein
++ Date Created: ???
++ Date Last Updated: 22 February 1993 (JHD/BMT)
++ Description:
++ CombinatorialOpsCategory is the category obtaining by adjoining
++ summations and products to the usual combinatorial operations;
CombinatorialOpsCategory(): Category ==
CombinatorialFunctionCategory with
factorials : $ -> $
++ factorials(f) rewrites the permutations and binomials in f
++ in terms of factorials;
factorials : ($, Symbol) -> $
++ factorials(f, x) rewrites the permutations and binomials in f
++ involving x in terms of factorials;
summation : ($, Symbol) -> $
++ summation(f(n), n) returns the formal sum S(n) which verifies
++ S(n+1) - S(n) = f(n);
summation : ($, SegmentBinding $) -> $
++ summation(f(n), n = a..b) returns f(a) + ... + f(b) as a
++ formal sum;
product : ($, Symbol) -> $
++ product(f(n), n) returns the formal product P(n) which verifies
++ P(n+1)/P(n) = f(n);
product : ($, SegmentBinding $) -> $
++ product(f(n), n = a..b) returns f(a) * ... * f(b) as a
++ formal product;
@
The latest change allows Axiom to reduce
\begin{verbatim}
sum(1/i,i=1..n)-sum(1/i,i=1..n)
\end{verbatim}
to reduce to zero.
<<package COMBF CombinatorialFunction>>=
)abbrev package COMBF CombinatorialFunction
++ Provides the usual combinatorial functions
++ Author: Manuel Bronstein, Martin Rubey
++ Date Created: 2 Aug 1988
++ Date Last Updated: 30 October 2005
++ Description:
++ Provides combinatorial functions over an integral domain.
++ Keywords: combinatorial, function, factorial.
++ Examples: )r COMBF INPUT
CombinatorialFunction(R, F): Exports == Implementation where
R: IntegralDomain
F: FunctionSpace R
OP ==> BasicOperator
K ==> Kernel F
SE ==> Symbol
O ==> OutputForm
SMP ==> SparseMultivariatePolynomial(R, K)
Z ==> Integer
POWER ==> '%power
OPEXP ==> 'exp
SPECIALDIFF ==> "%specialDiff"
SPECIALDISP ==> "%specialDisp"
SPECIALEQUAL ==> "%specialEqual"
Exports ==> with
belong? : OP -> Boolean
++ belong?(op) is true if op is a combinatorial operator;
operator : OP -> OP
++ operator(op) returns a copy of op with the domain-dependent
++ properties appropriate for F;
++ error if op is not a combinatorial operator;
** : (F, F) -> F
++ a ** b is the formal exponential a**b;
binomial : (F, F) -> F
++ binomial(n, r) returns the number of subsets of r objects
++ taken among n objects, i.e. n!/(r! * (n-r)!);
permutation: (F, F) -> F
++ permutation(n, r) returns the number of permutations of
++ n objects taken r at a time, i.e. n!/(n-r)!;
factorial : F -> F
++ factorial(n) returns the factorial of n, i.e. n!;
factorials : F -> F
++ factorials(f) rewrites the permutations and binomials in f
++ in terms of factorials;
factorials : (F, SE) -> F
++ factorials(f, x) rewrites the permutations and binomials in f
++ involving x in terms of factorials;
summation : (F, SE) -> F
++ summation(f(n), n) returns the formal sum S(n) which verifies
++ S(n+1) - S(n) = f(n);
summation : (F, SegmentBinding F) -> F
++ summation(f(n), n = a..b) returns f(a) + ... + f(b) as a
++ formal sum;
product : (F, SE) -> F
++ product(f(n), n) returns the formal product P(n) which verifies
++ P(n+1)/P(n) = f(n);
product : (F, SegmentBinding F) -> F
++ product(f(n), n = a..b) returns f(a) * ... * f(b) as a
++ formal product;
iifact : F -> F
++ iifact(x) should be local but conditional;
iibinom : List F -> F
++ iibinom(l) should be local but conditional;
iiperm : List F -> F
++ iiperm(l) should be local but conditional;
iipow : List F -> F
++ iipow(l) should be local but conditional;
iidsum : List F -> F
++ iidsum(l) should be local but conditional;
iidprod : List F -> F
++ iidprod(l) should be local but conditional;
ipow : List F -> F
++ ipow(l) should be local but conditional;
Implementation ==> add
ifact : F -> F
iiipow : List F -> F
iperm : List F -> F
ibinom : List F -> F
isum : List F -> F
idsum : List F -> F
iprod : List F -> F
idprod : List F -> F
dsum : List F -> O
ddsum : List F -> O
dprod : List F -> O
ddprod : List F -> O
equalsumprod : (K, K) -> Boolean
equaldsumprod : (K, K) -> Boolean
fourth : List F -> F
dvpow1 : List F -> F
dvpow2 : List F -> F
summand : List F -> F
dvsum : (List F, SE) -> F
dvdsum : (List F, SE) -> F
dvprod : (List F, SE) -> F
dvdprod : (List F, SE) -> F
facts : (F, List SE) -> F
K2fact : (K, List SE) -> F
smpfact : (SMP, List SE) -> F
dummy == new()$SE :: F
@
This macro will be used in [[product]] and [[summation]], both the $5$ and $3$
argument forms. It is used to introduce a dummy variable in place of the
summation index within the summands. This in turn is necessary to keep the
indexing variable local, circumventing problems, for example, with
differentiation.
This works if we don't accidently use such a symbol as a bound of summation or
product.
Note that up to [[patch--25]] this used to read
\begin{verbatim}
dummy := new()$SE :: F
\end{verbatim}
thus introducing the same dummy variable for all products and summations, which
caused nested products and summations to fail. (Issue~\#72)
<<package COMBF CombinatorialFunction>>=
opfact := operator('factorial)$CommonOperators
opperm := operator('permutation)$CommonOperators
opbinom := operator('binomial)$CommonOperators
opsum := operator('summation)$CommonOperators
opdsum := operator('%defsum)$CommonOperators
opprod := operator('product)$CommonOperators
opdprod := operator('%defprod)$CommonOperators
oppow := operator(POWER)$CommonOperators
factorial x == opfact x
binomial(x, y) == opbinom [x, y]
permutation(x, y) == opperm [x, y]
import F
import Kernel F
number?(x:F):Boolean ==
if R has RetractableTo(Z) then
ground?(x) or
((retractIfCan(x)@Union(Fraction(Z),"failed")) case Fraction(Z))
else
ground?(x)
x ** y ==
-- Do some basic simplifications
is?(x,POWER) =>
args : List F := argument first kernels x
not(#args = 2) => error "Too many arguments to **"
number?(first args) and number?(y) =>
oppow [first(args)**y, second args]
oppow [first args, (second args)* y]
-- Generic case
exp : Union(Record(val:F,exponent:Z),"failed") := isPower x
exp case Record(val:F,exponent:Z) =>
expr := exp::Record(val:F,exponent:Z)
oppow [expr.val, (expr.exponent)*y]
oppow [x, y]
belong? op == has?(op, 'comb)
fourth l == third rest l
dvpow1 l == second(l) * first(l) ** (second l - 1)
factorials x == facts(x, variables x)
factorials(x, v) == facts(x, [v])
facts(x, l) == smpfact(numer x, l) / smpfact(denom x, l)
summand l == eval(first l, retract(second l)@K, third l)
product(x:F, i:SE) ==
dm := dummy
opprod [eval(x, k := kernel(i)$K, dm), dm, k::F]
summation(x:F, i:SE) ==
dm := dummy
opsum [eval(x, k := kernel(i)$K, dm), dm, k::F]
@
These two operations return the product or the sum as unevaluated operators. A
dummy variable is introduced to make the indexing variable \lq local\rq.
<<package COMBF CombinatorialFunction>>=
dvsum(l, x) ==
opsum [differentiate(first l, x), second l, third l]
dvdsum(l, x) ==
x = retract(y := third l)@SE => 0
if member?(x, variables(h := third rest rest l)) or
member?(x, variables(g := third rest l)) then
error "a sum cannot be differentiated with respect to a bound"
else
opdsum [differentiate(first l, x), second l, y, g, h]
@
The above two operations implement differentiation of sums with and without
bounds. Note that the function
$$n\mapsto\sum_{k=1}^n f(k,n)$$
is well defined only for integral values of $n$ greater than or equal to zero.
There is not even consensus how to define this function for $n<0$. Thus, it is
not differentiable. Therefore, we need to check whether we erroneously are
differentiating with respect to the upper bound or the lower bound, where the
same reasoning holds.
Differentiating a sum with respect to its indexing variable correctly gives
zero. This is due to the introduction of dummy variables in the internal
representation of a sum: the operator [[%defsum]] takes 5 arguments, namely
\begin{enumerate}
\item the summands, where each occurrence of the indexing variable is replaced
by
\item the dummy variable,
\item the indexing variable,
\item the lower bound, and
\item the upper bound.
\end{enumerate}
Note that up to [[patch--40]] the following incorrect code was used, which
tried to parallel the known rules for integration: (Issue~\#180)
\begin{verbatim}
dvdsum(l, x) ==
x = retract(y := third l)@SE => 0
k := retract(d := second l)@K
differentiate(h := third rest rest l,x) * eval(f := first l, k, h)
- differentiate(g := third rest l, x) * eval(f, k, g)
+ opdsum [differentiate(f, x), d, y, g, h]
\end{verbatim}
Up to [[patch--45]] a similar mistake could be found in the code for
differentiation of formal sums, which read
\begin{verbatim}
dvsum(l, x) ==
k := retract(second l)@K
differentiate(third l, x) * summand l
+ opsum [differentiate(first l, x), second l, third l]
\end{verbatim}
<<package COMBF CombinatorialFunction>>=
dvprod(l, x) ==
dm := retract(dummy)@SE
f := eval(first l, retract(second l)@K, dm::F)
p := product(f, dm)
opsum [differentiate(first l, x)/first l * p, second l, third l]
dvdprod(l, x) ==
x = retract(y := third l)@SE => 0
if member?(x, variables(h := third rest rest l)) or
member?(x, variables(g := third rest l)) then
error "a product cannot be differentiated with respect to a bound"
else
opdsum cons(differentiate(first l, x)/first l, rest l) * opdprod l
@
The above two operations implement differentiation of products with and without
bounds. Note again, that we cannot even properly define products with bounds
that are not integral.
To differentiate the product, we use Leibniz rule:
$$\frac{d}{dx}\prod_{i=a}^b f(i,x) =
\sum_{i=a}^b \frac{\frac{d}{dx} f(i,x)}{f(i,x)}\prod_{i=a}^b f(i,x)
$$
There is one situation where this definition might produce wrong results,
namely when the product is zero, but axiom failed to recognize it: in this
case,
$$
\frac{d}{dx} f(i,x)/f(i,x)
$$
is undefined for some $i$. However, I was not able to come up with an
example. The alternative definition
$$
\frac{d}{dx}\prod_{i=a}^b f(i,x) =
\sum_{i=a}^b \left(\frac{d}{dx} f(i,x)\right)\prod_{j=a,j\neq i}^b f(j,x)
$$
has the slight (display) problem that we would have to come up with a new index
variable, which looks very ugly. Furthermore, it seems to me that more
simplifications will occur with the first definition.
<<TEST COMBF>>=
f := operator 'f
D(product(f(i,x),i=1..m),x)
@
Note that up to [[patch--45]] these functions did not exist and products were
differentiated according to the usual chain rule, which gave incorrect
results. (Issue~\#211)
<<package COMBF CombinatorialFunction>>=
dprod l ==
prod(summand(l)::O, third(l)::O)
ddprod l ==
prod(summand(l)::O, third(l)::O = fourth(l)::O, fourth(rest l)::O)
dsum l ==
sum(summand(l)::O, third(l)::O)
ddsum l ==
sum(summand(l)::O, third(l)::O = fourth(l)::O, fourth(rest l)::O)
@
These four operations handle the conversion of sums and products to
[[OutputForm]]. Note that up to [[patch--45]] the definitions for sums and
products without bounds were missing and output was illegible.
<<package COMBF CombinatorialFunction>>=
equalsumprod(s1, s2) ==
l1 := argument s1
l2 := argument s2
(eval(first l1, retract(second l1)@K, second l2) = first l2)
equaldsumprod(s1, s2) ==
l1 := argument s1
l2 := argument s2
((third rest l1 = third rest l2) and
(third rest rest l1 = third rest rest l2) and
(eval(first l1, retract(second l1)@K, second l2) = first l2))
@
The preceding two operations handle the testing for equality of sums and
products. This functionality was missing up to [[patch--45]]. (Issue~\#213) The
corresponding property [[%specialEqual]] set below is checked in
[[Kernel]]. Note that we can assume that the operators are equal, since this is
checked in [[Kernel]] itself.
<<package COMBF CombinatorialFunction>>=
product(x:F, s:SegmentBinding F) ==
k := kernel(variable s)$K
dm := dummy
opdprod [eval(x,k,dm), dm, k::F, lo segment s, hi segment s]
summation(x:F, s:SegmentBinding F) ==
k := kernel(variable s)$K
dm := dummy
opdsum [eval(x,k,dm), dm, k::F, lo segment s, hi segment s]
@
These two operations return the product or the sum as unevaluated operators. A
dummy variable is introduced to make the indexing variable \lq local\rq.
<<package COMBF CombinatorialFunction>>=
smpfact(p, l) ==
map(K2fact(#1, l), #1::F, p)$PolynomialCategoryLifting(
IndexedExponents K, K, R, SMP, F)
K2fact(k, l) ==
kf : F
empty? [v for v in variables(kf := k::F) | member?(v, l)] => kf
empty?(args:List F := [facts(a, l) for a in argument k]) => kf
is?(k, opperm) =>
factorial(n := first args) / factorial(n - second args)
is?(k, opbinom) =>
n := first args
p := second args
factorial(n) / (factorial(p) * factorial(n-p))
(operator k) args
operator op ==
is?(op,'factorial) => opfact
is?(op,'permutation) => opperm
is?(op,'binomial) => opbinom
is?(op,'summation) => opsum
is?(op,'%defsum) => opdsum
is?(op,'product) => opprod
is?(op,'%defprod) => opdprod
is?(op, POWER) => oppow
error "Not a combinatorial operator"
iprod l ==
zero? first l => 0
one? first l => 1
kernel(opprod, l)
isum l ==
zero? first l => 0
kernel(opsum, l)
idprod l ==
member?(retract(second l)@SE, variables first l) =>
kernel(opdprod, l)
first(l) ** (fourth rest l - fourth l + 1)
idsum l ==
member?(retract(second l)@SE, variables first l) =>
kernel(opdsum, l)
first(l) * (fourth rest l - fourth l + 1)
ifact x ==
zero? x or one? x => 1
kernel(opfact, x)
ibinom l ==
n := first l
((p := second l) = 0) or (p = n) => 1
one? p or (p = n - 1) => n
kernel(opbinom, l)
iperm l ==
zero? second l => 1
kernel(opperm, l)
if R has RetractableTo Z then
iidsum l ==
(r1:=retractIfCan(fourth l)@Union(Z,"failed"))
case "failed" or
(r2:=retractIfCan(fourth rest l)@Union(Z,"failed"))
case "failed" or
(k:=retractIfCan(second l)@Union(K,"failed")) case "failed"
=> idsum l
+/[eval(first l,k::K,i::F) for i in r1::Z .. r2::Z]
iidprod l ==
(r1:=retractIfCan(fourth l)@Union(Z,"failed"))
case "failed" or
(r2:=retractIfCan(fourth rest l)@Union(Z,"failed"))
case "failed" or
(k:=retractIfCan(second l)@Union(K,"failed")) case "failed"
=> idprod l
*/[eval(first l,k::K,i::F) for i in r1::Z .. r2::Z]
iiipow l ==
(u := isExpt(x := first l, OPEXP)) case "failed" => kernel(oppow, l)
rec := u::Record(var: K, exponent: Z)
y := first argument(rec.var)
(r := retractIfCan(y)@Union(Fraction Z, "failed")) case
"failed" => kernel(oppow, l)
(operator(rec.var)) (rec.exponent * y * second l)
if F has RadicalCategory then
ipow l ==
(r := retractIfCan(second l)@Union(Fraction Z,"failed"))
case "failed" => iiipow l
first(l) ** (r::Fraction(Z))
else
ipow l ==
(r := retractIfCan(second l)@Union(Z, "failed"))
case "failed" => iiipow l
first(l) ** (r::Z)
else
ipow l ==
zero?(x := first l) =>
zero? second l => error "0 ** 0"
0
one? x or zero?(n := second l) => 1
one? n => x
(u := isExpt(x, OPEXP)) case "failed" => kernel(oppow, l)
rec := u::Record(var: K, exponent: Z)
one?(y := first argument(rec.var)) or y = -1 =>
(operator(rec.var)) (rec.exponent * y * n)
kernel(oppow, l)
if R has CombinatorialFunctionCategory then
iifact x ==
(r:=retractIfCan(x)@Union(R,"failed")) case "failed" => ifact x
factorial(r::R)::F
iiperm l ==
(r1 := retractIfCan(first l)@Union(R,"failed")) case "failed" or
(r2 := retractIfCan(second l)@Union(R,"failed")) case "failed"
=> iperm l
permutation(r1::R, r2::R)::F
if R has RetractableTo(Z) and F has Algebra(Fraction(Z)) then
iibinom l ==
(s:=retractIfCan(first l-second l)@Union(R,"failed")) case R and
(t:=retractIfCan(s)@Union(Z,"failed")) case Z and positive? t =>
ans:=1::F
for i in 1..t repeat
ans:=ans*(second l+i::R::F)
(1/factorial t) * ans
(r1 := retractIfCan(first l)@Union(R,"failed")) case "failed" or
(r2 := retractIfCan(second l)@Union(R,"failed")) case "failed"
=> ibinom l
binomial(r1::R, r2::R)::F
else
iibinom l ==
(r1 := retractIfCan(first l)@Union(R,"failed")) case "failed" or
(r2 := retractIfCan(second l)@Union(R,"failed")) case "failed"
=> ibinom l
binomial(r1::R, r2::R)::F
else
iifact x == ifact x
iibinom l == ibinom l
iiperm l == iperm l
if R has ElementaryFunctionCategory then
iipow l ==
(r1:=retractIfCan(first l)@Union(R,"failed")) case "failed" or
(r2:=retractIfCan(second l)@Union(R,"failed")) case "failed"
=> ipow l
(r1::R ** r2::R)::F
else
iipow l == ipow l
if F has ElementaryFunctionCategory then
dvpow2 l == if zero?(first l) then
0
else
log(first l) * first(l) ** second(l)
@
This operation implements the differentiation of the power operator [[%power]]
with respect to its second argument, i.e., the exponent. It uses the formula
$$\frac{d}{dx} g(y)^x = \frac{d}{dx} e^{x\log g(y)} = \log g(y) g(y)^x.$$
If $g(y)$ equals zero, this formula is not valid, since the logarithm is not
defined there. Although strictly speaking $0^x$ is not differentiable at zero,
we return zero for convenience.
Note that up to [[patch--25]] this used to read
\begin{verbatim}
if F has ElementaryFunctionCategory then
dvpow2 l == log(first l) * first(l) ** second(l)
\end{verbatim}
which caused differentiating $0^x$ to fail. (Issue~\#19)
<<package COMBF CombinatorialFunction>>=
evaluate(opfact, iifact)$BasicOperatorFunctions1(F)
evaluate(oppow, iipow)
evaluate(opperm, iiperm)
evaluate(opbinom, iibinom)
evaluate(opsum, isum)
evaluate(opdsum, iidsum)
evaluate(opprod, iprod)
evaluate(opdprod, iidprod)
derivative(oppow, [dvpow1, dvpow2])
setProperty(opsum, SPECIALDIFF, dvsum@((List F, SE) -> F) pretend None)
setProperty(opdsum, SPECIALDIFF, dvdsum@((List F, SE)->F) pretend None)
setProperty(opprod, SPECIALDIFF, dvprod@((List F, SE)->F) pretend None)
setProperty(opdprod, SPECIALDIFF, dvdprod@((List F, SE)->F) pretend None)
@
The last four properties define special differentiation rules for sums and
products. Note that up to [[patch--45]] the rules for products were missing.
Thus products were differentiated according the usual chain-rule, which gave
incorrect results.
<<package COMBF CombinatorialFunction>>=
setProperty(opsum, SPECIALDISP, dsum@(List F -> O) pretend None)
setProperty(opdsum, SPECIALDISP, ddsum@(List F -> O) pretend None)
setProperty(opprod, SPECIALDISP, dprod@(List F -> O) pretend None)
setProperty(opdprod, SPECIALDISP, ddprod@(List F -> O) pretend None)
setProperty(opsum, SPECIALEQUAL, equalsumprod@((K,K) -> Boolean) pretend None)
setProperty(opdsum, SPECIALEQUAL, equaldsumprod@((K,K) -> Boolean) pretend None)
setProperty(opprod, SPECIALEQUAL, equalsumprod@((K,K) -> Boolean) pretend None)
setProperty(opdprod, SPECIALEQUAL, equaldsumprod@((K,K) -> Boolean) pretend None)
@
Finally, we set the properties for displaying sums and products and testing for
equality.
\section{package FSPECF FunctionalSpecialFunction}
<<package FSPECF FunctionalSpecialFunction>>=
)abbrev package FSPECF FunctionalSpecialFunction
++ Provides the special functions
++ Author: Manuel Bronstein
++ Date Created: 18 Apr 1989
++ Date Last Updated: 4 October 1993
++ Description: Provides some special functions over an integral domain.
++ Keywords: special, function.
FunctionalSpecialFunction(R, F): Exports == Implementation where
R: IntegralDomain
F: FunctionSpace R
OP ==> BasicOperator
K ==> Kernel F
SE ==> Symbol
Exports ==> with
belong? : OP -> Boolean
++ belong?(op) is true if op is a special function operator;
operator: OP -> OP
++ operator(op) returns a copy of op with the domain-dependent
++ properties appropriate for F;
++ error if op is not a special function operator
abs : F -> F
++ abs(f) returns the absolute value operator applied to f
Gamma : F -> F
++ Gamma(f) returns the formal Gamma function applied to f
Gamma : (F,F) -> F
++ Gamma(a,x) returns the incomplete Gamma function applied to a and x
Beta: (F,F) -> F
++ Beta(x,y) returns the beta function applied to x and y
digamma: F->F
++ digamma(x) returns the digamma function applied to x
polygamma: (F,F) ->F
++ polygamma(x,y) returns the polygamma function applied to x and y
besselJ: (F,F) -> F
++ besselJ(x,y) returns the besselj function applied to x and y
besselY: (F,F) -> F
++ besselY(x,y) returns the bessely function applied to x and y
besselI: (F,F) -> F
++ besselI(x,y) returns the besseli function applied to x and y
besselK: (F,F) -> F
++ besselK(x,y) returns the besselk function applied to x and y
airyAi: F -> F
++ airyAi(x) returns the airyai function applied to x
airyBi: F -> F
++ airyBi(x) returns the airybi function applied to x
iiGamma : F -> F
++ iiGamma(x) should be local but conditional;
iiabs : F -> F
++ iiabs(x) should be local but conditional;
Implementation ==> add
iabs : F -> F
iGamma: F -> F
opabs := operator('abs)$CommonOperators
opGamma := operator('Gamma)$CommonOperators
opGamma2 := operator('Gamma2)$CommonOperators
opBeta := operator('Beta)$CommonOperators
opdigamma := operator('digamma)$CommonOperators
oppolygamma := operator('polygamma)$CommonOperators
opBesselJ := operator('besselJ)$CommonOperators
opBesselY := operator('besselY)$CommonOperators
opBesselI := operator('besselI)$CommonOperators
opBesselK := operator('besselK)$CommonOperators
opAiryAi := operator('airyAi)$CommonOperators
opAiryBi := operator('airyBi)$CommonOperators
abs x == opabs x
Gamma(x) == opGamma(x)
Gamma(a,x) == opGamma2(a,x)
Beta(x,y) == opBeta(x,y)
digamma x == opdigamma(x)
polygamma(k,x)== oppolygamma(k,x)
besselJ(a,x) == opBesselJ(a,x)
besselY(a,x) == opBesselY(a,x)
besselI(a,x) == opBesselI(a,x)
besselK(a,x) == opBesselK(a,x)
airyAi(x) == opAiryAi(x)
airyBi(x) == opAiryBi(x)
belong? op == has?(op, 'special)
operator op ==
is?(op,'abs) => opabs
is?(op,'Gamma) => opGamma
is?(op,'Gamma2) => opGamma2
is?(op,'Beta) => opBeta
is?(op,'digamma) => opdigamma
is?(op,'polygamma)=> oppolygamma
is?(op,'besselJ) => opBesselJ
is?(op,'besselY) => opBesselY
is?(op,'besselI) => opBesselI
is?(op,'besselK) => opBesselK
is?(op,'airyAi) => opAiryAi
is?(op,'airyBi) => opAiryBi
error "Not a special operator"
-- Could put more unconditional special rules for other functions here
iGamma x ==
one? x => x
kernel(opGamma, x)
iabs x ==
zero? x => 0
is?(x, opabs) => x
before?(x,0) => kernel(opabs, -x)
kernel(opabs, x)
-- Could put more conditional special rules for other functions here
if R has abs : R -> R then
iiabs x ==
(r := retractIfCan(x)@Union(Fraction Polynomial R, "failed"))
case "failed" => iabs x
f := r::Fraction Polynomial R
(a := retractIfCan(numer f)@Union(R, "failed")) case "failed" or
(b := retractIfCan(denom f)@Union(R,"failed")) case "failed" => iabs x
abs(a::R)::F / abs(b::R)::F
else iiabs x == iabs x
if R has SpecialFunctionCategory then
iiGamma x ==
(r:=retractIfCan(x)@Union(R,"failed")) case "failed" => iGamma x
Gamma(r::R)::F
else
if R has RetractableTo Integer then
iiGamma x ==
(r := retractIfCan(x)@Union(Integer, "failed")) case Integer
and (r::Integer >= 1) => factorial(r::Integer - 1)::F
iGamma x
else
iiGamma x == iGamma x
-- Default behaviour is to build a kernel
evaluate(opGamma, iiGamma)$BasicOperatorFunctions1(F)
evaluate(opabs, iiabs)$BasicOperatorFunctions1(F)
import Fraction Integer
ahalf: F := recip(2::F)::F
athird: F := recip(2::F)::F
twothirds: F := 2*recip(3::F)::F
lzero(l: List F): F == 0
iBesselJGrad(l: List F): F ==
n := first l; x := second l
ahalf * (besselJ (n-1,x) - besselJ (n+1,x))
iBesselYGrad(l: List F): F ==
n := first l; x := second l
ahalf * (besselY (n-1,x) - besselY (n+1,x))
iBesselIGrad(l: List F): F ==
n := first l; x := second l
ahalf * (besselI (n-1,x) + besselI (n+1,x))
iBesselKGrad(l: List F): F ==
n := first l; x := second l
ahalf * (-besselK (n-1,x) - besselK (n+1,x))
ipolygammaGrad(l: List F): F ==
n := first l; x := second l
polygamma(n+1, x)
iBetaGrad1(l: List F): F ==
x := first l; y := second l
Beta(x,y)*(digamma x - digamma(x+y))
iBetaGrad2(l: List F): F ==
x := first l; y := second l
Beta(x,y)*(digamma y - digamma(x+y))
if F has ElementaryFunctionCategory then
iGamma2Grad(l: List F):F ==
a := first l; x := second l
- x ** (a - 1) * exp(-x)
derivative(opGamma2, [lzero, iGamma2Grad])
derivative(opabs, abs(#1) * inv(#1))
derivative(opGamma, digamma #1 * Gamma #1)
derivative(opBeta, [iBetaGrad1, iBetaGrad2])
derivative(opdigamma, polygamma(1, #1))
derivative(oppolygamma, [lzero, ipolygammaGrad])
derivative(opBesselJ, [lzero, iBesselJGrad])
derivative(opBesselY, [lzero, iBesselYGrad])
derivative(opBesselI, [lzero, iBesselIGrad])
derivative(opBesselK, [lzero, iBesselKGrad])
@
\section{package SUMFS FunctionSpaceSum}
<<package SUMFS FunctionSpaceSum>>=
)abbrev package SUMFS FunctionSpaceSum
++ Top-level sum function
++ Author: Manuel Bronstein
++ Date Created: ???
++ Date Last Updated: 19 April 1991
++ Description: computes sums of top-level expressions;
FunctionSpaceSum(R, F): Exports == Implementation where
R: Join(IntegralDomain,
RetractableTo Integer, LinearlyExplicitRingOver Integer)
F: Join(FunctionSpace R, CombinatorialOpsCategory,
AlgebraicallyClosedField, TranscendentalFunctionCategory)
SE ==> Symbol
K ==> Kernel F
Exports ==> with
sum: (F, SE) -> F
++ sum(a(n), n) returns A(n) such that A(n+1) - A(n) = a(n);
sum: (F, SegmentBinding F) -> F
++ sum(f(n), n = a..b) returns f(a) + f(a+1) + ... + f(b);
Implementation ==> add
import ElementaryFunctionStructurePackage(R, F)
import GosperSummationMethod(IndexedExponents K, K, R,
SparseMultivariatePolynomial(R, K), F)
innersum: (F, K) -> Union(F, "failed")
notRF? : (F, K) -> Boolean
newk : () -> K
newk() == kernel(new()$SE)
sum(x:F, s:SegmentBinding F) ==
k := kernel(variable s)@K
(u := innersum(x, k)) case "failed" => summation(x, s)
eval(u::F, k, 1 + hi segment s) - eval(u::F, k, lo segment s)
sum(x:F, v:SE) ==
(u := innersum(x, kernel(v)@K)) case "failed" => summation(x,v)
u::F
notRF?(f, k) ==
for kk in tower f repeat
member?(k, tower(kk::F)) and (symbolIfCan(kk) case "failed") =>
return true
false
innersum(x, k) ==
zero? x => 0
notRF?(f := normalize(x / (x1 := eval(x, k, k::F - 1))), k) =>
"failed"
(u := GospersMethod(f, k, newk)) case "failed" => "failed"
x1 * eval(u::F, k, k::F - 1)
@
\section{License}
<<license>>=
--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
--All rights reserved.
--Copyright (C) 2007-2009, Gabriel Dos Reis.
--All rights reserved.
--
--Redistribution and use in source and binary forms, with or without
--modification, are permitted provided that the following conditions are
--met:
--
-- - Redistributions of source code must retain the above copyright
-- notice, this list of conditions and the following disclaimer.
--
-- - Redistributions in binary form must reproduce the above copyright
-- notice, this list of conditions and the following disclaimer in
-- the documentation and/or other materials provided with the
-- distribution.
--
-- - Neither the name of The Numerical ALgorithms Group Ltd. nor the
-- names of its contributors may be used to endorse or promote products
-- derived from this software without specific prior written permission.
--
--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
@
<<*>>=
<<license>>
-- SPAD files for the functional world should be compiled in the
-- following order:
--
-- op kl function funcpkgs manip algfunc
-- elemntry constant funceval COMBFUNC fe
<<category COMBOPC CombinatorialOpsCategory>>
<<package COMBF CombinatorialFunction>>
<<package FSPECF FunctionalSpecialFunction>>
<<package SUMFS FunctionSpaceSum>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}