axiom-test 20170501-3 / usr / share / axiom-20170501 / input / combinatorics.input

)set break resume
)sys rm -f combinatorics.output
)spool combinatorics.output
)set message test on
)set message auto off
)clear all

powerSet(S:Set Any):Set Set Any==
++ This produces the power set (set of all subsets) of a given set, using a 
++ simple recursive algorithm.
  card:=cardinality(S)
  if card=1 then
     return(set [set[],S])
  else
     x:=members(S).card
     S1:=remove(x,S)
     S1P:=powerSet(S1)
     S1Px:Set Set Any:=set [union(set[x],s) for s in members(S1P)]
     return(union(S1P,S1Px))

choose(S:Set ANY, n:NNI):Set Set Any==
++ choose(S,n) produces all subsets of S with cardinality n.  Again, a simple
++ recursive algorithm is used
  card:=cardinality(S)
  if n=0 then
    return(set[set[]])
  if card=n then
    return(set[S])
  else
    x:=members(S).1
    S1:=remove(x,S)
    S1P:=choose(S1,n-1)
    S1Px:Set Set Any:=set [union(set[x],s) for s in members(S1P)]
    return(union(choose(S1,n),S1Px))

nextperm(L:List INT,n:INT):List INT==
++ The driver of Algorithm L: given any permutation, produces the next one
++ in lexicographical order.  This algorithm is implemented in "lexperm" below
  a:=L
  for i in 1..n-1 repeat if a.i<a.(i+1) then k:=i
  for i in k+1..n repeat if a.k<a.i then j:=i
  temp:=a.k
  a.k:=a.j
  a.j:=temp
  next:List INT:=concat([a.i for i in 1..k],reverse([a.i for i in k+1..n]))
  return(next)

lexperm(n:INT):List List INT==
++ This is an implementation of Knuth's "Algorithm L" for listing 
++ permutations of [1..n] in lexicographical order.  
  a:List INT:=[k for k in 1..n]
  P:List List INT:=[]
  repeat
    ++print(a)
    P:=append(P,[copy(a)])
    if #[k for k in 1..n-1 | a.k<a.(k+1)]=0 then
      return(P)
    else
      a:=nextperm(a,n)

lexperm2(L:List INT):List List INT==
++ A variation of lexperm which applies the Algorithm L to any list of integers
++ It works well for lists which include repeated elements
  a:=L
  n:=#L
  P:List List INT:=[]
  repeat
    ++print(a)
    P:=append(P,[copy(a)])
    if #[k for k in 1..n-1 | a.k<a.(k+1)]=0 then
      return(P)
    else
      a:=nextperm(a,n)

listPermutations(L:List Any):List List Any==
++ Another variant of lexperm, this time for any list
  E:=removeDuplicates(L)
  EN:=sort([position(x,E) for x in L])
  P:=lexperm2(EN)
  Q:List List Any:=[[E.i for i in p] for p in P]
  return(Q)

randomPermutation(L:List Any):List Any==
++ This is an implementation of the Knuth Shuffle
  M:=copy(L)
  n:=#L
  for k in 1..n-1 repeat
    p:=random(n)+1
    temp:=M.p
    M.p:=M.k
    M.k:=temp
  return(M)

derangements(n:INT):INT==
++ returns the number of derangements D(n) of n objects, using the recursion 
++ D(1)=0,D(2)=1,D(n)=(n-1)*(D(n-1)+D(n-2))
  if n<3 then
    return(n-1)
  else
    a:=0
    b:=1
    for i in 3..n repeat
      c:=(i-1)*(b+a)
      a:=b
      b:=c
    return(c)

fixed?(L1:List Any,L2:List Any):Boolean==
++ This is a helper program for computing derangements.  It returns true if the
++ input lists have at least one element in the same place in each list. 
  if #L1~=#L2 then
    error("Lists must be of the same length")
  else
    return(reduce(\/, [(L1.i=L2.i)@Boolean for i in 1..#L1]))

listDerangements(L:List Any):List List Any==
++ List all the derangements (permutations with no fixed points) by the 
++ simple method of first listing all permutations, and keeping only 
++ those with no fixed points.
  P:=listPermutations(L)
  Ds:List List Any:=[]
  for X in P repeat
    if ~fixed?(L,X) then
       Ds:=append(Ds,[X])
  return(Ds)

countDerangements(L:List Any):INT==
++ Returns the number of derangements of a list possibly with repeated
++ elements, using McMahon's method of generating functions
  E:=removeDuplicates(L)
  a:=[count(x,L) for x in E]
  if reduce(max,a)=1 then
    return(derangements(#a))
  else
    vs:=[subscript('x::Symbol,[i])@Symbol for i in 1..#a]
    S:=reduce(+,vs)
    P:DMP(vs,FRAC INT):=reduce(*,[(S-vs.i)^a.i for i in 1..#a])
    c:=coefficient(P,vs.1,a.1)
    for i in 2..#a repeat
      c:=coefficient(c,vs.i,a.i)
    return(c)

randomDerangement(L:List Any):List Any==
++ Produces a random derangment by choosing permutations at random until 
++ one of them is a derangement.  This is not in fact guaranteed to work, 
++ but the probability of it failing is negligible.  This should be replaced 
++ with a better program.
  if countDerangements(L)=0 then
    error("The list has zero derangements")
  else
    repeat
      X:=randomPermutation(L)
      if ~fixed?(X,L) then
        return(X)

listStringPermutations(S:String):List String==
  E:=removeDuplicates(S)
  L:=sort([position(x,parts(E)) for x in parts(S)])
  P:=lexperm2(L)
  Q:List String:=[reduce(concat,[(E.x)::String for x in p]) for p in P]
  return(Q)

stringFixed?(S1:String,S2:String):Boolean==
  if (#S1~=#S2) then
     error("Strings must have the same length")
  else
     return(reduce(\/, [(S1.i=S2.i)@Boolean for i in 1..#S1]))

listStringDerangements(S:String):List String==
  P:=listStringPermutations(S)
  Ds:List String:=[]
  for X in P repeat
    if ~stringFixed?(S,X) then
       Ds:=append(Ds,[X])
  return(Ds)

listPartitions(n:INT):List List INT==
++ This list all the set partitions in terms of "codes": for a set S with 
++ n elements a partition code C is a list of n numbers 1,2,..,k with 
++ k<=n, such that the i-th element of S belongs to subset C.i.  For example,
++ if S={2,4,6,8,10} and C=[1,2,1,3,1] then the corresponding subset 
++ partition is {{2,6,10},{4},{8}}
  p:List List INT:=[[1]]
  if n=1 then
    return(p)
  else
    for i in 2..n repeat
      q:List List INT:=[]
      for x in p repeat
        c:=reduce(max,x)+1
        np:=[append(x,[j]) for j in 1..c]
        q:=concat(q,np)
      p:=copy(q)
    return(p)

setPartitions(S:Set Any):Set Set Any==
++ For a given set S, this returns the subset list generated from the codes 
++ produced by the listPartitions function
  n:=#S
  P:List List INT:=listPartitions(n)
  Q:Set Set Any:=set[]
  for x in P repeat
    R:List Set Any:=[set[] for i in 1..reduce(max,x)]
    for i in 1..n repeat R.(x.i):=union(R.(x.i),set[parts(S).i])
    Q:=union(Q,set R)
  return(Q)

listSizePartitions(n:INT,k:INT):List List INT==
++ This lists the codes of all the set partitions whose maximum value is k.  
++ These codes correspond to partitions of exactly k subsets.
  p:List List INT:=[[1]]
  if n=1 then
    return(p)
  else
    for i in 2..n repeat
      q:List List INT:=[]
      for x in p repeat
        c:=reduce(max,x)+1
        np:=[append(x,[j]) for j in 1..min(k,c)]
        q:=concat(q,np)
      p:=copy(q)
    for x in p repeat
      if reduce(max,x)<k then p:=remove(x,p)
    return(p)

setSizePartitions(S:Set Any,k:INT):Set Set Any==
++ For a given set S, this returns all partitions of S into exactly k 
++ disjoint subsets
  n:=#S
  P:List List INT:=listSizePartitions(n,k)
  Q:Set Set Any:=set[]
  for x in P repeat
    R:List Set Any:=[set[] for i in 1..k]
    for i in 1..n repeat R.(x.i):=union(R.(x.i),set[parts(S).i])
    Q:=union(Q,set R)
  return(Q)

bell(n:NNI):PI==
++ The Bell numbers B(n): the number of ways a set of n elements can be 
++ partitioned into disjoint subsets.  Bell numbers are in fact sums of 
++ Stirling numbers of the second kind so that we could simply define
++ bell(n)==reduce(+,[stirling2(n,k) for k in 0..n])
  if n<2 then
    return(1)
  else
    B:=[1,1]
    for i in 2..n repeat
      N:=reduce(+,[binomial(i-1,k)*B.(k+1) for k in 0..i-1])
      B:=append(B,[N])
    return(B.last)

touchard(n:PI,z:Symbol):POLY FRAC INT==
++ The Touchard polynomials T(n,x) are generating functions for Stirling 
++ numbers of the second kind: the coefficient of x^k in T(n,x) is the 
++ value stirling2(n,k)
  if n=0 then
    return(1)
  else
    Tc:POLY FRAC INT:=1
    for i in 1..n repeat
      Tc:=z*(Tc+D(Tc,z))
    return(Tc)

randomStringPermutation(S:String):String==
  L:List String:=parts(S)
  P:=randomPermutation(L)
  return(reduce(concat,P))

randomStringDerangement(S:String):String==
  E:=removeDuplicates(S)
  a:=[count(x,S) for x in parts(E)]
  if reduce(\/,[(2*z>#S)@Boolean for z in a]) then
    error("The string has zero derangements")
  else
    repeat
      X:=randomStringPermutation(S)
      if ~stringFixed?(X,S) then
        return(X)

--S 1 of 21
S := set[4,14,46,5]
--R 
--R
--R   (24)  {4,5,14,46}
--R                                                   Type: Set(PositiveInteger)
--E 1

--S 2 of 21
powerSet(S)
--R 
--R   Compiling function powerSet with type Set(Any) -> Set(Set(Any)) 
--R
--R   (25)
--R   {{}, {4}, {5}, {5,4}, {14}, {14,4}, {14,5}, {14,5,4}, {46}, {46,4}, {46,5},
--R    {46,5,4}, {46,14}, {46,14,4}, {46,14,5}, {46,14,5,4}}
--R                                                          Type: Set(Set(Any))
--E 2

--S 3 of 21
S:=set["A",2.5,vector[1,2,3]]
--R 
--R
--R   (26)  {"A",2.5,[1,2,3]}
--R                                                               Type: Set(Any)
--E 3

--S 4 of 21
powerSet(S)
--R 
--R
--R   (27)
--R   {{}, {"A"}, {2.5}, {2.5,"A"}, {[1,2,3]}, {[1,2,3],"A"}, {[1,2,3],2.5},
--R    {[1,2,3],2.5,"A"}}
--R                                                          Type: Set(Set(Any))
--E 4

--S 5 of 21
S:=set["cat","dog","fly","eel"]
--R 
--R
--R   (28)  {"cat","dog","eel","fly"}
--R                                                            Type: Set(String)
--E 5

--S 6 of 21
choose(S,2)
--R 
--R   There are no library operations named choose
--R      Use HyperDoc Browse or issue
--R                               )what op choose
--R      to learn if there is any operation containing choose in its name.
--R   Cannot find a definition or applicable library operation named 
--R      choose with argument type(s)
--R                                  Set(Any)
--R                                   Integer
--R      Perhaps you should use @ to indicate the required return type, or
--R      $ to specify which version of the function you need.
--R   Axiom will attempt to step through and interpret the code.
--R   Compiling function choose with type (Set(Any),NonNegativeInteger)
--R       -> Set(Set(Any)) 
--I   Compiling function G2391 with type Integer -> Boolean 
--R
--R   (29)
--R   {{"eel","fly"}, {"dog","fly"}, {"dog","eel"}, {"cat","fly"}, {"cat","eel"},
--R    {"cat","dog"}}
--R                                                          Type: Set(Set(Any))
--E 6

--S 7 of 21
L:=[1,1,2,3,3]
--R 
--R
--R   (30)  [1,1,2,3,3]
--R                                                  Type: List(PositiveInteger)
--E 7

--S 8 of 21
listPermutations(L)
--R 
--R   Compiling function nextperm with type (List(Integer),Integer) -> 
--R      List(Integer) 
--R   Compiling function lexperm2 with type List(Integer) -> List(List(
--R      Integer)) 
--R   Compiling function listPermutations with type List(Any) -> List(List
--R      (Any)) 
--R
--R   (31)
--R   [[1,1,2,3,3], [1,1,3,2,3], [1,1,3,3,2], [1,2,1,3,3], [1,2,3,1,3],
--R    [1,2,3,3,1], [1,3,1,2,3], [1,3,1,3,2], [1,3,2,1,3], [1,3,2,3,1],
--R    [1,3,3,1,2], [1,3,3,2,1], [2,1,1,3,3], [2,1,3,1,3], [2,1,3,3,1],
--R    [2,3,1,1,3], [2,3,1,3,1], [2,3,3,1,1], [3,1,1,2,3], [3,1,1,3,2],
--R    [3,1,2,1,3], [3,1,2,3,1], [3,1,3,1,2], [3,1,3,2,1], [3,2,1,1,3],
--R    [3,2,1,3,1], [3,2,3,1,1], [3,3,1,1,2], [3,3,1,2,1], [3,3,2,1,1]]
--R                                                        Type: List(List(Any))
--E 8

--S 9 of 21
listDerangements(L)
--R 
--R   Compiling function fixed? with type (List(Any),List(Any)) -> Boolean
--R      
--R   Compiling function listDerangements with type List(Any) -> List(List
--R      (Any)) 
--R
--R   (32)  [[2,3,3,1,1],[3,2,3,1,1],[3,3,1,1,2],[3,3,1,2,1]]
--R                                                        Type: List(List(Any))
--E 9

--S 10 of 21
listStringPermutations("EERIE")
--R 
--R   Compiling function listStringPermutations with type String -> List(
--R      String) 
--R
--R   (33)
--R   ["EEERI", "EEEIR", "EEREI", "EERIE", "EEIER", "EEIRE", "EREEI", "EREIE",
--R    "ERIEE", "EIEER", "EIERE", "EIREE", "REEEI", "REEIE", "REIEE", "RIEEE",
--R    "IEEER", "IEERE", "IEREE", "IREEE"]
--R                                                           Type: List(String)
--E 10

--S 11 of 21
listStringDerangements("banana")
--R 
--R   Compiling function stringFixed? with type (String,String) -> Boolean
--R      
--R   Compiling function listStringDerangements with type String -> List(
--R      String) 
--R
--R   (34)  ["abanan","anaban","ananab"]
--R                                                           Type: List(String)
--E 11

--S 12 of 21
stirling2(3,2)
--R 
--R
--R   (35)  3
--R                                                        Type: PositiveInteger
--E 12

--S 13 of 21
listPartitions(3)
--R 
--R   Compiling function listPartitions with type Integer -> List(List(
--R      Integer)) 
--R
--R   (36)  [[1,1,1],[1,1,2],[1,2,1],[1,2,2],[1,2,3]]
--R                                                    Type: List(List(Integer))
--E 13

--S 14 of 21
G:Set Any:=set[2.5,'x+'y,"A"]
--R 
--R
--R   (37)  {2.5,y + x,"A"}
--R                                                               Type: Set(Any)
--E 14

--S 15 of 21
setPartitions(G)
--R 
--R   Compiling function setPartitions with type Set(Any) -> Set(Set(Any))
--R      
--R
--R   (38)
--R   {{{2.5,y + x,"A"}}, {{2.5,y + x},{"A"}}, {{2.5,"A"},{y + x}},
--R    {{2.5},{y + x,"A"}}, {{2.5},{y + x},{"A"}}}
--R                                                          Type: Set(Set(Any))
--E 15

--S 16 of 21
P:=listSizePartitions(5,3)
--R   Compiling function listSizePartitions with type (Integer,Integer)
--R       -> List(List(Integer)) 
--R
--R   (39)
--R   [[1,1,1,2,3], [1,1,2,1,3], [1,1,2,2,3], [1,1,2,3,1], [1,1,2,3,2],
--R    [1,1,2,3,3], [1,2,1,1,3], [1,2,1,2,3], [1,2,1,3,1], [1,2,1,3,2],
--R    [1,2,1,3,3], [1,2,2,1,3], [1,2,2,2,3], [1,2,2,3,1], [1,2,2,3,2],
--R    [1,2,2,3,3], [1,2,3,1,1], [1,2,3,1,2], [1,2,3,1,3], [1,2,3,2,1],
--R    [1,2,3,2,2], [1,2,3,2,3], [1,2,3,3,1], [1,2,3,3,2], [1,2,3,3,3]]
--R                                                    Type: List(List(Integer))
--E 16

--S 17 of 21
#P
--R 
--R
--R   (40)  25
--R                                                        Type: PositiveInteger
--E 17

--S 18 of 21
stirling2(5,3)
--R 
--R
--R   (41)  25
--R                                                        Type: PositiveInteger
--E 18

--S 19 of 21
p := touchard(7,x)
--R 
--R   Compiling function touchard with type (PositiveInteger,Symbol) -> 
--R      Polynomial(Fraction(Integer)) 
--R
--R          7      6       5       4       3      2
--R   (42)  x  + 21x  + 140x  + 350x  + 301x  + 63x  + x
--R                                          Type: Polynomial(Fraction(Integer))
--E 19

--S 20 of 21
[coefficient(p,x,k) for k in 0..7]
--R 
--R
--R   (43)  [0,1,63,301,350,140,21,1]
--R                                    Type: List(Polynomial(Fraction(Integer)))
--E 20

--S 21 of 21
[stirling2(7,k) for k in 0..7]
--R 
--R
--R   (44)  [0,1,63,301,350,140,21,1]
--R                                                          Type: List(Integer)
--E 21


)spool