/usr/share/axiom-20170501/src/algebra/COMBINAT.spad is in axiom-source 20170501-3.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 | )abbrev package COMBINAT IntegerCombinatoricFunctions
++ Authors: Martin Brock, Robert Sutor, Michael Monagan
++ Date Created: June 1987
++ Description:
++ The \spadtype{IntegerCombinatoricFunctions} package provides some
++ standard functions in combinatorics.
IntegerCombinatoricFunctions(I) : SIG == CODE where
I:IntegerNumberSystem
Z ==> Integer
N ==> NonNegativeInteger
SUP ==> SparseUnivariatePolynomial
SIG ==> 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.
++
++X [binomial(5,i) for i in 0..5]
factorial : I -> I
++ \spad{factorial(n)} returns \spad{n!}. this is the product of all
++ integers between 1 and n (inclusive).
++ Note that \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]}.
CODE ==> 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
n < 0 => 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
n < 0 => 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) ==
s,b:I
n < 0 or m < 0 or m > n => 0
m = 0 => 1
n < 2*m => binomial(n, n-m)
(s,b) := (0,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,b) := (0,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 t < 0 => return 0
n < _+/m => 0
s:I := 1
for t in m repeat s := s * factorial t
factorial n quo s
permutation(n, m) ==
t:I
m < 0 or n < m => 0
m := n-m
p:I := 1
for k in convert(m+1)@Z .. convert(n)@Z by 2 repeat
if k::I = n then t := n else t := (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.
n < 0 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
n < 0 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
|