\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{\$SPAD/src/algebra combinat.spad}
\author{Martin Brock, Robert Sutor, Michael Monagan}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{package COMBINAT IntegerCombinatoricFunctions}
<<package COMBINAT IntegerCombinatoricFunctions>>=
)abbrev package COMBINAT IntegerCombinatoricFunctions
++ Authors: Martin Brock, Robert Sutor, Michael Monagan
++ Date Created: June 1987
++ Date Last Updated:
++ Basic Operations:
++ Related Domains:
++ Also See:
++ AMS Classifications:
++ Keywords: integer, combinatoric function
++ Examples:
++ References:
++ Description:
++   The \spadtype{IntegerCombinatoricFunctions} package provides some
++   standard functions in combinatorics.
Z   ==> Integer
N   ==> NonNegativeInteger
SUP ==> SparseUnivariatePolynomial

IntegerCombinatoricFunctions(I:IntegerNumberSystem): with
   binomial: (I, I) -> I
      ++ \spad{binomial(n,r)} returns the binomial coefficient
      ++ \spad{C(n,r) = n!/(r! (n-r)!)}, where \spad{n >= r >= 0}.
      ++ This is the number of combinations of n objects taken r at a time.
   factorial: I -> I
      ++ \spad{factorial(n)} returns \spad{n!}. this is the product of all
      ++ integers between 1 and n (inclusive). 
      ++ Note: \spad{0!} is defined to be 1.
   multinomial: (I, List I) -> I
      ++ \spad{multinomial(n,[m1,m2,...,mk])} returns the multinomial
      ++ coefficient \spad{n!/(m1! m2! ... mk!)}.
   partition: I -> I
      ++ \spad{partition(n)} returns the number of partitions of the integer n.
      ++ This is the number of distinct ways that n can be written as
      ++ a sum of positive integers.
   permutation: (I, I) -> I
      ++ \spad{permutation(n)} returns \spad{!P(n,r) = n!/(n-r)!}. This is
      ++ the number of permutations of n objects taken r at a time.
   stirling1: (I, I) -> I
      ++ \spad{stirling1(n,m)} returns the Stirling number of the first kind
      ++ denoted \spad{S[n,m]}.
   stirling2: (I, I) -> I
      ++ \spad{stirling2(n,m)} returns the Stirling number of the second kind
      ++ denoted \spad{SS[n,m]}.
 == add
   F : Record(Fn:I, Fv:I) := [0,1]
   B : Record(Bn:I, Bm:I, Bv:I) := [0,0,0]
   S : Record(Sn:I, Sp:SUP I) := [0,0]
   P : IndexedFlexibleArray(I,0) := new(1,1)$IndexedFlexibleArray(I,0)
 
   partition n ==
      -- This is the number of ways of expressing n as a sum of positive
      -- integers, without regard to order.  For example partition 5 = 7
      -- since 5 = 1+1+1+1+1 = 1+1+1+2 = 1+2+2 = 1+1+3 = 1+4 = 2+3 = 5 .
      -- Uses O(sqrt n) term recurrence from Abramowitz & Stegun pp. 825
      -- p(n) = sum (-1)**k p(n-j) where 0 < j := (3*k**2+-k) quo 2 <= n
      minIndex(P) ~= 0 => error "Partition: must have minIndex of 0"
      m := #P
      negative? n => error "partition is not defined for negative integers"
      n < m::I => P(convert(n)@Z)
      concat!(P, new((convert(n+1)@Z - m)::N,0)$IndexedFlexibleArray(I,0))
      for i in m..convert(n)@Z repeat
         s: I := 1
         t: I := 0
         for k in 1.. repeat
            l := (3*k*k-k) quo 2
            l > i => leave
            u := l+k
            t := t + s * P(convert(i-l)@Z)
            u > i => leave
            t := t + s * P(convert(i-u)@Z)
            s := -s
         P.i := t
      P(convert(n)@Z)
 
   factorial n ==
      s,f,t : I
      negative? n => error "factorial not defined for negative integers"
      if n <= F.Fn then s := f := 1 else (s, f) := F
      for k in convert(s+1)@Z .. convert(n)@Z by 2 repeat
         if k::I = n then t := n else t := k::I * (k+1)::I
         f := t * f
      F.Fn := n
      F.Fv := f
 
   binomial(n, m) ==
      negative? n or negative? m or m > n => 0
      m = 0 => 1
      n < 2*m => binomial(n, n-m)
      s: I := 0
      b: I := 1
      if B.Bn = n then
         B.Bm = m+1 =>
            b := (B.Bv * (m+1)) quo (n-m)
            B.Bn := n
            B.Bm := m
            return(B.Bv := b)
         if m >= B.Bm then
           s := B.Bm
           b := B.Bv
         else
           s := 0
           b := 1
      for k in convert(s+1)@Z .. convert(m)@Z repeat
        b := (b*(n-k::I+1)) quo k::I
      B.Bn := n
      B.Bm := m
      B.Bv := b
 
   multinomial(n, m) ==
      for t in m repeat negative? t => return 0
      n < _+/m => 0
      s:I := 1
      for t in m repeat s := s * factorial t
      factorial n quo s
 
   permutation(n, m) ==
      negative? m or n < m => 0
      m := n-m
      p:I := 1
      for k in convert(m+1)@Z .. convert(n)@Z by 2 repeat
         t: I :=
           k::I = n => n
           (k*(k+1))::I
         p := p * t
      p
 
   stirling1(n, m) ==
      -- Definition: (-1)**(n-m) S[n,m] is the number of
      -- permutations of n symbols which have m cycles.
      negative? n or m < 1 or m > n => 0
      m = n => 1
      S.Sn = n => coefficient(S.Sp, convert(m)@Z :: N)
      x := monomial(1, 1)$SUP(I)
      S.Sn := n
      S.Sp := x
      for k in 1 .. convert(n-1)@Z repeat S.Sp := S.Sp * (x - k::SUP(I))
      coefficient(S.Sp, convert(m)@Z :: N)
 
   stirling2(n, m) ==
      -- definition: SS[n,m] is the number of ways of partitioning
      -- a set of n elements into m non-empty subsets
      negative? n or m < 1 or m > n => 0
      m = 1 or n = m => 1
      s:I := if odd? m then -1 else 1
      t:I := 0
      for k in 1..convert(m)@Z repeat
         s := -s
         t := t + s * binomial(m, k::I) * k::I ** (convert(n)@Z :: N)
      t quo factorial m

@
\section{License}
<<license>>=
--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
--All rights reserved.
--
--Redistribution and use in source and binary forms, with or without
--modification, are permitted provided that the following conditions are
--met:
--
--    - Redistributions of source code must retain the above copyright
--      notice, this list of conditions and the following disclaimer.
--
--    - Redistributions in binary form must reproduce the above copyright
--      notice, this list of conditions and the following disclaimer in
--      the documentation and/or other materials provided with the
--      distribution.
--
--    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
--      names of its contributors may be used to endorse or promote products
--      derived from this software without specific prior written permission.
--
--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
@
<<*>>=
<<license>>

<<package COMBINAT IntegerCombinatoricFunctions>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}