\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} <>= )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. <>= )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) <>= 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. <>= 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} <>= 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. <>= 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) <>= 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. <>= 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. <>= 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. <>= 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) <>= 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. <>= 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} <>= )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} <>= )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} <>= --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. @ <<*>>= <> -- SPAD files for the functional world should be compiled in the -- following order: -- -- op kl function funcpkgs manip algfunc -- elemntry constant funceval COMBFUNC fe <> <> <> <> @ \eject \begin{thebibliography}{99} \bibitem{1} nothing \end{thebibliography} \end{document}