open-axiom repository from github
\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}