Edit detail for SandBoxCombfunc revision 5 of 5

Add conjugate

(1) -> )library BOP BOP1

   BasicOperator is now explicitly exposed in frame initial 
   BasicOperator will be automatically loaded when needed from 
      /var/aw/var/LatexWiki/BOP.NRLIB/BOP
   BasicOperatorFunctions1 is now explicitly exposed in frame initial 
   BasicOperatorFunctions1 will be automatically loaded when needed 
      from /var/aw/var/LatexWiki/BOP1.NRLIB/BOP1

)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 : Join(Comparable, 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

    COMB := 'comb

    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

    -- Must be a macro to produce fresh symbol for each use
    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::Symbol)$CommonOperators

    factorial x          == opfact x
    binomial(x, y)       == opbinom [x, y]
    permutation(x, y)    == opperm [x, y]

    import from F
    import from 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(x +-> K2fact(x, l), x +-> x::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)   => 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
      (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

         -- Try to explicitly evaluete binomial coefficient
         -- We currently simplify binomial coefficients for non-negative
         -- integral second argument, using the formula
         -- $$ \binom{n}{k}=\frac{1}{k!}\prod_{i=0..k-1} (n-i), $$
         iibinom l ==
           (s := retractIfCan(second l)@Union(R,"failed")) case R and
              (t := retractIfCan(s)@Union(Z,"failed")) case Z and t>0 =>
                ans := 1::F
                for i in 0..t-1 repeat
                    ans := ans*(first l - i::R::F)
                (1/factorial t) * ans
           (s := retractIfCan(first l-second l)@Union(R,"failed")) case R and
             (t := retractIfCan(s)@Union(Z,"failed")) case Z and t>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)

    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
)boot $tryRecompileArguments := nil
++ Provides the special functions
++ Author: Manuel Bronstein
++ Date Created: 18 Apr 1989
++ Description: Provides some special functions over an integral domain.
++ Keywords: special, function.
++ Date Last Updated: 2 October 2014
++   Added: 'conjugate'
++   Modied: 'abs'
FunctionalSpecialFunction(R, F) : Exports == Implementation where
  R : Join(Comparable, IntegralDomain)
  F : FunctionSpace R

  OP  ==> BasicOperator
  K   ==> Kernel F
  SE  ==> Symbol
  SPECIALDIFF  ==> '%specialDiff

  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
    conjugate: F -> F
      ++ conjugate(f) returns the conjugate 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 Airy Ai function applied to x
    airyAiPrime :    F -> F
      ++ airyAiPrime(x) returns the derivative of Airy Ai function applied to x
    airyBi :    F -> F
      ++ airyBi(x) returns the Airy Bi function applied to x
    airyBiPrime :    F -> F
      ++ airyBiPrime(x) returns the derivative of Airy Bi function applied to x
    lambertW :  F -> F
      ++ lambertW(x) is the Lambert W function at x
    polylog : (F, F) -> F
      ++ polylog(s, x) is the polylogarithm of order s at x
    weierstrassP : (F, F, F) -> F
      ++ weierstrassP(g2, g3, x)
    weierstrassPPrime : (F, F, F) -> F
      ++ weierstrassPPrime(g2, g3, x)
    weierstrassSigma : (F, F, F) -> F
      ++ weierstrassSigma(g2, g3, x)
    weierstrassZeta : (F, F, F) -> F
      ++ weierstrassZeta(g2, g3, x)
    -- weierstrassPInverse : (F, F, F) -> F
    --    ++ weierstrassPInverse(g2, g3, z) is the inverse of Weierstrass
    --    ++ P function, defined by the formula
    --    ++ \spad{weierstrassP(g2, g3, weierstrassPInverse(g2, g3, z)) = z}
    whittakerM : (F, F, F) -> F
        ++ whittakerM(k, m, z) is the Whittaker M function
    whittakerW : (F, F, F) -> F
        ++ whittakerW(k, m, z) is the Whittaker W function
    angerJ : (F, F) -> F
        ++ angerJ(v, z) is the Anger J function
    weberE : (F, F) -> F
        ++ weberE(v, z) is the Weber E function
    struveH : (F, F) -> F
        ++ struveH(v, z) is the Struve H function
    struveL : (F, F) -> F
        ++ struveL(v, z) is the Struve L function defined by the formula
        ++ \spad{struveL(v, z) = -%i^exp(-v*%pi*%i/2)*struveH(v, %i*z)}
    hankelH1 : (F, F) -> F
        ++ hankelH1(v, z) is first Hankel function (Bessel function of
        ++ the third kind)
    hankelH2 : (F, F) -> F
        ++ hankelH2(v, z) is the second Hankel function (Bessel function of
        ++ the third kind)
    lommelS1 : (F, F, F) -> F
        ++ lommelS1(mu, nu, z) is the Lommel s function
    lommelS2 : (F, F, F) -> F
        ++ lommelS2(mu, nu, z) is the Lommel S function
    kummerM : (F, F, F) -> F
        ++ kummerM(a, b, z) is the Kummer M function
    kummerU : (F, F, F) -> F
        ++ kummerU(a, b, z) is the Kummer U function
    legendreP : (F, F, F) -> F
        ++ legendreP(nu, mu, z) is the Legendre P function
    legendreQ : (F, F, F) -> F
        ++ legendreQ(nu, mu, z) is the Legendre Q function
    kelvinBei : (F, F) -> F
        ++ kelvinBei(v, z) is the Kelvin bei function defined by equality
        ++ \spad{kelvinBei(v, z) = imag(besselJ(v, exp(3*%pi*%i/4)*z))}
        ++ for z and v real
    kelvinBer : (F, F) -> F
        ++ kelvinBer(v, z) is the Kelvin ber function defined by equality
        ++ \spad{kelvinBer(v, z) = real(besselJ(v, exp(3*%pi*%i/4)*z))}
        ++ for z and v real
    kelvinKei : (F, F) -> F
        ++ kelvinKei(v, z) is the Kelvin kei function defined by equality
        ++ \spad{kelvinKei(v, z) =
        ++ imag(exp(-v*%pi*%i/2)*besselK(v, exp(%pi*%i/4)*z))}
        ++ for z and v real
    kelvinKer : (F, F) -> F
        ++ kelvinKer(v, z) is the Kelvin kei function defined by equality
        ++ \spad{kelvinKer(v, z) =
        ++ real(exp(-v*%pi*%i/2)*besselK(v, exp(%pi*%i/4)*z))}
        ++ for z and v real
    ellipticK : F -> F
        ++ ellipticK(m) is the complete elliptic integral of the
        ++ first kind: \spad{ellipticK(m) =
        ++ integrate(1/sqrt((1-t^2)*(1-m*t^2)), t = 0..1)}
    ellipticE : F -> F
        ++ ellipticE(m) is the complete elliptic integral of the
        ++ second kind: \spad{ellipticE(m) =
        ++ integrate(sqrt(1-m*t^2)/sqrt(1-t^2), t = 0..1)}
    ellipticE : (F, F) -> F
        ++ ellipticE(z, m) is the incomplete elliptic integral of the
        ++ second kind: \spad{ellipticE(z, m) =
        ++ integrate(sqrt(1-m*t^2)/sqrt(1-t^2), t = 0..z)}
    ellipticF : (F, F) -> F
        ++ ellipticF(z, m) is the incomplete elliptic integral of the
        ++ first kind : \spad{ellipticF(z, m) =
        ++ integrate(1/sqrt((1-t^2)*(1-m*t^2)), t = 0..z)}
    ellipticPi : (F, F, F) -> F
        ++ ellipticPi(z, n, m) is the incomplete elliptic integral of
        ++ the third kind: \spad{ellipticPi(z, n, m) =
        ++ integrate(1/((1-n*t^2)*sqrt((1-t^2)*(1-m*t^2))), t = 0..z)}
    jacobiSn : (F, F) -> F
        ++ jacobiSn(z, m) is the Jacobi elliptic sn function, defined
        ++ by the formula \spad{jacobiSn(ellipticF(z, m), m) = z}
    jacobiCn : (F, F) -> F
        ++ jacobiCn(z, m) is the Jacobi elliptic cn function, defined
        ++ by \spad{jacobiCn(z, m)^2 + jacobiSn(z, m)^2 = 1} and
        ++ \spad{jacobiCn(0, m) = 1}
    jacobiDn : (F, F) -> F
        ++ jacobiDn(z, m) is the Jacobi elliptic dn function, defined
        ++ by \spad{jacobiDn(z, m)^2 + m*jacobiSn(z, m)^2 = 1} and
        ++ \spad{jacobiDn(0, m) = 1}
    jacobiZeta : (F, F) -> F
        ++ jacobiZeta(z, m) is the Jacobi elliptic zeta function, defined
        ++ by \spad{D(jacobiZeta(z, m), z) =
        ++ jacobiDn(z, m)^2 - ellipticE(m)/ellipticK(m)} and
        ++ \spad{jacobiZeta(0, m) = 0}.
    jacobiTheta : (F, F) -> F
        ++ jacobiTheta(q, z) is the third Jacobi Theta function
    lerchPhi : (F, F, F) -> F
        ++ lerchPhi(z, s, a) is the Lerch Phi function
    riemannZeta : F -> F
        ++ riemannZeta(z) is the Riemann Zeta function
    charlierC : (F, F, F) -> F
        ++ charlierC(n, a, z) is the Charlier polynomial
    hermiteH : (F, F) -> F
        ++ hermiteH(n, z) is the Hermite polynomial
    jacobiP : (F, F, F, F) -> F
        ++ jacobiP(n, a, b, z) is the Jacobi polynomial
    laguerreL: (F, F, F) -> F
        ++ laguerreL(n, a, z) is the Laguerre polynomial
    meixnerM : (F, F, F, F) -> F
        ++ meixnerM(n, b, c, z) is the Meixner polynomial
    if F has RetractableTo(Integer) then
        hypergeometricF : (List F, List F, F) -> F
          ++ hypergeometricF(la, lb, z) is the generalized hypergeometric
          ++ function
        meijerG : (List F, List F, List F, List F, F) -> F
          ++ meijerG(la, lb, lc, ld, z) is the meijerG function
    -- Functions below should be local but conditional
    iiGamma : F -> F
      ++ iiGamma(x) should be local but conditional;
    iiabs     : F -> F
      ++ iiabs(x) should be local but conditional;
    iiconjugate: F -> F
      ++ iiconjugate(x) should be local but conditional;
    iiBeta     : List F -> F
      ++ iiBeta(x) should be local but conditional;
    iidigamma  : F -> F
      ++ iidigamma(x) should be local but conditional;
    iipolygamma : List F -> F
      ++ iipolygamma(x) should be local but conditional;
    iiBesselJ  : List F -> F
      ++ iiBesselJ(x) should be local but conditional;
    iiBesselY  : List F -> F
      ++ iiBesselY(x) should be local but conditional;
    iiBesselI  : List F -> F
      ++ iiBesselI(x) should be local but conditional;
    iiBesselK  : List F -> F
      ++ iiBesselK(x) should be local but conditional;
    iiAiryAi   : F -> F
      ++ iiAiryAi(x) should be local but conditional;
    iiAiryAiPrime : F -> F
      ++ iiAiryAiPrime(x) should be local but conditional;
    iiAiryBi   : F -> F
      ++ iiAiryBi(x) should be local but conditional;
    iiAiryBiPrime : F -> F
      ++ iiAiryBiPrime(x) should be local but conditional;
    iAiryAi   : F -> F
      ++ iAiryAi(x) should be local but conditional;
    iAiryAiPrime : F -> F
      ++ iAiryAiPrime(x) should be local but conditional;
    iAiryBi   : F -> F
      ++ iAiryBi(x) should be local but conditional;
    iAiryBiPrime : F -> F
      ++ iAiryBiPrime(x) should be local but conditional;
    iiHypergeometricF : List F -> F
      ++ iiHypergeometricF(l) should be local but conditional;
    iiPolylog : (F, F) -> F
      ++ iiPolylog(x, s) should be local but conditional;
    iLambertW : F -> F
      ++ iLambertW(x) should be local but conditional;

  Implementation ==> add

    SPECIAL := 'special

    INP ==> InputForm

    SPECIALINPUT ==> '%specialInput

    iabs      : F -> F
    iGamma    : F -> F
    iBeta     : (F, F) -> F
    idigamma  : F -> F
    iiipolygamma : (F, F) -> F
    iiiBesselJ  : (F, F) -> F
    iiiBesselY  : (F, F) -> F
    iiiBesselI  : (F, F) -> F
    iiiBesselK  : (F, F) -> F
    iPolylog : List F -> F

    iWeierstrassP : (F, F, F) -> F
    iWeierstrassPPrime : (F, F, F) -> F
    iWeierstrassSigma : (F, F, F) -> F
    iWeierstrassZeta : (F, F, F) -> F

    iiWeierstrassP : List F -> F
    iiWeierstrassPPrime : List F -> F
    iiWeierstrassSigma : List F -> F
    iiWeierstrassZeta : List F -> F
    iiMeijerG : List F -> F

    opabs       := operator('abs)$CommonOperators
    opconjugate := operator('conjugate)$CommonOperators
    opsqrt      := operator('sqrt)$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
    opAiryAiPrime := operator('airyAiPrime)$CommonOperators
    opAiryBi    := operator('airyBi)$CommonOperators
    opAiryBiPrime := operator('airyBiPrime)$CommonOperators
    opLambertW := operator('lambertW)$CommonOperators
    opPolylog := operator('polylog)$CommonOperators
    opWeierstrassP := operator('weierstrassP)$CommonOperators
    opWeierstrassPPrime := operator('weierstrassPPrime)$CommonOperators
    opWeierstrassSigma := operator('weierstrassSigma)$CommonOperators
    opWeierstrassZeta := operator('weierstrassZeta)$CommonOperators
    opHypergeometricF := operator('hypergeometricF)$CommonOperators
    opMeijerG := operator('meijerG)$CommonOperators
    opCharlierC := operator('charlierC)$CommonOperators
    opHermiteH := operator('hermiteH)$CommonOperators
    opJacobiP := operator('jacobiP)$CommonOperators
    opLaguerreL := operator('laguerreL)$CommonOperators
    opMeixnerM := operator('meixnerM)$CommonOperators
    op_log_gamma := operator('%logGamma)$CommonOperators
    op_eis := operator('%eis)$CommonOperators
    op_erfs := operator('%erfs)$CommonOperators
    op_erfis := operator('%erfis)$CommonOperators

    abs x         == opabs x
    conjugate x   == opconjugate 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)
    airyAiPrime(x) == opAiryAiPrime(x)
    airyBi(x)     == opAiryBi(x)
    airyBiPrime(x) == opAiryBiPrime(x)
    lambertW(x) == opLambertW(x)
    polylog(s, x) == opPolylog(s, x)
    weierstrassP(g2, g3, x) == opWeierstrassP(g2, g3, x)
    weierstrassPPrime(g2, g3, x) == opWeierstrassPPrime(g2, g3, x)
    weierstrassSigma(g2, g3, x) == opWeierstrassSigma(g2, g3, x)
    weierstrassZeta(g2, g3, x) == opWeierstrassZeta(g2, g3, x)
    if F has RetractableTo(Integer) then
        hypergeometricF(a, b, z) ==
            nai := #a
            nbi := #b
            z = 0 and nai <= nbi + 1 => 1
            p := (#a)::F
            q := (#b)::F
            opHypergeometricF concat(concat(a, concat(b, [z])), [p, q])

        meijerG(a, b, c, d, z) ==
            n1 := (#a)::F
            n2 := (#b)::F
            m1 := (#c)::F
            m2 := (#d)::F
            opMeijerG concat(concat(a, concat(b,
                         concat(c, concat(d, [z])))), [n1, n2, m1, m2])

    import from List Kernel(F)
    opdiff := operator(operator('%diff)$CommonOperators)$F
    dummy ==> new()$SE :: F

    ahalf : F := recip(2::F)::F
    athird : F    := recip(3::F)::F
    afourth : F   := recip(4::F)::F
    asixth : F := recip(6::F)::F
    twothirds : F := 2*athird
    threehalfs : F := 3*ahalf

    -- Helpers for partially defined derivatives

    grad2(l : List F, t : SE, op : OP, d2 : (F, F) -> F ) : F ==
        x1 := l(1)
        x2 := l(2)
        dm := dummy
        differentiate(x1, t)*kernel(opdiff, [op [dm, x2], dm, x1])
             + differentiate(x2, t)*d2(x1, x2)

    grad3(l : List F, t : SE, op : OP, d3 : (F, F, F) -> F ) : F ==
        x1 := l(1)
        x2 := l(2)
        x3 := l(3)
        dm1 := dummy
        dm2 := dummy
        differentiate(x1, t)*kernel(opdiff, [op [dm1, x2, x3], dm1, x1])
           + differentiate(x2, t)*kernel(opdiff, [op [x1, dm2, x3], dm2, x2])
             + differentiate(x3, t)*d3(x1, x2, x3)

    grad4(l : List F, t : SE, op : OP, d4 : (F, F, F, F) -> F ) : F ==
        x1 := l(1)
        x2 := l(2)
        x3 := l(3)
        x4 := l(4)
        dm1 := dummy
        dm2 := dummy
        dm3 := dummy
        kd1 := kernel(opdiff, [op [dm1, x2, x3, x4], dm1, x1])
        kd2 := kernel(opdiff, [op [x1, dm2, x3, x4], dm2, x2])
        kd3 := kernel(opdiff, [op [x1, x2, dm3, x4], dm3, x3])
        differentiate(x1, t)*kd1 + differentiate(x2, t)*kd2 +
          differentiate(x3, t)*kd3 +
            differentiate(x4, t)*d4(x1, x2, x3, x4)

    -- handle WeierstrassPInverse

)if false

    opWeierstrassPInverse := operator('weierstrassPInverse)$CommonOperators

    weierstrassPInverse(g2, g3, z) == opWeierstrassPInverse(g2, g3, z)

    eWeierstrassPInverse(g2 : F, g3 : F, z : F) : F ==
        kernel(opWeierstrassPInverse, [g2, g3, z])

    elWeierstrassPInverse(l : List F) : F == eWeierstrassPInverse(l(1), l(2), l(3))
    evaluate(opWeierstrassPInverse, elWeierstrassPInverse)$BasicOperatorFunctions1(F)

    if F has RadicalCategory then

        eWeierstrassPInverseGrad_g2(l : List F) : F ==
            g2 := l(1)
            g3 := l(2)
            z := l(3)
            error "unimplemented"

        eWeierstrassPInverseGrad_g3(l : List F) : F ==
            g2 := l(1)
            g3 := l(2)
            z := l(3)
            error "unimplemented"

        eWeierstrassPInverseGrad_z(l : List F) : F ==
            g2 := l(1)
            g3 := l(2)
            z := l(3)
            1/sqrt(4*z^3 - g2*z - g3)

        derivative(opWeierstrassPInverse, [eWeierstrassPInverseGrad_g2,
               eWeierstrassPInverseGrad_g3, eWeierstrassPInverseGrad_z])

)endif

    -- handle WhittakerM

    opWhittakerM := operator('whittakerM)$CommonOperators

    whittakerM(k, m, z) == opWhittakerM(k, m, z)

    eWhittakerM(k : F, m : F, z : F) : F ==
        kernel(opWhittakerM, [k, m, z])

    elWhittakerM(l : List F) : F == eWhittakerM(l(1), l(2), l(3))
    evaluate(opWhittakerM, elWhittakerM)$BasicOperatorFunctions1(F)

    eWhittakerMGrad_z(k : F, m : F, z : F) : F ==
        (ahalf - k/z)*whittakerM(k, m, z) +
            (ahalf + k + m)*whittakerM(k + 1, m, z)/z

    dWhittakerM(l : List F, t : SE) : F ==
        grad3(l, t, opWhittakerM, eWhittakerMGrad_z)

    setProperty(opWhittakerM, SPECIALDIFF, dWhittakerM@((List F, SE)->F)
                                                 pretend None)

    -- handle WhittakerW

    opWhittakerW := operator('whittakerW)$CommonOperators

    whittakerW(k, m, z) == opWhittakerW(k, m, z)

    eWhittakerW(k : F, m : F, z : F) : F ==
        kernel(opWhittakerW, [k, m, z])

    elWhittakerW(l : List F) : F == eWhittakerW(l(1), l(2), l(3))
    evaluate(opWhittakerW, elWhittakerW)$BasicOperatorFunctions1(F)

    eWhittakerWGrad_z(k : F, m : F, z : F) : F ==
        (ahalf - k/z)*whittakerW(k, m, z) - whittakerW(k + 1, m, z)/z

    dWhittakerW(l : List F, t : SE) : F ==
        grad3(l, t, opWhittakerW, eWhittakerWGrad_z)

    setProperty(opWhittakerW, SPECIALDIFF, dWhittakerW@((List F, SE)->F)
                                                 pretend None)

    -- handle AngerJ

    opAngerJ := operator('angerJ)$CommonOperators

    angerJ(v, z) == opAngerJ(v, z)

    if F has TranscendentalFunctionCategory then

        eAngerJ(v : F, z : F) : F ==
            z = 0 => sin(v*pi())/(v*pi())
            kernel(opAngerJ, [v, z])

        elAngerJ(l : List F) : F == eAngerJ(l(1), l(2))
        evaluate(opAngerJ, elAngerJ)$BasicOperatorFunctions1(F)

        eAngerJGrad_z(v : F, z : F) : F ==
            -angerJ(v + 1, z) + v*angerJ(v, z)/z - sin(v*pi())/(pi()*z)

        dAngerJ(l : List F, t : SE) : F ==
            grad2(l, t, opAngerJ, eAngerJGrad_z)

        setProperty(opAngerJ, SPECIALDIFF, dAngerJ@((List F, SE)->F)
                                                 pretend None)

    else

        eeAngerJ(l : List F) : F == kernel(opAngerJ, l)
        evaluate(opAngerJ, eeAngerJ)$BasicOperatorFunctions1(F)

    -- handle WeberE

    opWeberE := operator('weberE)$CommonOperators

    weberE(v, z) == opWeberE(v, z)

    if F has TranscendentalFunctionCategory then

        eWeberE(v : F, z : F) : F ==
            z = 0 => 2*sin(ahalf*v*pi())^2/(v*pi())
            kernel(opWeberE, [v, z])

        elWeberE(l : List F) : F == eWeberE(l(1), l(2))
        evaluate(opWeberE, elWeberE)$BasicOperatorFunctions1(F)

        eWeberEGrad_z(v : F, z : F) : F ==
            -weberE(v + 1, z) + v*weberE(v, z)/z - (1 - cos(v*pi()))/(pi()*z)

        dWeberE(l : List F, t : SE) : F ==
            grad2(l, t, opWeberE, eWeberEGrad_z)

        setProperty(opWeberE, SPECIALDIFF, dWeberE@((List F, SE)->F)
                                                 pretend None)

    else

        eeWeberE(l :  List F) : F == kernel(opWeberE, l)
        evaluate(opWeberE, eeWeberE)$BasicOperatorFunctions1(F)

    -- handle StruveH

    opStruveH := operator('struveH)$CommonOperators

    struveH(v, z) == opStruveH(v, z)

    eStruveH(v : F, z : F) : F ==
        kernel(opStruveH, [v, z])

    elStruveH(l : List F) : F == eStruveH(l(1), l(2))
    evaluate(opStruveH, elStruveH)$BasicOperatorFunctions1(F)

    if F has TranscendentalFunctionCategory
       and F has RadicalCategory then

        eStruveHGrad_z(v : F, z : F) : F ==
            -struveH(v + 1, z) + v*struveH(v, z)/z +
               (ahalf*z)^v/(sqrt(pi())*Gamma(v + threehalfs))

        dStruveH(l : List F, t : SE) : F ==
            grad2(l, t, opStruveH, eStruveHGrad_z)

        setProperty(opStruveH, SPECIALDIFF, dStruveH@((List F, SE)->F)
                                                 pretend None)

    -- handle StruveL

    opStruveL := operator('struveL)$CommonOperators

    struveL(v, z) == opStruveL(v, z)

    eStruveL(v : F, z : F) : F ==
        kernel(opStruveL, [v, z])

    elStruveL(l : List F) : F == eStruveL(l(1), l(2))
    evaluate(opStruveL, elStruveL)$BasicOperatorFunctions1(F)

    if F has TranscendentalFunctionCategory
       and F has RadicalCategory then

        eStruveLGrad_z(v : F, z : F) : F ==
            struveL(v + 1, z) + v*struveL(v, z)/z +
               (ahalf*z)^v/(sqrt(pi())*Gamma(v + threehalfs))

        dStruveL(l : List F, t : SE) : F ==
            grad2(l, t, opStruveL, eStruveLGrad_z)

        setProperty(opStruveL, SPECIALDIFF, dStruveL@((List F, SE)->F)
                                                 pretend None)

    -- handle HankelH1

    opHankelH1 := operator('hankelH1)$CommonOperators

    hankelH1(v, z) == opHankelH1(v, z)

    eHankelH1(v : F, z : F) : F ==
        kernel(opHankelH1, [v, z])

    elHankelH1(l : List F) : F == eHankelH1(l(1), l(2))
    evaluate(opHankelH1, elHankelH1)$BasicOperatorFunctions1(F)

    eHankelH1Grad_z(v : F, z : F) : F ==
        -hankelH1(v + 1, z) + v*hankelH1(v, z)/z

    dHankelH1(l : List F, t : SE) : F ==
        grad2(l, t, opHankelH1, eHankelH1Grad_z)

    setProperty(opHankelH1, SPECIALDIFF, dHankelH1@((List F, SE)->F)
                                                 pretend None)

    -- handle HankelH2

    opHankelH2 := operator('hankelH2)$CommonOperators

    hankelH2(v, z) == opHankelH2(v, z)

    eHankelH2(v : F, z : F) : F ==
        kernel(opHankelH2, [v, z])

    elHankelH2(l : List F) : F == eHankelH2(l(1), l(2))
    evaluate(opHankelH2, elHankelH2)$BasicOperatorFunctions1(F)

    eHankelH2Grad_z(v : F, z : F) : F ==
        -hankelH2(v + 1, z) + v*hankelH2(v, z)/z

    dHankelH2(l : List F, t : SE) : F ==
        grad2(l, t, opHankelH2, eHankelH2Grad_z)

    setProperty(opHankelH2, SPECIALDIFF, dHankelH2@((List F, SE)->F)
                                                 pretend None)

    -- handle LommelS1

    opLommelS1 := operator('lommelS1)$CommonOperators

    lommelS1(m, v, z) == opLommelS1(m, v, z)

    eLommelS1(m : F, v : F, z : F) : F ==
        kernel(opLommelS1, [m, v, z])

    elLommelS1(l : List F) : F == eLommelS1(l(1), l(2), l(3))
    evaluate(opLommelS1, elLommelS1)$BasicOperatorFunctions1(F)

    eLommelS1Grad_z(m : F, v : F, z : F) : F ==
        -v*lommelS1(m, v, z)/z + (m + v - 1)*lommelS1(m - 1, v - 1, z)

    dLommelS1(l : List F, t : SE) : F ==
        grad3(l, t, opLommelS1, eLommelS1Grad_z)

    setProperty(opLommelS1, SPECIALDIFF, dLommelS1@((List F, SE)->F)
                                                 pretend None)

    -- handle LommelS2

    opLommelS2 := operator('lommelS2)$CommonOperators

    lommelS2(mu, nu, z) == opLommelS2(mu, nu, z)

    eLommelS2(mu : F, nu : F, z : F) : F ==
        kernel(opLommelS2, [mu, nu, z])

    elLommelS2(l : List F) : F == eLommelS2(l(1), l(2), l(3))
    evaluate(opLommelS2, elLommelS2)$BasicOperatorFunctions1(F)

    eLommelS2Grad_z(m : F, v : F, z : F) : F ==
        -v*lommelS2(m, v, z)/z + (m + v - 1)*lommelS2(m - 1, v - 1, z)

    dLommelS2(l : List F, t : SE) : F ==
        grad3(l, t, opLommelS2, eLommelS2Grad_z)

    setProperty(opLommelS2, SPECIALDIFF, dLommelS2@((List F, SE)->F)
                                                 pretend None)

    -- handle KummerM

    opKummerM := operator('kummerM)$CommonOperators

    kummerM(mu, nu, z) == opKummerM(mu, nu, z)

    eKummerM(a : F, b : F, z : F) : F ==
        z = 0 => 1
        kernel(opKummerM, [a, b, z])

    elKummerM(l : List F) : F == eKummerM(l(1), l(2), l(3))
    evaluate(opKummerM, elKummerM)$BasicOperatorFunctions1(F)

    eKummerMGrad_z(a : F, b : F, z : F) : F ==
        ((z + a - b)*kummerM(a, b, z)+(b - a)*kummerM(a - 1, b, z))/z

    dKummerM(l : List F, t : SE) : F ==
        grad3(l, t, opKummerM, eKummerMGrad_z)

    setProperty(opKummerM, SPECIALDIFF, dKummerM@((List F, SE)->F)
                                                 pretend None)

    -- handle KummerU

    opKummerU := operator('kummerU)$CommonOperators

    kummerU(a, b, z) == opKummerU(a, b, z)

    eKummerU(a : F, b : F, z : F) : F ==
        kernel(opKummerU, [a, b, z])

    elKummerU(l : List F) : F == eKummerU(l(1), l(2), l(3))
    evaluate(opKummerU, elKummerU)$BasicOperatorFunctions1(F)

    eKummerUGrad_z(a : F, b : F, z : F) : F ==
        ((z + a - b)*kummerU(a, b, z) - kummerU(a - 1, b, z))/z

    dKummerU(l : List F, t : SE) : F ==
        grad3(l, t, opKummerU, eKummerUGrad_z)

    setProperty(opKummerU, SPECIALDIFF, dKummerU@((List F, SE)->F)
                                                 pretend None)

    -- handle LegendreP

    opLegendreP := operator('legendreP)$CommonOperators

    legendreP(nu, mu, z) == opLegendreP(nu, mu, z)

    eLegendreP(nu : F, mu : F, z : F) : F ==
        kernel(opLegendreP, [nu, mu, z])

    elLegendreP(l : List F) : F == eLegendreP(l(1), l(2), l(3))
    evaluate(opLegendreP, elLegendreP)$BasicOperatorFunctions1(F)

    eLegendrePGrad_z(nu : F, mu : F, z : F) : F ==
        (nu - mu + 1)*legendreP(nu + 1, mu, z) -
          (nu + 1)*z*legendreP(nu, mu, z)

    dLegendreP(l : List F, t : SE) : F ==
        grad3(l, t, opLegendreP, eLegendrePGrad_z)

    setProperty(opLegendreP, SPECIALDIFF, dLegendreP@((List F, SE)->F)
                                                 pretend None)

    -- handle LegendreQ

    opLegendreQ := operator('legendreQ)$CommonOperators

    legendreQ(nu, mu, z) == opLegendreQ(nu, mu, z)

    eLegendreQ(nu : F, mu : F, z : F) : F ==
        kernel(opLegendreQ, [nu, mu, z])

    elLegendreQ(l : List F) : F == eLegendreQ(l(1), l(2), l(3))
    evaluate(opLegendreQ, elLegendreQ)$BasicOperatorFunctions1(F)

    eLegendreQGrad_z(nu : F, mu : F, z : F) : F ==
        (nu - mu + 1)*legendreQ(nu + 1, mu, z) -
          (nu + 1)*z*legendreQ(nu, mu, z)

    dLegendreQ(l : List F, t : SE) : F ==
        grad3(l, t, opLegendreQ, eLegendreQGrad_z)

    setProperty(opLegendreQ, SPECIALDIFF, dLegendreQ@((List F, SE)->F)
                                                 pretend None)

    -- handle KelvinBei

    opKelvinBei := operator('kelvinBei)$CommonOperators

    kelvinBei(v, z) == opKelvinBei(v, z)

    eKelvinBei(v : F, z : F) : F ==
        kernel(opKelvinBei, [v, z])

    elKelvinBei(l : List F) : F == eKelvinBei(l(1), l(2))
    evaluate(opKelvinBei, elKelvinBei)$BasicOperatorFunctions1(F)

    if F has RadicalCategory then

        eKelvinBeiGrad_z(v : F, z : F) : F ==
            ahalf*sqrt(2::F)*(kelvinBei(v + 1, z) - kelvinBer(v + 1, z)) +
              v*kelvinBei(v, z)/z

        dKelvinBei(l : List F, t : SE) : F ==
            grad2(l, t, opKelvinBei, eKelvinBeiGrad_z)

        setProperty(opKelvinBei, SPECIALDIFF, dKelvinBei@((List F, SE)->F)
                                                 pretend None)

    -- handle KelvinBer

    opKelvinBer := operator('kelvinBer)$CommonOperators

    kelvinBer(v, z) == opKelvinBer(v, z)

    eKelvinBer(v : F, z : F) : F ==
        kernel(opKelvinBer, [v, z])

    elKelvinBer(l : List F) : F == eKelvinBer(l(1), l(2))
    evaluate(opKelvinBer, elKelvinBer)$BasicOperatorFunctions1(F)

    if F has RadicalCategory then

        eKelvinBerGrad_z(v : F, z : F) : F ==
            ahalf*sqrt(2::F)*(kelvinBer(v + 1, z) + kelvinBei(v + 1, z)) +
              v*kelvinBer(v, z)/z

        dKelvinBer(l : List F, t : SE) : F ==
            grad2(l, t, opKelvinBer, eKelvinBerGrad_z)

        setProperty(opKelvinBer, SPECIALDIFF, dKelvinBer@((List F, SE)->F)
                                                 pretend None)

    -- handle KelvinKei

    opKelvinKei := operator('kelvinKei)$CommonOperators

    kelvinKei(v, z) == opKelvinKei(v, z)

    eKelvinKei(v : F, z : F) : F ==
        kernel(opKelvinKei, [v, z])

    elKelvinKei(l : List F) : F == eKelvinKei(l(1), l(2))
    evaluate(opKelvinKei, elKelvinKei)$BasicOperatorFunctions1(F)

    if F has RadicalCategory then

        eKelvinKeiGrad_z(v : F, z : F) : F ==
            ahalf*sqrt(2::F)*(kelvinKei(v + 1, z) - kelvinKer(v + 1, z)) +
              v*kelvinKei(v, z)/z

        dKelvinKei(l : List F, t : SE) : F ==
            grad2(l, t, opKelvinKei, eKelvinKeiGrad_z)

        setProperty(opKelvinKei, SPECIALDIFF, dKelvinKei@((List F, SE)->F)
                                                 pretend None)

    -- handle KelvinKer

    opKelvinKer := operator('kelvinKer)$CommonOperators

    kelvinKer(v, z) == opKelvinKer(v, z)

    eKelvinKer(v : F, z : F) : F ==
        kernel(opKelvinKer, [v, z])

    elKelvinKer(l : List F) : F == eKelvinKer(l(1), l(2))
    evaluate(opKelvinKer, elKelvinKer)$BasicOperatorFunctions1(F)

    if F has RadicalCategory then

        eKelvinKerGrad_z(v : F, z : F) : F ==
            ahalf*sqrt(2::F)*(kelvinKer(v + 1, z) + kelvinKei(v + 1, z)) +
              v*kelvinKer(v, z)/z

        dKelvinKer(l : List F, t : SE) : F ==
            grad2(l, t, opKelvinKer, eKelvinKerGrad_z)

        setProperty(opKelvinKer, SPECIALDIFF, dKelvinKer@((List F, SE)->F)
                                                 pretend None)

    -- handle EllipticK

    opEllipticK := operator('ellipticK)$CommonOperators

    ellipticK(m) == opEllipticK(m)

    eEllipticK(m : F) : F ==
        kernel(opEllipticK, [m])

    elEllipticK(l : List F) : F == eEllipticK(l(1))
    evaluate(opEllipticK, elEllipticK)$BasicOperatorFunctions1(F)

    dEllipticK(m : F) : F ==
        ahalf*(ellipticE(m) - (1 - m)*ellipticK(m))/(m*(1 - m))

    derivative(opEllipticK, dEllipticK)

    -- handle one argument EllipticE

    opEllipticE := operator('ellipticE)$CommonOperators

    ellipticE(m) == opEllipticE(m)

    eEllipticE(m : F) : F ==
        kernel(opEllipticE, [m])

    elEllipticE(l : List F) : F == eEllipticE(l(1))
    evaluate(opEllipticE, elEllipticE)$BasicOperatorFunctions1(F)

    dEllipticE(m : F) : F ==
        ahalf*(ellipticE(m) - ellipticK(m))/m

    derivative(opEllipticE, dEllipticE)

    -- handle two argument EllipticE

    opEllipticE2 := operator('ellipticE2)$CommonOperators

    ellipticE(z, m) == opEllipticE2(z, m)

    eEllipticE2(z : F, m : F) : F ==
        z = 0 => 0
        z = 1 => eEllipticE(m)
        kernel(opEllipticE2, [z, m])

    elEllipticE2(l : List F) : F == eEllipticE2(l(1), l(2))
    evaluate(opEllipticE2, elEllipticE2)$BasicOperatorFunctions1(F)

    if F has RadicalCategory then

        eEllipticE2Grad_z(l : List F) : F ==
            z := l(1)
            m := l(2)
            sqrt(1 - m*z^2)/sqrt(1 - z^2)

        eEllipticE2Grad_m(l : List F) : F ==
            z := l(1)
            m := l(2)
            ahalf*(ellipticE(z, m) - ellipticF(z, m))/m

        derivative(opEllipticE2, [eEllipticE2Grad_z, eEllipticE2Grad_m])

    inEllipticE2(li : List INP) : INP ==
        convert cons(convert('ellipticE), li)

    input(opEllipticE2, inEllipticE2@((List INP) -> INP))

    -- handle EllipticF

    opEllipticF := operator('ellipticF)$CommonOperators

    ellipticF(z, m) == opEllipticF(z, m)

    eEllipticF(z : F, m : F) : F ==
        z = 0 => 0
        z = 1 => ellipticK(m)
        kernel(opEllipticF, [z, m])

    elEllipticF(l : List F) : F == eEllipticF(l(1), l(2))
    evaluate(opEllipticF, elEllipticF)$BasicOperatorFunctions1(F)

    if F has RadicalCategory then

        eEllipticFGrad_z(l : List F) : F ==
            z := l(1)
            m := l(2)
            1/(sqrt(1 - m*z^2)*sqrt(1 - z^2))

        eEllipticFGrad_m(l : List F) : F ==
            z := l(1)
            m := l(2)
            ahalf*((ellipticE(z, m) - (1 - m)*ellipticF(z, m))/m -
              z*sqrt(1 - z^2)/sqrt(1 - m*z^2))/(1 - m)

        derivative(opEllipticF, [eEllipticFGrad_z, eEllipticFGrad_m])

    -- handle EllipticPi

    opEllipticPi := operator('ellipticPi)$CommonOperators

    ellipticPi(z, n, m) == opEllipticPi(z, n, m)

    eEllipticPi(z : F, n : F, m : F) : F ==
        z = 0 => 0
        kernel(opEllipticPi, [z, n, m])

    elEllipticPi(l : List F) : F == eEllipticPi(l(1), l(2), l(3))
    evaluate(opEllipticPi, elEllipticPi)$BasicOperatorFunctions1(F)

    if F has RadicalCategory then

        eEllipticPiGrad_z(l : List F) : F ==
            z := l(1)
            n := l(2)
            m := l(3)
            1/((1 - n*z^2)*sqrt(1 - m*z^2)*sqrt(1 - z^2))

        eEllipticPiGrad_n(l : List F) : F ==
            z := l(1)
            n := l(2)
            m := l(3)
            t1 := -(n^2 - m)*ellipticPi(z, n, m)/((n - 1)*(n - m)*n)
            t2 := ellipticF(z, m)/((n - 1)*n)
            t3 := -ellipticE(z, m)/((n - 1)*(n - m))
            t4 := n*z*sqrt(1 - m*z^2)*sqrt(1 - z^2)/
                   ((1 - n*z^2)*(n - 1)*(n - m))
            ahalf*(t1 + t2 + t3 + t4)

        eEllipticPiGrad_m(l : List F) : F ==
            z := l(1)
            n := l(2)
            m := l(3)
            t1 := m*z*sqrt(1 - z^2)/sqrt(1 - m*z^2)
            t2 := (-ellipticE(z, m) + t1)/(1 - m)
            ahalf*(ellipticPi(z, n, m) + t2)/(n - m)

        derivative(opEllipticPi, [eEllipticPiGrad_z, eEllipticPiGrad_n,
                                  eEllipticPiGrad_m])

    -- handle JacobiSn

    opJacobiSn := operator('jacobiSn)$CommonOperators

    jacobiSn(z, m) == opJacobiSn(z, m)

    eJacobiSn(z : F, m : F) : F ==
        z = 0 => 0
        if is?(z, opEllipticF) then
            args := argument(retract(z)@K)
            m = args(2) => return args(1)
        kernel(opJacobiSn, [z, m])

    elJacobiSn : List F -> F
    elJacobiSn(l : List F) : F == eJacobiSn(l(1), l(2))
    evaluate(opJacobiSn, elJacobiSn)$BasicOperatorFunctions1(F)

    jacobiGradHelper(z : F, m : F) : F ==
        (z - ellipticE(jacobiSn(z, m), m)/(1 - m))/m

    eJacobiSnGrad_z(l : List F) : F ==
        z := l(1)
        m := l(2)
        jacobiCn(z, m)*jacobiDn(z, m)

    eJacobiSnGrad_m(l : List F) : F ==
        z := l(1)
        m := l(2)
        ahalf*(eJacobiSnGrad_z(l)*jacobiGradHelper(z, m) +
           jacobiSn(z, m)*jacobiCn(z, m)^2/(1 - m))

    derivative(opJacobiSn, [eJacobiSnGrad_z, eJacobiSnGrad_m])

    -- handle JacobiCn

    opJacobiCn := operator('jacobiCn)$CommonOperators

    jacobiCn(z, m) == opJacobiCn(z, m)

    eJacobiCn(z : F, m : F) : F ==
        z = 0 => 1
        kernel(opJacobiCn, [z, m])

    elJacobiCn(l : List F) : F == eJacobiCn(l(1), l(2))
    evaluate(opJacobiCn, elJacobiCn)$BasicOperatorFunctions1(F)

    eJacobiCnGrad_z(l : List F) : F ==
        z := l(1)
        m := l(2)
        -jacobiSn(z, m)*jacobiDn(z, m)

    eJacobiCnGrad_m(l : List F) : F ==
        z := l(1)
        m := l(2)
        ahalf*(eJacobiCnGrad_z(l)*jacobiGradHelper(z, m) -
           jacobiSn(z, m)^2*jacobiCn(z, m)/(1 - m))

    derivative(opJacobiCn, [eJacobiCnGrad_z, eJacobiCnGrad_m])

    -- handle JacobiDn

    opJacobiDn := operator('jacobiDn)$CommonOperators

    jacobiDn(z, m) == opJacobiDn(z, m)

    eJacobiDn(z : F, m : F) : F ==
        z = 0 => 1
        kernel(opJacobiDn, [z, m])

    elJacobiDn(l : List F) : F == eJacobiDn(l(1), l(2))
    evaluate(opJacobiDn, elJacobiDn)$BasicOperatorFunctions1(F)

    eJacobiDnGrad_z(l : List F) : F ==
        z := l(1)
        m := l(2)
        -m*jacobiSn(z, m)*jacobiCn(z, m)

    eJacobiDnGrad_m(l : List F) : F ==
        z := l(1)
        m := l(2)
        ahalf*(eJacobiDnGrad_z(l)*jacobiGradHelper(z, m) -
           jacobiSn(z, m)^2*jacobiDn(z, m)/(1 - m))

    derivative(opJacobiDn, [eJacobiDnGrad_z, eJacobiDnGrad_m])

    -- handle JacobiZeta

    opJacobiZeta := operator('jacobiZeta)$CommonOperators

    jacobiZeta(z, m) == opJacobiZeta(z, m)

    eJacobiZeta(z : F, m : F) : F ==
        z = 0 => 0
        kernel(opJacobiZeta, [z, m])

    elJacobiZeta(l : List F) : F == eJacobiZeta(l(1), l(2))
    evaluate(opJacobiZeta, elJacobiZeta)$BasicOperatorFunctions1(F)

    eJacobiZetaGrad_z(l : List F) : F ==
        z := l(1)
        m := l(2)
        dn := jacobiDn(z, m)
        dn*dn - ellipticE(m)/ellipticK(m)

    eJacobiZetaGrad_m(l : List F) : F ==
        z := l(1)
        m := l(2)
        ek := ellipticK(m)
        ee := ellipticE(m)
        er := ee/ek
        dn := jacobiDn(z, m)
        res1 := (dn*dn + m - 1)*jacobiZeta(z, m)
        res2 := res1 + (m - 1)*z*dn*dn
        res3 := res2 - m*jacobiCn(z, m)*jacobiDn(z, m)*jacobiSn(z, m)
        res4 := res3 + z*(1 - m + dn*dn)*er
        ahalf*(res4 - z*er*er)/(m*m - m)

    derivative(opJacobiZeta, [eJacobiZetaGrad_z, eJacobiZetaGrad_m])

    -- handle JacobiTheta

    opJacobiTheta := operator('jacobiTheta)$CommonOperators

    jacobiTheta(q, z) == opJacobiTheta(q, z)

    eJacobiTheta(q : F, z : F) : F ==
        kernel(opJacobiTheta, [q, z])

    elJacobiTheta(l : List F) : F == eJacobiTheta(l(1), l(2))
    evaluate(opJacobiTheta, elJacobiTheta)$BasicOperatorFunctions1(F)

    -- handle LerchPhi

    opLerchPhi := operator('lerchPhi)$CommonOperators

    lerchPhi(z, s, a) == opLerchPhi(z, s, a)

    eLerchPhi(z : F, s : F, a : F) : F ==
        -- z = 0 => 1/a^s
        a = 1 => polylog(s, z)/z
        kernel(opLerchPhi, [z, s, a])

    elLerchPhi(l : List F) : F == eLerchPhi(l(1), l(2), l(3))
    evaluate(opLerchPhi, elLerchPhi)$BasicOperatorFunctions1(F)

    dLerchPhi(l : List F, t : SE) : F ==
        z := l(1)
        s := l(2)
        a := l(3)
        dz := differentiate(z, t)*(lerchPhi(z, s - 1, a) -
                a*lerchPhi(z, s, a))/z
        da := -differentiate(a, t)*s*lerchPhi(z, s + 1, a)
        dm := dummy
        differentiate(s, t)*kernel(opdiff, [opLerchPhi [z, dm, a], dm, s])
           + dz + da

    setProperty(opLerchPhi, SPECIALDIFF, dLerchPhi@((List F, SE)->F)
                                                 pretend None)

    -- handle RiemannZeta

    opRiemannZeta := operator('riemannZeta)$CommonOperators

    riemannZeta(z) == opRiemannZeta(z)

    eRiemannZeta(z : F) : F ==
        kernel(opRiemannZeta, [z])

    elRiemannZeta(l : List F) : F == eRiemannZeta(l(1))
    evaluate(opRiemannZeta, elRiemannZeta)$BasicOperatorFunctions1(F)

    -- orthogonal polynomials

    charlierC(n : F, a : F, z : F) : F == opCharlierC(n, a, z)

    eCharlierC(n : F, a : F, z : F) : F ==
        n = 0 => 1
        n = 1 => (z - a)/a
        kernel(opCharlierC, [n, a, z])

    elCharlierC(l : List F) : F == eCharlierC(l(1), l(2), l(3))

    evaluate(opCharlierC, elCharlierC)$BasicOperatorFunctions1(F)

    hermiteH(n : F, z: F) : F == opHermiteH(n, z)

    eHermiteH(n : F, z: F) : F ==
        n = -1 => 0
        n = 0 => 1
        n = 1 => (2::F)*z
        kernel(opHermiteH, [n, z])

    elHermiteH(l : List F) : F == eHermiteH(l(1), l(2))

    evaluate(opHermiteH, elHermiteH)$BasicOperatorFunctions1(F)

    eHermiteHGrad_z(n : F, z : F) : F == (2::F)*n*hermiteH(n - 1, z)

    dHermiteH(l : List F, t : SE) : F ==
        grad2(l, t, opHermiteH, eHermiteHGrad_z)

    setProperty(opHermiteH, SPECIALDIFF, dHermiteH@((List F, SE)->F)
                                                 pretend None)

    jacobiP(n : F, a : F, b : F, z : F) : F == opJacobiP(n, a, b, z)

    eJacobiP(n : F, a : F, b : F, z : F) : F ==
        n = -1 => 0
        n = 0 => 1
        n = 1 => ahalf*(a - b) + (1 + ahalf*(a + b))*z
        kernel(opJacobiP, [n, a, b, z])

    elJacobiP(l : List F) : F == eJacobiP(l(1), l(2), l(3), l(4))

    evaluate(opJacobiP, elJacobiP)$BasicOperatorFunctions1(F)

    eJacobiPGrad_z(n : F, a : F, b : F, z : F) : F ==
        ahalf*(a + b + n + 1)*jacobiP(n - 1, a + 1, b + 1, z)

    dJacobiP(l : List F, t : SE) : F ==
        grad4(l, t, opJacobiP, eJacobiPGrad_z)

    setProperty(opJacobiP, SPECIALDIFF, dJacobiP@((List F, SE)->F)
                                                 pretend None)

    laguerreL(n : F, a : F, z : F) : F == opLaguerreL(n, a, z)

    eLaguerreL(n : F, a : F, z : F) : F ==
        n = -1 => 0
        n = 0 => 1
        n = 1 => (1 + a - z)
        kernel(opLaguerreL, [n, a, z])

    elLaguerreL(l : List F) : F == eLaguerreL(l(1), l(2), l(3))

    evaluate(opLaguerreL, elLaguerreL)$BasicOperatorFunctions1(F)

    eLaguerreLGrad_z(n : F, a : F, z : F) : F ==
        laguerreL(n - 1, a + 1, z)

    dLaguerreL(l : List F, t : SE) : F ==
        grad3(l, t, opLaguerreL, eLaguerreLGrad_z)

    setProperty(opLaguerreL, SPECIALDIFF, dLaguerreL@((List F, SE)->F)
                                                 pretend None)

    meixnerM(n : F, b : F, c : F, z : F) : F == opMeixnerM(n, b, c, z)

    eMeixnerM(n : F, b : F, c : F, z : F) : F ==
        n = 0 => 1
        n = 1 => (c - 1)*z/(c*b) + 1
        kernel(opMeixnerM, [n, b, c, z])

    elMeixnerM(l : List F) : F == eMeixnerM(l(1), l(2), l(3), l(4))

    evaluate(opMeixnerM, elMeixnerM)$BasicOperatorFunctions1(F)

    --
    belong? op == has?(op, SPECIAL)

    operator op ==
      is?(op, 'abs)      => opabs
      is?(op, 'conjugate)=> opconjugate
      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, 'airyAiPrime) => opAiryAiPrime
      is?(op, 'airyBi)   => opAiryBi
      is?(op, 'airyBiPrime) => opAiryBiPrime
      is?(op, 'lambertW) => opLambertW
      is?(op, 'polylog) => opPolylog
      is?(op, 'weierstrassP)   => opWeierstrassP
      is?(op, 'weierstrassPPrime)   => opWeierstrassPPrime
      is?(op, 'weierstrassSigma)   => opWeierstrassSigma
      is?(op, 'weierstrassZeta)   => opWeierstrassZeta
      is?(op, 'hypergeometricF)  => opHypergeometricF
      is?(op, 'meijerG)  => opMeijerG
      -- is?(op, 'weierstrassPInverse) => opWeierstrassPInverse
      is?(op, 'whittakerM) => opWhittakerM
      is?(op, 'whittakerW) => opWhittakerW
      is?(op, 'angerJ) => opAngerJ
      is?(op, 'weberE) => opWeberE
      is?(op, 'struveH) => opStruveH
      is?(op, 'struveL) => opStruveL
      is?(op, 'hankelH1) => opHankelH1
      is?(op, 'hankelH2) => opHankelH2
      is?(op, 'lommelS1) => opLommelS1
      is?(op, 'lommelS2) => opLommelS2
      is?(op, 'kummerM) => opKummerM
      is?(op, 'kummerU) => opKummerU
      is?(op, 'legendreP) => opLegendreP
      is?(op, 'legendreQ) => opLegendreQ
      is?(op, 'kelvinBei) => opKelvinBei
      is?(op, 'kelvinBer) => opKelvinBer
      is?(op, 'kelvinKei) => opKelvinKei
      is?(op, 'kelvinKer) => opKelvinKer
      is?(op, 'ellipticK) => opEllipticK
      is?(op, 'ellipticE) => opEllipticE
      is?(op, 'ellipticE2) => opEllipticE2
      is?(op, 'ellipticF) => opEllipticF
      is?(op, 'ellipticPi) => opEllipticPi
      is?(op, 'jacobiSn) => opJacobiSn
      is?(op, 'jacobiCn) => opJacobiCn
      is?(op, 'jacobiDn) => opJacobiDn
      is?(op, 'jacobiZeta) => opJacobiZeta
      is?(op, 'jacobiTheta) => opJacobiTheta
      is?(op, 'lerchPhi) => opLerchPhi
      is?(op, 'riemannZeta) => opRiemannZeta
      is?(op, 'charlierC) => opCharlierC
      is?(op, 'hermiteH) => opHermiteH
      is?(op, 'jacobiP) => opJacobiP
      is?(op, 'laguerreL) => opLaguerreL
      is?(op, 'meixnerM) => opMeixnerM
      is?(op, '%logGamma) => op_log_gamma
      is?(op, '%eis) => op_eis
      is?(op, '%erfs) => op_erfs
      is?(op, '%erfis) => op_erfis

      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
      one? x => 1
      is?(x, opabs) => x
      is?(x, opconjugate) => kernel(opabs, argument(retract(x)@K)(1))
      smaller?(x, 0) => kernel(opabs, -x)
      kernel(opabs, x)

    if R has abs : R -> R then
      import from Polynomial R
      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
      if F has RadicalCategory and R has conjugate : R -> R then
        iiabs x ==
          (r := retractIfCan(x)@Union(R, "failed"))
            case "failed" => iabs x
          sqrt( (r::R*conjugate(r::R))::F)
      else iiabs x == iabs x

    iconjugate(k:K):F ==
      --output("in: ",k::OutputForm)$OutputPackage
      if symbolIfCan(k) case Symbol then
        x:=kernel(opconjugate, k::F)
      else
        op:=operator(operator(name(k))$CommonOperators)$Expression(R)
        --output("op: ",properties(op)::OutputForm)$OutputPackage
        --output("conj op: ",properties(conjugate op)::OutputForm)$OutputPackage
        x := conjugate(op) argument(k)
      --output("out: ",x::OutputForm)$OutputPackage
      return x

    iiconjugate(x:F):F ==
      --output("iin: ",x::OutputForm)$OutputPackage
      if (s:=isPlus(x)) case List F then
        --output("isPlus: ",s::OutputForm)$OutputPackage
        x := reduce(_+$F,map(iiconjugate,s))
      else if (s:=isTimes(x)) case List F then
        --output("isTimes: ",s::OutputForm)$OutputPackage
        x:= reduce(_*$F,map(iiconjugate,s))
      else if #(ks:List K:=kernels(x))>0 then
        --output("kernels: ",ks::OutputForm)$OutputPackage
        x := eval(x,ks, _
          map((k:Kernel F):F +-> (height(k)=0 =>k::F;iconjugate k), _
            ks)$ListFunctions2(Kernel F,F))
      else if R has conjugate:R->R and F has RetractableTo(R) then
        if (r:=retractIfCan(x)@Union(R,"failed")) case R then
          x:=conjugate(r::R)::F
      --output("iout: ",x::OutputForm)$OutputPackage
      return x

    iBeta(x, y) == kernel(opBeta, [x, y])
    idigamma x == kernel(opdigamma, x)
    iiipolygamma(n, x) == kernel(oppolygamma, [n, x])
    iiiBesselJ(x, y) == kernel(opBesselJ, [x, y])
    iiiBesselY(x, y) == kernel(opBesselY, [x, y])
    iiiBesselI(x, y) == kernel(opBesselI, [x, y])
    iiiBesselK(x, y) == kernel(opBesselK, [x, y])

    import from Fraction(Integer)

    if F has ElementaryFunctionCategory then
        iAiryAi x ==
            zero?(x) => 1::F/((3::F)^twothirds*Gamma(twothirds))
            kernel(opAiryAi, x)
        iAiryAiPrime x ==
            zero?(x) => -1::F/((3::F)^athird*Gamma(athird))
            kernel(opAiryAiPrime, x)
        iAiryBi x ==
            zero?(x) => 1::F/((3::F)^asixth*Gamma(twothirds))
            kernel(opAiryBi, x)
        iAiryBiPrime x ==
            zero?(x) => (3::F)^asixth/Gamma(athird)
            kernel(opAiryBiPrime, x)
    else
        iAiryAi x == kernel(opAiryAi, x)
        iAiryAiPrime x == kernel(opAiryAiPrime, x)
        iAiryBi x == kernel(opAiryBi, x)
        iAiryBiPrime x == kernel(opAiryBiPrime, x)

    if F has ElementaryFunctionCategory then
        iLambertW(x) ==
            zero?(x) => 0
            x = exp(1$F) => 1$F
            x = -exp(-1$F) => -1$F
            kernel(opLambertW, x)
    else
        iLambertW(x) ==
            zero?(x) => 0
            kernel(opLambertW, x)

    if F has ElementaryFunctionCategory then
        if F has LiouvillianFunctionCategory then
            iiPolylog(s, x) ==
                s = 1 => -log(1 - x)
                s = 2::F => dilog(1 - x)
                kernel(opPolylog, [s, x])
        else
            iiPolylog(s, x) ==
                s = 1 => -log(1 - x)
                kernel(opPolylog, [s, x])
    else
        iiPolylog(s, x) == kernel(opPolylog, [s, x])

    iPolylog(l) == iiPolylog(first l, second l)

    iWeierstrassP(g2, g3, x) == kernel(opWeierstrassP, [g2, g3, x])
    iWeierstrassPPrime(g2, g3, x) == kernel(opWeierstrassPPrime, [g2, g3, x])

    iWeierstrassSigma(g2, g3, x) ==
        x = 0 => 0
        kernel(opWeierstrassSigma, [g2, g3, x])

    iWeierstrassZeta(g2, g3, x) == kernel(opWeierstrassZeta, [g2, g3, x])

    -- Could put more conditional special rules for other functions here

    if R has SpecialFunctionCategory then
      iiGamma x ==
        (r := retractIfCan(x)@Union(R,"failed")) case "failed" => iGamma x
        Gamma(r::R)::F

      iiBeta l ==
        (r := retractIfCan(first l)@Union(R,"failed")) case "failed" or _
        (s := retractIfCan(second l)@Union(R,"failed")) case "failed" _
            => iBeta(first l, second l)
        Beta(r::R, s::R)::F

      iidigamma x ==
        (r := retractIfCan(x)@Union(R,"failed")) case "failed" => idigamma x
        digamma(r::R)::F

      iipolygamma l ==
        (s := retractIfCan(first l)@Union(R,"failed")) case "failed" or _
        (r := retractIfCan(second l)@Union(R,"failed")) case "failed" _
            => iiipolygamma(first l, second l)
        polygamma(s::R, r::R)::F

      iiBesselJ l ==
        (r := retractIfCan(first l)@Union(R,"failed")) case "failed" or _
        (s := retractIfCan(second l)@Union(R,"failed")) case "failed" _
            => iiiBesselJ(first l, second l)
        besselJ(r::R, s::R)::F

      iiBesselY l ==
        (r := retractIfCan(first l)@Union(R,"failed")) case "failed" or _
        (s := retractIfCan(second l)@Union(R,"failed")) case "failed" _
            => iiiBesselY(first l, second l)
        besselY(r::R, s::R)::F

      iiBesselI l ==
        (r := retractIfCan(first l)@Union(R,"failed")) case "failed" or _
        (s := retractIfCan(second l)@Union(R,"failed")) case "failed" _
            => iiiBesselI(first l, second l)
        besselI(r::R, s::R)::F

      iiBesselK l ==
        (r := retractIfCan(first l)@Union(R,"failed")) case "failed" or _
        (s := retractIfCan(second l)@Union(R,"failed")) case "failed" _
            => iiiBesselK(first l, second l)
        besselK(r::R, s::R)::F

      iiAiryAi x ==
        (r := retractIfCan(x)@Union(R,"failed")) case "failed" => iAiryAi x
        airyAi(r::R)::F

      iiAiryAiPrime x ==
        (r := retractIfCan(x)@Union(R,"failed")) case "failed" => iAiryAiPrime x
        airyAiPrime(r::R)::F

      iiAiryBi x ==
        (r := retractIfCan(x)@Union(R,"failed")) case "failed" => iAiryBi x
        airyBi(r::R)::F

      iiAiryBi x ==
        (r := retractIfCan(x)@Union(R,"failed")) case "failed" => iAiryBiPrime x
        airyBiPrime(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

      iiBeta l == iBeta(first l, second l)
      iidigamma x == idigamma x
      iipolygamma l == iiipolygamma(first l, second l)
      iiBesselJ l == iiiBesselJ(first l, second l)
      iiBesselY l == iiiBesselY(first l, second l)
      iiBesselI l == iiiBesselI(first l, second l)
      iiBesselK l == iiiBesselK(first l, second l)
      iiAiryAi x == iAiryAi x
      iiAiryAiPrime x == iAiryAiPrime x
      iiAiryBi x == iAiryBi x
      iiAiryBiPrime x == iAiryBiPrime x

    iiWeierstrassP l == iWeierstrassP(first l, second l, third l)
    iiWeierstrassPPrime l == iWeierstrassPPrime(first l, second l, third l)
    iiWeierstrassSigma l == iWeierstrassSigma(first l, second l, third l)
    iiWeierstrassZeta l == iWeierstrassZeta(first l, second l, third l)

    -- Default behaviour is to build a kernel
    evaluate(opGamma, iiGamma)$BasicOperatorFunctions1(F)
    evaluate(opabs, iiabs)$BasicOperatorFunctions1(F)
    evaluate(conjugate(opabs), (x:F):F +-> iiabs(x))$BasicOperatorFunctions1(F)
    evaluate(opconjugate, iiconjugate)$BasicOperatorFunctions1(F)
    evaluate(conjugate(opconjugate), (x:F):F +-> x)$BasicOperatorFunctions1(F)
--    evaluate(opGamma2    , iiGamma2   )$BasicOperatorFunctions1(F)
    evaluate(opBeta      , iiBeta     )$BasicOperatorFunctions1(F)
    evaluate(opdigamma   , iidigamma  )$BasicOperatorFunctions1(F)
    evaluate(oppolygamma , iipolygamma)$BasicOperatorFunctions1(F)
    evaluate(opBesselJ   , iiBesselJ  )$BasicOperatorFunctions1(F)
    evaluate(opBesselY   , iiBesselY  )$BasicOperatorFunctions1(F)
    evaluate(opBesselI   , iiBesselI  )$BasicOperatorFunctions1(F)
    evaluate(opBesselK   , iiBesselK  )$BasicOperatorFunctions1(F)
    evaluate(opAiryAi    , iiAiryAi   )$BasicOperatorFunctions1(F)
    evaluate(opAiryAiPrime, iiAiryAiPrime)$BasicOperatorFunctions1(F)
    evaluate(opAiryBi    , iiAiryBi   )$BasicOperatorFunctions1(F)
    evaluate(opAiryBiPrime, iiAiryBiPrime)$BasicOperatorFunctions1(F)
    evaluate(opLambertW, iLambertW)$BasicOperatorFunctions1(F)
    evaluate(opPolylog, iPolylog)$BasicOperatorFunctions1(F)
    evaluate(opWeierstrassP, iiWeierstrassP)$BasicOperatorFunctions1(F)
    evaluate(opWeierstrassPPrime,
             iiWeierstrassPPrime)$BasicOperatorFunctions1(F)
    evaluate(opWeierstrassSigma, iiWeierstrassSigma)$BasicOperatorFunctions1(F)
    evaluate(opWeierstrassZeta, iiWeierstrassZeta)$BasicOperatorFunctions1(F)
    evaluate(opHypergeometricF, iiHypergeometricF)$BasicOperatorFunctions1(F)
    evaluate(opMeijerG, iiMeijerG)$BasicOperatorFunctions1(F)

    diff1(op : OP, n : F, x : F) : F ==
        dm := dummy
        kernel(opdiff, [op [dm, x], dm, n])

    iBesselJ(l : List F, t : SE) : F ==
        n := first l; x := second l
        differentiate(n, t)*diff1(opBesselJ, n, x)
          + differentiate(x, t) * ahalf * (besselJ (n-1, x) - besselJ (n+1, x))

    iBesselY(l : List F, t : SE) : F ==
        n := first l; x := second l
        differentiate(n, t)*diff1(opBesselY, n, x)
          + differentiate(x, t) * ahalf * (besselY (n-1, x) - besselY (n+1, x))

    iBesselI(l : List F, t : SE) : F ==
        n := first l; x := second l
        differentiate(n, t)*diff1(opBesselI, n, x)
          + differentiate(x, t)* ahalf * (besselI (n-1, x) + besselI (n+1, x))

    iBesselK(l : List F, t : SE) : F ==
        n := first l; x := second l
        differentiate(n, t)*diff1(opBesselK, n, x)
          - differentiate(x, t)* ahalf * (besselK (n-1, x) + besselK (n+1, x))

    dPolylog(l : List F, t : SE) : F ==
        s := first l; x := second l
        differentiate(s, t)*diff1(opPolylog, s, x)
           + differentiate(x, t)*polylog(s-1, x)/x

    ipolygamma(l : List F, x : SE) : F ==
        import from List(Symbol)
        member?(x, variables first l) =>
            error "cannot differentiate polygamma with respect to the first argument"
        n := first l; y := second l
        differentiate(y, x)*polygamma(n+1, y)
    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
      iGamma2(l : List F, t : SE) : F ==
        a := first l; x := second l
        differentiate(a, t)*diff1(opGamma2, a, x)
          - differentiate(x, t)* x ^ (a - 1) * exp(-x)
      setProperty(opGamma2, SPECIALDIFF, iGamma2@((List F, SE)->F)
                                                 pretend None)

    inGamma2(li : List INP) : INP ==
        convert cons(convert('Gamma), li)

    input(opGamma2, inGamma2@((List INP) -> INP))

    dLambertW(x : F) : F ==
        lw := lambertW(x)
        lw/(x*(1+lw))

    iWeierstrassPGrad1(l : List F) : F ==
        g2 := first l
        g3 := second l
        x := third l
        delta := g2^3 - 27*g3^2
        wp := weierstrassP(g2, g3, x)
        (weierstrassPPrime(g2, g3, x)*(-9*ahalf*g3
          *weierstrassZeta(g2, g3, x) + afourth*g2^2*x)
            - 9*g3*wp^2 + ahalf*g2^2*wp + 3*ahalf*g2*g3)/delta

    iWeierstrassPGrad2(l : List F) : F ==
        g2 := first l
        g3 := second l
        x := third l
        delta := g2^3 - 27*g3^2
        wp := weierstrassP(g2, g3, x)
        (weierstrassPPrime(g2, g3, x)*(3*g2*weierstrassZeta(g2, g3, x)
          - 9*ahalf*g3*x) + 6*g2*wp^2 - 9*g3*wp-g2^2)/delta

    iWeierstrassPGrad3(l : List F) : F ==
        weierstrassPPrime(first l, second l, third l)

    iWeierstrassPPrimeGrad1(l : List F) : F ==
        g2 := first l
        g3 := second l
        x := third l
        delta := g2^3 - 27*g3^2
        wp := weierstrassP(g2, g3, x)
        wpp := weierstrassPPrime(g2, g3, x)
        wpp2 := 6*wp^2 - ahalf*g2
        (wpp2*(-9*ahalf*g3*weierstrassZeta(g2, g3, x) + afourth*g2^2*x)
           + wpp*(9*ahalf*g3*wp + afourth*g2^2) - 18*g3*wp*wpp
             + ahalf*g2^2*wpp)/delta

    iWeierstrassPPrimeGrad2(l : List F) : F ==
        g2 := first l
        g3 := second l
        x := third l
        delta := g2^3 - 27*g3^2
        wp := weierstrassP(g2, g3, x)
        wpp := weierstrassPPrime(g2, g3, x)
        wpp2 := 6*wp^2 - ahalf*g2
        (wpp2*(3*g2*weierstrassZeta(g2, g3, x) - 9*ahalf*g3*x)
          + wpp*(-3*g2*wp - 9*ahalf*g3) + 12*g2*wp*wpp - 9*g3*wpp)/delta

    iWeierstrassPPrimeGrad3(l : List F) : F ==
        g2 := first l
        6*weierstrassP(g2, second l, third l)^2 - ahalf*g2

    iWeierstrassSigmaGrad1(l : List F) : F ==
        g2 := first l
        g3 := second l
        x := third l
        delta := g2^3 - 27*g3^2
        ws := weierstrassSigma(g2, g3, x)
        wz := weierstrassZeta(g2, g3, x)
        wsp := wz*ws
        wsp2 := - weierstrassP(g2, g3, x)*ws + wz^2*ws
        afourth*(-9*g3*wsp2 - g2^2*ws
                  - 3*afourth*g2*g3*x^2*ws + g2^2*x*wsp)/delta

    iWeierstrassSigmaGrad2(l : List F) : F ==
        g2 := first l
        g3 := second l
        x := third l
        delta := g2^3 - 27*g3^2
        ws := weierstrassSigma(g2, g3, x)
        wz := weierstrassZeta(g2, g3, x)
        wsp := wz*ws
        wsp2 := - weierstrassP(g2, g3, x)*ws + wz^2*ws
        ahalf*(3*g2*wsp2 + 9*g3*ws
                 + afourth*g2^2*x^2*ws - 9*g3*x*wsp)/delta

    iWeierstrassSigmaGrad3(l : List F) : F ==
        g2 := first l
        g3 := second l
        x := third l
        weierstrassZeta(g2, g3, x)*weierstrassSigma(g2, g3, x)

    iWeierstrassZetaGrad1(l : List F) : F ==
        g2 := first l
        g3 := second l
        x := third l
        delta := g2^3 - 27*g3^2
        wp := weierstrassP(g2, g3, x)
        (ahalf*weierstrassZeta(g2, g3, x)*(9*g3*wp + ahalf*g2^2)
             - ahalf*g2*x*(ahalf*g2*wp+3*afourth*g3)
               + 9*afourth*g3*weierstrassPPrime(g2, g3, x))/delta

    iWeierstrassZetaGrad2(l : List F) : F ==
        g2 := first l
        g3 := second l
        x := third l
        delta := g2^3 - 27*g3^2
        wp := weierstrassP(g2, g3, x)
        (-3*weierstrassZeta(g2, g3, x)*(g2*wp + 3*ahalf*g3) +
                     ahalf*x*(9*g3*wp + ahalf*g2^2)
                       - 3*ahalf*g2*weierstrassPPrime(g2, g3, x))/delta

    iWeierstrassZetaGrad3(l : List F) : F ==
        -weierstrassP(first l, second l, third l)

    OF ==> OutputForm
    SEX ==> SExpression
    NNI ==> NonNegativeInteger

    dconjugate(lo : List OF) : OF == overbar lo.1
    display(opconjugate,dconjugate)

    if F has RetractableTo(Integer) then
        get_int_listf : List F -> List Integer
        get_int_listo : (Integer, List OF) ->  List Integer
        get_int_listi : (Integer, List INP) ->  List Integer

        get_int_listf(lf : List F) : List Integer ==
            map(z +-> retract(z)@Integer, lf)$ListFunctions2(F, Integer)

        replace_i(lp : List F, v : F, i : NNI) : List F ==
            concat(first(lp, (i - 1)::NNI), cons(v, rest(lp, i)))

        iiHypergeometricF(l) ==
            n := #l
            z := l(n-2)
            if z = 0 then
                nn := (n - 2)::NNI
                pq := rest(l, nn)
                pqi := get_int_listf(pq)
                p := first(pqi)
                q := first(rest(pqi))
                p <= q + 1 => return 1
            kernel(opHypergeometricF, l)

        idvsum(op : BasicOperator, n : Integer, l : List F, x : Symbol) : F  ==
            res : F := 0
            for i in 1..n for a in l repeat
                dm := dummy
                nl := replace_i(l, dm, i)
                res := res + differentiate(a, x)*kernel(opdiff, [op nl, dm, a])
            res

        dvhypergeom(l : List F, x : Symbol) : F ==
            n := #l
            nn := (n - 2)::NNI
            pq := rest(l, nn)
            pqi := get_int_listf(pq)
            ol := l
            l := first(l, nn)
            l1 := reverse(l)
            z := first(l1)
            p := first(pqi)
            q := first(rest(pqi))
            aprod := 1@F
            nl := []@(List F)
            for i in 1..p repeat
                a := first(l)
                nl := cons(a + 1, nl)
                aprod := aprod * a
                l := rest(l)
            bprod := 1@F
            for i in 1..q repeat
                b := first(l)
                nl := cons(b + 1, nl)
                bprod := bprod * b
                l := rest(l)
            nl0 := reverse!(nl)
            nl1 := cons(z, pq)
            nl := concat(nl0, nl1)
            aprod := aprod/bprod
            idvsum(opHypergeometricF, nn - 1, ol, x) +
                 differentiate(z, x)*aprod*opHypergeometricF(nl)

        add_pairs_to_list(lp : List List F, l : List F) : List F ==
            for p in lp repeat
                #p ~= 2 => error "not a list of pairs"
                l := cons(p(2), cons(p(1), l))
            l

        dvmeijer(l : List F, x : Symbol) : F ==
            n := #l
            nn := (n - 4)::NNI
            l0 := l
            nl := rest(l, nn)
            nli := get_int_listf(nl)
            l := first(l, nn)
            l1 := reverse(l)
            z := first(l1)
            n1 := first(nli)
            n2 := nli(2)
            a := first l
            sign : F := 1
            if n1 > 0 or n2 > 0 then
                 na := a - 1
                 if n1 = 0 then sign := -1
                 l2 := cons(na, rest l)
            else
                 na := a
                 if nli(3) > 0 then sign := -1
                 l2 := cons(a + 1, rest l)
            nm : F := opMeijerG(concat(l2, nl))
            om : F := opMeijerG(l0)
            idvsum(opMeijerG, nn - 1, l0, x) +
                differentiate(z, x)*(sign*nm + na*om)/z

        get_if_list(n : Integer, lf : List INP) : List List INP ==
            a := []@(List INP)
            for i in 1..n repeat
                a := cons(first(lf), a)
                lf := rest(lf)
            a := cons(convert('construct), reverse!(a))
            [a, lf]

        get_if_lists(ln : List Integer, lf : List INP) : List List INP ==
            rl := []@(List List INP)
            for n in ln repeat
                al := get_if_list(n, lf)
                rl := cons(first(al), rl)
                lf := first(rest(al))
            rl := reverse!(rl)
            cons(lf, rl)

        get_int_listi(n : Integer, lo : List INP) : List Integer ==
            n0 := (#lo - n)::NNI
            lo := rest(lo, n0)
            rl := []@(List Integer)
            for i in 1..n repeat
                p := integer(first(lo) pretend SEX)$SEX
                rl := cons(p, rl)
                lo := rest(lo)
            rl := reverse!(rl)
            rl

        get_of_list(n : Integer, lo : List OF) : List List OF ==
            a := []@(List OF)
            for i in 1..n repeat
                a := cons(first(lo), a)
                lo := rest(lo)
            a := reverse!(a)
            [a, lo]

        get_of_lists(ln : List Integer, lo : List OF) : List List OF ==
            rl := []@(List List OF)
            for n in ln repeat
                al := get_of_list(n, lo)
                rl := cons(first(al), rl)
                lo := first(rest(al))
            rl := reverse!(rl)
            cons(lo, rl)

        get_int_listo(n : Integer, lo : List OF) : List Integer ==
            n0 := (#lo - n)::NNI
            lo := rest(lo, n0)
            rl := []@(List Integer)
            for i in 1..n repeat
                p := integer(first(lo) pretend SEX)$SEX
                rl := cons(p, rl)
                lo := rest(lo)
            rl := reverse!(rl)
            rl

        dhyper0(op : OF, lo : List OF) : OF ==
            n0 := (#lo - 2)::NNI
            pql := get_int_listo(2, lo)
            lo := first(lo, n0)
            al := get_of_lists(pql, lo)
            lo := first(al)
            al := rest(al)
            a := first al
            b := first(rest(al))
            z := first(lo)
            prefix(op, [bracket a, bracket b, z])

        dhyper(lo : List OF) : OF ==
            dhyper0("hypergeometricF"::Symbol::OF, lo)

        ddhyper(lo : List OF) : OF ==
            dhyper0(first lo, rest lo)

        dmeijer0(op : OF, lo : List OF) : OF ==
            n0 := (#lo - 4)::NNI
            nl := get_int_listo(4, lo)
            lo := first(lo, n0)
            al := get_of_lists(nl, lo)
            lo := first(al)
            al := rest(al)
            z := first(lo)
            prefix(op, concat(
                         map(bracket, al)$ListFunctions2(List OF, OF), [z]))

        dmeijer(lo : List OF) : OF ==
            dmeijer0('meijerG::OF, lo)

        ddmeijer(lo : List OF) : OF ==
            dmeijer0(first lo, rest lo)

        setProperty(opHypergeometricF, '%diffDisp,
                     ddhyper@(List OF -> OF) pretend None)
        setProperty(opMeijerG, '%diffDisp,
                     ddmeijer@(List OF -> OF) pretend None)
        display(opHypergeometricF, dhyper)
        display(opMeijerG, dmeijer)
        setProperty(opHypergeometricF, SPECIALDIFF,
                 dvhypergeom@((List F, Symbol)->F) pretend None)
        setProperty(opMeijerG, SPECIALDIFF, dvmeijer@((List F, Symbol)->F)
                                            pretend None)

        inhyper(lf : List INP) : INP ==
            pqi := get_int_listi(2, lf)
            al := get_if_lists(pqi, lf)
            lf := first(al)
            al := rest(al)
            a := first al
            ai : INP := convert(a)
            b := first(rest(al))
            bi : INP := convert(b)
            zi := first(lf)
            li : List INP := [convert('hypergeometricF), ai, bi, zi]
            convert(li)

        input(opHypergeometricF, inhyper@((List INP) -> INP))

        inmeijer(lf : List INP) : INP ==
            pqi := get_int_listi(4, lf)
            al := get_if_lists(pqi, lf)
            lf := first(al)
            al := rest(al)
            a := first al
            ai : INP := convert(a)
            al := rest(al)
            b := first(al)
            bi : INP := convert(b)
            al := rest(al)
            c := first(al)
            ci : INP := convert(c)
            al := rest(al)
            d := first(al)
            di : INP := convert(d)
            zi := first(lf)
            li : List INP := [convert('meijerG), ai, bi, ci, di, zi]
            convert(li)

        input(opMeijerG, inmeijer@((List INP) -> INP))

    else
        iiHypergeometricF(l) == kernel(opHypergeometricF, l)

    iiMeijerG(l) == kernel(opMeijerG, l)

    d_eis(x : F) : F == -kernel(op_eis, x) + 1/x

    if F has TranscendentalFunctionCategory
       and F has RadicalCategory then

        d_erfs(x : F) : F == 2*x*kernel(op_erfs, x) - 2::F/sqrt(pi())

        d_erfis(x : F) : F == -2*x*kernel(op_erfis, x) + 2::F/sqrt(pi())

        derivative(op_erfs,     d_erfs)
        derivative(op_erfis,    d_erfis)

    -- differentiate(abs(g(z)),z)
    --dvabs     : (List F, SE) -> F
    --dvabs(arg, z) ==
    --  g := first arg
    --  conjugate(g)*inv(2*abs(g))*differentiate(g,z)+g*inv(2*abs(g))*differentiate(conjugate(g),z)
    --setProperty(opabs,  SPECIALDIFF, dvabs@((List F, SE)->F) pretend None)
    derivative(opabs,       (x : F) : F +-> conjugate(x)*inv(2*abs(x)))

    -- Wirtinger derivate of non-holomorphic functions
    derivative(opconjugate, (x : F) : F +-> 0)
    derivative(conjugate opconjugate, (x : F) : F +-> 1)
    derivative(opGamma,     (x : F) : F +-> digamma(x)*Gamma(x))
    derivative(op_log_gamma, (x : F) : F +-> digamma(x))
    derivative(opBeta,      [iBetaGrad1, iBetaGrad2])
    derivative(opdigamma,   (x : F) : F +-> polygamma(1, x))
    derivative(op_eis,      d_eis)
    derivative(opAiryAi,    (x : F) : F +-> airyAiPrime(x))
    derivative(opAiryAiPrime,    (x : F) : F +-> x*airyAi(x))
    derivative(opAiryBi,    (x : F) : F +-> airyBiPrime(x))
    derivative(opAiryBiPrime,    (x : F) : F +-> x*airyBi(x))
    derivative(opLambertW,  dLambertW)
    derivative(opWeierstrassP, [iWeierstrassPGrad1, iWeierstrassPGrad2,
                                iWeierstrassPGrad3])
    derivative(opWeierstrassPPrime, [iWeierstrassPPrimeGrad1,
                   iWeierstrassPPrimeGrad2, iWeierstrassPPrimeGrad3])
    derivative(opWeierstrassSigma, [iWeierstrassSigmaGrad1,
                   iWeierstrassSigmaGrad2, iWeierstrassSigmaGrad3])
    derivative(opWeierstrassZeta, [iWeierstrassZetaGrad1,
                   iWeierstrassZetaGrad2, iWeierstrassZetaGrad3])

    setProperty(oppolygamma, SPECIALDIFF, ipolygamma@((List F, SE)->F)
                                                     pretend None)
    setProperty(opBesselJ, SPECIALDIFF, iBesselJ@((List F, SE)->F)
                                                 pretend None)
    setProperty(opBesselY, SPECIALDIFF, iBesselY@((List F, SE)->F)
                                                 pretend None)
    setProperty(opBesselI, SPECIALDIFF, iBesselI@((List F, SE)->F)
                                                 pretend None)
    setProperty(opBesselK, SPECIALDIFF, iBesselK@((List F, SE)->F)
                                                 pretend None)
    setProperty(opPolylog, SPECIALDIFF, dPolylog@((List F, SE)->F)
                                                 pretend None)

)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, Comparable,
          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 from ElementaryFunctionStructurePackage(R, F)
    import from GosperSummationMethod(IndexedExponents K, K, R,
                                 SparseMultivariatePolynomial(R, K), F)
    import from Segment(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)

--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.

-- SPAD files for the functional world should be compiled in the
-- following order:
--
--   op  kl  function  funcpkgs  manip  algfunc
--   elemntry  constant  funceval  COMBFUNC  fe

   Compiling FriCAS source code from file 
      /var/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/554697800083835704-25px002.spad
      using old system compiler.
   COMBOPC abbreviates category CombinatorialOpsCategory 
------------------------------------------------------------------------
   initializing NRLIB COMBOPC for CombinatorialOpsCategory 
   compiling into NRLIB COMBOPC 

;;;     ***       |CombinatorialOpsCategory| REDEFINED
Time: 0 SEC.

   finalizing NRLIB COMBOPC 
   Processing CombinatorialOpsCategory for Browser database:
--------constructor---------
--------(factorials (% %))---------
--------(factorials (% % (Symbol)))---------
--------(summation (% % (Symbol)))---------
--------(summation (% % (SegmentBinding %)))---------
--------(product (% % (Symbol)))---------
--------(product (% % (SegmentBinding %)))---------
; compiling file "/var/aw/var/LatexWiki/COMBOPC.NRLIB/COMBOPC.lsp" (written 07 OCT 2024 08:03:59 PM):

; wrote /var/aw/var/LatexWiki/COMBOPC.NRLIB/COMBOPC.fasl
; compilation finished in 0:00:00.004
------------------------------------------------------------------------
   CombinatorialOpsCategory is now explicitly exposed in frame initial 
   CombinatorialOpsCategory will be automatically loaded when needed 
      from /var/aw/var/LatexWiki/COMBOPC.NRLIB/COMBOPC

   COMBF abbreviates package CombinatorialFunction 
------------------------------------------------------------------------
   initializing NRLIB COMBF for CombinatorialFunction 
   compiling into NRLIB COMBF 
****** Domain: R already in scope
   processing macro definition dummy ==> ::((Sel (Symbol) new),F) 
   compiling exported factorial : F -> F
Time: 0.02 SEC.

   compiling exported binomial : (F,F) -> F
Time: 0 SEC.

   compiling exported permutation : (F,F) -> F
Time: 0 SEC.

   importing F
   importing Kernel F
   compiling local number? : F -> Boolean
****** Domain: R already in scope
augmenting R: (RetractableTo (Integer))
Time: 0 SEC.

   compiling exported ^ : (F,F) -> F
Time: 0.02 SEC.

   compiling exported belong? : BasicOperator -> Boolean
****** comp fails at level 2 with expression: ******
error in function belong? 

(|has?| | << op >> | COMB)
****** level 2  ******
$x:= op
$m:= $
$f:=
((((|op| # #) (|number?| #)) ((|number?| #) (|/\\| #) (< #) (<= #) ...)))

   >> Apparent user error:
   Cannot coerce op 
      of mode (BasicOperator) 
      to mode $

SandBoxCombfuncDiscussion