Initial population.

1 files changed, 913 insertions, 0 deletions

diff --git a/src/algebra/combfunc.spad.pamphlet b/src/algebra/combfunc.spad.pamphlet
new file mode 100644
index 00000000..ffbf5887
--- /dev/null
+++ b/src/algebra/combfunc.spad.pamphlet
@@ -0,0 +1,913 @@
+\documentclass{article}
+\usepackage{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: Join(OrderedSet, IntegralDomain)
+  F: FunctionSpace R
+
+  OP  ==> BasicOperator
+  K   ==> Kernel F
+  SE  ==> Symbol
+  O   ==> OutputForm
+  SMP ==> SparseMultivariatePolynomial(R, K)
+  Z   ==> Integer
+
+  POWER        ==> "%power"::Symbol
+  OPEXP        ==> "exp"::Symbol
+  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"::Symbol)$CommonOperators
+    opperm  := operator("permutation"::Symbol)$CommonOperators
+    opbinom := operator("binomial"::Symbol)$CommonOperators
+    opsum   := operator("summation"::Symbol)$CommonOperators
+    opdsum  := operator("%defsum"::Symbol)$CommonOperators
+    opprod  := operator("product"::Symbol)$CommonOperators
+    opdprod := operator("%defprod"::Symbol)$CommonOperators
+    oppow   := operator(POWER::Symbol)$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) ==
+      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"::Symbol)   => opfact
+      is?(op, "permutation"::Symbol) => opperm
+      is?(op, "binomial"::Symbol)    => opbinom
+      is?(op, "summation"::Symbol)   => opsum
+      is?(op, "%defsum"::Symbol)     => opdsum
+      is?(op, "product"::Symbol)     => opprod
+      is?(op, "%defprod"::Symbol)    => opdprod
+      is?(op, POWER)                 => oppow
+      error "Not a combinatorial operator"
+
+    iprod l ==
+      zero? first l => 0
+--      one? first l => 1
+      (first l = 1) => 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
+      zero? x or (x = 1) => 1
+      kernel(opfact, x)
+
+    ibinom l ==
+      n := first l
+      ((p := second l) = 0) or (p = n) => 1
+--      one? p or (p = n - 1) => n
+      (p = 1) 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
+        (x = 1) or zero?(n: F := second l) => 1
+--        one? n => x
+        (n = 1) => 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 =>
+        ((y := first argument(rec.var))=1) 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 s>0=>
+              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: Join(OrderedSet, 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"::Symbol)$CommonOperators
+    opGamma     := operator("Gamma"::Symbol)$CommonOperators
+    opGamma2    := operator("Gamma2"::Symbol)$CommonOperators
+    opBeta      := operator("Beta"::Symbol)$CommonOperators
+    opdigamma   := operator("digamma"::Symbol)$CommonOperators
+    oppolygamma := operator("polygamma"::Symbol)$CommonOperators
+    opBesselJ   := operator("besselJ"::Symbol)$CommonOperators
+    opBesselY   := operator("besselY"::Symbol)$CommonOperators
+    opBesselI   := operator("besselI"::Symbol)$CommonOperators
+    opBesselK   := operator("besselK"::Symbol)$CommonOperators
+    opAiryAi    := operator("airyAi"::Symbol)$CommonOperators
+    opAiryBi    := operator("airyBi"::Symbol)$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"::Symbol)      => opabs
+      is?(op, "Gamma"::Symbol)    => opGamma
+      is?(op, "Gamma2"::Symbol)   => opGamma2
+      is?(op, "Beta"::Symbol)     => opBeta
+      is?(op, "digamma"::Symbol)  => opdigamma
+      is?(op, "polygamma"::Symbol)=> oppolygamma
+      is?(op, "besselJ"::Symbol)  => opBesselJ
+      is?(op, "besselY"::Symbol)  => opBesselY
+      is?(op, "besselI"::Symbol)  => opBesselI
+      is?(op, "besselK"::Symbol)  => opBesselK
+      is?(op, "airyAi"::Symbol)   => opAiryAi
+      is?(op, "airyBi"::Symbol)   => opAiryBi
+
+      error "Not a special operator"
+
+    -- Could put more unconditional special rules for other functions here
+    iGamma x ==
+--      one? x => x
+      (x = 1) => x
+      kernel(opGamma, x)
+
+    iabs x ==
+      zero? x => 0
+      is?(x, opabs) => x
+      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, OrderedSet,
+          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.
+--
+--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}