axiom-source 20170501-3 / usr / share / axiom-20170501 / src / algebra / COMBF.spad

)abbrev package COMBF CombinatorialFunction
++ Author: Manuel Bronstein, Martin Rubey
++ Date Created: 2 Aug 1988
++ Date Last Updated: 30 October 2005
++ Description:
++ Provides combinatorial functions over an integral domain.

CombinatorialFunction(R, F) : SIG == CODE where
  R : Join(OrderedSet, IntegralDomain)
  F : FunctionSpace R

  OP  ==> BasicOperator
  K   ==> Kernel F
  SE  ==> Symbol
  O   ==> OutputForm
  SMP ==> SparseMultivariatePolynomial(R, K)
  Z   ==> Integer

  POWER        ==> "%power"::Symbol
  OPEXP        ==> "exp"::Symbol
  SPECIALDIFF  ==> "%specialDiff"
  SPECIALDISP  ==> "%specialDisp"
  SPECIALEQUAL ==> "%specialEqual"

  SIG ==> 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, n!/(r! * (n-r)!);
      ++
      ++X [binomial(5,i) for i in 0..5]

    permutation : (F, F) -> F
      ++ permutation(n, r) returns the number of permutations of
      ++ n objects taken r at a time, n!/(n-r)!;

    factorial : F -> F
      ++ factorial(n) returns the factorial of n, 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;

  CODE ==> add

    ifact     : F -> F
    iiipow    : List F -> F
    iperm     : List F -> F
    ibinom    : List F -> F
    isum      : List F -> F
    idsum     : List F -> F
    iprod     : List F -> F
    idprod    : List F -> F
    dsum      : List F -> O
    ddsum     : List F -> O
    dprod     : List F -> O
    ddprod    : List F -> O
    equalsumprod  : (K, K) -> Boolean 
    equaldsumprod : (K, K) -> Boolean 
    fourth    : List F -> F
    dvpow1    : List F -> F
    dvpow2    : List F -> F
    summand   : List F -> F
    dvsum     : (List F, SE) -> F
    dvdsum    : (List F, SE) -> F
    dvprod    : (List F, SE) -> F
    dvdprod   : (List F, SE) -> F
    facts     : (F, List SE) -> F
    K2fact    : (K, List SE) -> F
    smpfact   : (SMP, List SE) -> F

-- This macro will be used in product and summation, both the 5 and 3
-- argument forms. It is used to introduce a dummy variable in place of the
-- summation index within the summands. This in turn is necessary to keep the
-- indexing variable local, circumventing problems, for example, with
-- differentiation.

    dummy == new()$SE :: F

    opfact  := operator("factorial"::Symbol)$CommonOperators

    opperm  := operator("permutation"::Symbol)$CommonOperators

    opbinom := operator("binomial"::Symbol)$CommonOperators

    opsum   := operator("summation"::Symbol)$CommonOperators

    opdsum  := operator("%defsum"::Symbol)$CommonOperators

    opprod  := operator("product"::Symbol)$CommonOperators

    opdprod := operator("%defprod"::Symbol)$CommonOperators

    oppow   := operator(POWER::Symbol)$CommonOperators

    factorial x          == opfact x

    binomial(x, y)       == opbinom [x, y]

    permutation(x, y)    == opperm [x, y]

    import F
    import Kernel F

    number?(x:F):Boolean ==
      if R has RetractableTo(Z) then
        ground?(x) or
         ((retractIfCan(x)@Union(Fraction(Z),"failed")) case Fraction(Z))
      else
        ground?(x)

    x ** y               == 
      -- Do some basic simplifications
      is?(x,POWER) =>
        args : List F := argument first kernels x
        not(#args = 2) => error "Too many arguments to **"
        number?(first args) and number?(y) =>
          oppow [first(args)**y, second args]
        oppow [first args, (second args)* y]
      -- Generic case
      exp : Union(Record(val:F,exponent:Z),"failed") := isPower x
      exp case Record(val:F,exponent:Z) =>
        expr := exp::Record(val:F,exponent:Z)
        oppow [expr.val, (expr.exponent)*y]
      oppow [x, y]

    belong? op           == has?(op, "comb")

    fourth l             == third rest l

    dvpow1 l             == second(l) * first(l) ** (second l - 1)

    factorials x         == facts(x, variables x)

    factorials(x, v)     == facts(x, [v])

    facts(x, l)          == smpfact(numer x, l) / smpfact(denom x, l)

    summand l            == eval(first l, retract(second l)@K, third l)

    product(x:F, i:SE) ==
      dm := dummy
      opprod [eval(x, k := kernel(i)$K, dm), dm, k::F]

    summation(x:F, i:SE) ==
      dm := dummy
      opsum [eval(x, k := kernel(i)$K, dm), dm, k::F]

-- These two operations return the product or the sum as unevaluated operators
-- A dummy variable is introduced to make the indexing variable local.

    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 

-- These four operations handle the conversion of sums and products to
-- OutputForm

    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)

-- The two operations handle the testing for equality of sums and products.
-- The corresponding property \verb|%specialEqual| set below is checked in
-- Kernel. Note that we can assume that the operators are equal, since this is
-- checked in Kernel itself.

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

-- These two operations return the product or the sum as unevaluated operators
-- A dummy variable is introduced to make the indexing variable local.

    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), y+->y::F, p)_
        $PolynomialCategoryLifting(IndexedExponents K, K, R, SMP, F)

    K2fact(k, l) ==
      empty? [v for v in variables(kf := k::F) | member?(v, l)] => kf
      empty?(args:List F := [facts(a, l) for a in argument k]) => kf
      is?(k, opperm) =>
        factorial(n := first args) / factorial(n - second args)
      is?(k, opbinom) =>
        n := first args
        p := second args
        factorial(n) / (factorial(p) * factorial(n-p))
      (operator k) args

    operator op ==
      is?(op, "factorial"::Symbol)   => opfact
      is?(op, "permutation"::Symbol) => opperm
      is?(op, "binomial"::Symbol)    => opbinom
      is?(op, "summation"::Symbol)   => opsum
      is?(op, "%defsum"::Symbol)     => opdsum
      is?(op, "product"::Symbol)     => opprod
      is?(op, "%defprod"::Symbol)    => opdprod
      is?(op, POWER)                 => oppow
      error "Not a combinatorial operator"

    iprod l ==
      zero? first l => 0
      (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 (x = 1) => 1
      kernel(opfact, x)

    ibinom l ==
      n := first l
      ((p := second l) = 0) or (p = n) => 1
      (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
        (x = 1) or zero?(n: F := second l) => 1
        (n = 1) => x
        (u := isExpt(x, OPEXP)) case "failed" => kernel(oppow, l)
        rec := u::Record(var: K, exponent: Z)
        ((y := first argument(rec.var))=1) or y = -1 =>
            (operator(rec.var)) (rec.exponent * y * n)
        kernel(oppow, l)

    if R has CombinatorialFunctionCategory then

      iifact x ==
        (r:=retractIfCan(x)@Union(R,"failed")) case "failed" => ifact x
        factorial(r::R)::F

      iiperm l ==
        (r1 := retractIfCan(first l)@Union(R,"failed")) case "failed" or
          (r2 := retractIfCan(second l)@Union(R,"failed")) case "failed"
            => iperm l
        permutation(r1::R, r2::R)::F

      if R has RetractableTo(Z) and F has Algebra(Fraction(Z)) then

         iibinom l ==
           (s:=retractIfCan(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

-- iibinom checks those cases in which the binomial coefficient may
-- be evaluated explicitly. Currently, the naive iterative algorithm is
-- used to calculate the coefficient, there is room for improvement here.

      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])

-- These four properties define special differentiation rules for sums and 
-- products. 

    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)

-- Set the properties for displaying sums and products and testing for
-- equality.

    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)