#!/usr/bin/env maple

interface(quiet=true):
# Maple program to produce the Table II
# of Rossky & Karplus' enumeration of Goldstone diagrams
# with n time labels and k external potentials.
# Richard J. Mathar, 2019-10-30


# Generate lists  compatible with (16)(1) 2nd equation
# return [[a1,a2,...],[a1,a2,...]]
Alists := proc(n::integer)
	option remember;
	local ai,jsum,p,ais,aiss;
	ai := [] ;
	# j <=n so the sum of the ai<=n and the compositions sum <=2*n
	for jsum from 0 to 2*n do
		for p in combinat[composition](jsum,n-1) do
			ais := [seq(op(i,p)-1,i=1..nops(p))] ;
			aiss := add(i*op(i,ais),i=1..n-1) ;
			if aiss = n then
				ai := [op(ai),ais] ;
			end if;
		end do:
	end do:
	ai ;
end proc:

# given an index i and constant ai, produce lists [b0i,b1i,...bii]
bisets := proc(i::integer,ai::integer)
	local bi,bis,p ;
	bi := [] ;
	for p in combinat[composition](ai+i+1,i+1) do
		bis := [seq(op(j,p)-1,j=1..nops(p))] ;
		bi := [op(bi),bis] ;
	end do:
	bi ;
end proc:

# Argument is a [a1,a2,...] list of ai.
# Return is [ [[[b01,[b11],],[b01,b11],,..], [[b02,b12,b22,],[b02,b12,b22],..]
# so the i-th component is a list of possible lists [b0i, b1i, b2i,..bii]
Bilists := proc(ailist::list)
	option remember;
	local bil,ai,n,i,thisbi ;
	bil := [] ;
	for i from 1 to nops(ailist) do
		ai := op(i,ailist) ;
		thisbi := bisets(i,ai) ;
		bil := [op(bil),thisbi] ;
	end do:
	bil ;
end proc:

Expabi := proc(bicols::list)
	local fin,pp,l,isum,ii,vlamx,i ;
	fin := [] ;
	l := nops(bicols) ;
	for isum from l to add(nops(op(i,bicols)),i=1..l) do
		for ii in combinat[composition](isum,l) do
			vlamx := true ;
			for i from 1 to l do
				if op(i,ii) > nops(op(i,bicols)) then
					vlamx := false;
					break ;
				end if;
			end do:
			if vlamx then
				fin := [op(fin), [seq( op(op(i,ii),op(i,bicols)),i=1..l)]] ;
			end if;
		end do:
	end do:
	fin ;
end proc:

# values C(n,k) of connected diagrams, n>=0, 0<=k<=n
# A328921
CKarplus := proc(n::integer,k::integer)
	option remember;
	local ailists,resu,j,ailist,bcols,bmats,bflat,thisk,lam,f,bli,i ;
	if k <0 or k > n then
		return 0 ;
	elif n = 0 then
		return 1 ;
	end if;
	resu := 0 ;
	ailists := Alists(n) ;
	for j from 2 to n do
		for ailist in ailists do
			if add(a,a=ailist) = j then
				bcols := Bilists(ailist) ;
				# build all possible lines of [b00,b01,b02,b12,b22,...]
				# by selecting one element of the first, one of the 2nd,
				# one of the third,.. lists of bcols.
				bmats := Expabi(bcols) ;
				# loop over the rows of the {{b lambda i}} of Table 1
				for bflat in bmats do
					thisk := 0 ;
					for i from 1 to nops(bflat) do
						thisk := thisk+ add(lam*op(lam+1,op(i,bflat)),lam=0..i) ;
					end do:
					if thisk = k then
						f := 1 ;
						for i from 1 to n-1 do
							f := f/(i!)^op(i,ailist) ;
							for lam from 0 to i do
								bli := op(lam+1,op(i,bflat)) ;
								f := f*procname(i,lam)^bli/bli! ;
							end do:
						end do:
						resu := resu+f ;
					end if ;
				end do:
			end if;
		end do:
	end do:
	# equation (18)
	binomial(n,k)*(2*n-k)!-n!*resu ;
end proc:


# D(n,k)
# A328922
DKarplus := proc(n::integer,k::integer)
	binomial(n,k)*(2*n-k)!-CKarplus(n,k) ;
end proc:

# T(n,k)
# A328923
TKarplus := proc(n::integer,k::integer)
	if n = 0 then
		1;
	elif k <> 0 then
		2^(k-n)*CKarplus(n,k) ;
	else
		CKarplus(n,k)/2^n+(n-1)! ;
	fi;
end proc:

# Teff(n,k)
# A328924
TeffKarplus := proc(n::integer,k::integer)
	add((-1)^(k+kpr)*2^(kpr-k)*binomial(kpr,k)*TKarplus(n,kpr),kpr=k..n) ;
end proc:

# print C(n,k) of Table II
for n from 0 to 8 do
	for k from 0 to n do
		printf("%11d ",CKarplus(n,k)) ;
	end do:
	printf("\n") ;
end do:

# print D(n,k) of Table II
for n from 0 to 8 do
	for k from 0 to n do
		printf("%10d ",DKarplus(n,k)) ;
	end do:
	printf("\n") ;
end do:

# print T(n,k) of Table II
for n from 0 to 8 do
	for k from 0 to n do
		printf("%10d ",TKarplus(n,k)) ;
	end do:
	printf("\n") ;
end do:

# print Teff(n,k) of Table II
for n from 0 to 8 do
	for k from 0 to n do
		printf("%10d ",TeffKarplus(n,k)) ;
	end do:
	printf("\n") ;
end do: