build open-axiom
--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
)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;
)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
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]
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]
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
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)
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))
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]
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)
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)
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)
)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])
)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)