# The partition-product and combinatorial triangles. # Author: Peter Luschny, written 2009-02-02, Maple V R5 # More information at: # www.luschny.de/math/seq/CountingWithPartitions.html restart; readlib(queue): # Maple puts an constant strain on the user with # inadequate definitions of the power and the # pochhammer function. So let's first correct this, # in accordance to Knuth's 'Two notes on notation'. pow := proc(n,k) if n=0 and k=0 then 1 else n^k fi end: # -- Generalized rising factorial GeneralizedRisingFactorial := proc(x,m,s) local j; if m >= 0 then mul((x+s*j),j=0..m-1) else 1/mul((x+s*j),j=m..-1) fi end: GeneratePartitions := proc(n, visit) local a,b,nextLeft,nextRight,count,normalize, partitions,putPart,nextPart,least,length,next; putPart := (x) -> queue[enqueue](partitions, x); nextPart := () -> queue[dequeue](partitions); least := a -> op(2,a); length := a -> op(3,a); normalize := proc(a) local s,t,i; t:=[]; s:=op(1,a); for i from 1 to length(a) do t := [op(t),s[i]] od; t end; next := proc(b) count := count+1; visit(normalize(b)); if least(b) <= length(b) then putPart(b) fi; end; nextLeft := proc() local u,la; la := least(a); if (la = length(a)) or (a[1][la] > 1) then RETURN(NULL) fi; u := normalize(a); u[la] := u[la]+1; [u,la+1,length(a)-1]; end; nextRight := proc() local u,la; la := least(a); if (la = 1) or ((la <> 2) and (a[1][la-1] = a[1][la-2])) then RETURN(NULL) fi; u := normalize(a); u[la-1] := u[la-1]+1; [u,la,length(a)-1]; end; partitions := queue[new](); count := 0; next([[seq(1,i=1..n)],1,n]); while not queue[empty](partitions) do a := nextPart(); b := nextLeft(); if b <> NULL then next(b); fi; b := nextRight(); if b <> NULL then next(b); fi; od; count end: PrintPartitionTree := proc(n) local R,len,visit; visit := proc(L) local i; if nops(L) = len then R := [op(R), L] else len := len - 1; print(R); R := [L] fi; end; R := []; len := n; GeneratePartitions(n,visit); print(R); end: for i from 3 to 7 do print("--",i,"--"); PrintPartitionTree(i); od: PartitionProduct := proc(n,f,g,statistic) local j,R,visit; visit := proc(L) local i,A,S,U,V,W; A := add(pow(x,L[i]), i=1..nops(L)); S := [seq(coeff(A,x,i), i=1..n)]; V := mul(pow(f(i),S[i]),i=1..n); W := mul(pow(g(i),S[i])*S[i]!,i=1..n); U := abs(n!*V/W); if statistic = "sum" then R := R + U elif statistic = "part" then R := [op(R), U] elif statistic = "len" then R[nops(L)] := R[nops(L)] + U elif statistic = "big" then R[L[1]] := R[L[1]] + U ; fi; end; if n = 0 then if statistic = "sum" then RETURN(1) else RETURN([1]) fi fi; if statistic = "sum" then R := 0 elif statistic = "part" then R := [] else R := [seq(0,j=1..n)] fi; GeneratePartitions(n,visit); R end: Stirling1Triangles := proc(n,k,statistic) local f; f := proc(n) GeneralizedRisingFactorial(k-n+2,n-1,1) end: PartitionProduct(n,f,factorial,statistic) end: Stirling2Triangles := proc(n,k,statistic) local f; f := proc(n) GeneralizedRisingFactorial(-1,n,k+1) end: PartitionProduct(n,f,factorial,statistic) end: BesselTriangles := proc(n,k,statistic) local f; f := proc(n) GeneralizedRisingFactorial(k-n+2,n-1,1) end: PartitionProduct(n,f,n->n,statistic) end: PrintTriangles := proc(T,k,len) local i,n; print(seq(T(n,k,"sum"),n=0..len)); seq(print(T(n,k,"part")),n=0..len); seq(print(T(n,k,"len")),n=0..len); seq(print(T(n,k,"big")),n=0..len); end: PrintForOEIS := proc(T,k,len,statistic) local n; seq(op(T(n,k,statistic)),n=0..len) end: for k from -6 to 5 do print("Stirling1Triangle",k); PrintTriangles(Stirling1Triangles,k,5) od; for k from -6 to 5 do print("Stirling2Triangle",k); PrintTriangles(Stirling2Triangles,k,5) od; BesselPartitions := proc(n,k) local f; f := proc(n) GeneralizedRisingFactorial(k-n+2,n-1,1) end: PartitionProduct(n,f,n->n,"part") end: for n from 0 to 7 do print(BesselPartitions(n,1)) od; for k from -6 to 5 do print("BesselTriangle",k); PrintTriangles(BesselTriangles,k,5) od;