path: root/src/algebra/combfunc.spad.pamphlet

\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}