######################################################################
## IntParts.txt: Save this file as  IntParts.txt to use it.          #
## In Maple, make sure that your working directory is where you      #
## saved the file.                                                   #
## Then type:  read `IntParts.txt` <Enter>                           #
## Then follow the instructions given there                          #
##                                                                   #
## Written by AJ Bu and Robert Dougherty-Bliss, Rutgers University   #
######################################################################

print(`IntParts.txt: A Maple package for statistical analysis of integer partitions`):
print(`First written: May 2021`):
print(`Current Version: May 2021`):
print(`This Maple package accompanies AJ Bu's and Robert Dougherty-Bliss's article`):
print(`A Statistical Analysis of Integer Partitions (fix name later)`):
print(``):
print(`The current version is available at:`):
print(`sites.math.rutgers.edu/~ab1854/Maple%20Packages/IntParts.txt`):
print(`For a list of the MAIN functions, type`):
print(`"Help();"`):
print(` For information on a specific function type`):
print(`"WhatIs(procedure_name);" `):
print(`For a list of the supporting functions type: `):
print(`"Help1():`):

Help:=proc():
	if args=NULL then
		print(`EstStatAnalk(n,k,K,N,rv,m,C,DI)`, `EstStatAnal(n,K,N,rv,m,C,DI)`):
		print(`partitionGFk(f, n, t, m, C, DI)`, `partitionGF(f, n, t, m, C, DI)`, 
			`partitionPGFk(f, n, t, k, m, C, DI)`, `partitionPGF(f, n, t, m, C, DI)`):
		print(`Sgf(k,t)`,  `NiceSeq(k,N)`, `Qgf(k, m, t)`, `Q(n, k, m, t)`):
		print(`QExp(n, k, m)`, `QMom(n, k, m)`, `QClosed(n, K, m)`, `f_PnkRest(k,m,C,DI,x)`):
		print(`For defined options for the r.v., type`):
		print(`"Help(rv)"`):
	elif args=rv then
		print(`numberParts for number of parts`):
		print(`largestPart for largest part`):
		print(`smallestPart for smallest part`):
		print(`distinctParts for number of distinct parts`):
	fi:
end:

Help1:=proc():
	print(`RandRestrictedPnk(n, k, m, C, DI)`, `RandRestrictedPn(n, m, C, DI)`):
	print(`numberParts(pi)`, `largestPart(pi)`, `smallestPart(pi)`, `distinctParts(pi)`):
	print(`NuRestrictedPnk(n, k, m, C, DI)`, `NuRestrictedPn(n, m, C, DI)`):
	print(`restrictedPnk(n, k, m, C, DI)`, `restrictedPn(n, m, C, DI)`):
	print(`LD(L)`, `Ave(X)`, `Mom(X,k)`, `SA(X,K)`):
	print( `Zgf(k,t)`, `powerDiff(f,q,m)`,`Child_PnkRest(k,m,C,DI,x)`):
end:

WhatIs:=proc(procedure_name):
	if procedure_name=EstStatAnalk then
		print(`EstStatAnalk(n,k,K,N,rv,m,C,DI)`):
		print(`estimates the expectation, variance, and the 3rd- through K-th scaled moments`):
		print(`of the random variable rv, of partitions of n with largest part k`):
		print(`s.t. all parts are congruent to a member of C mod m`):
		print(`and the difference between two consecutive parts is NOT in DI`):
		print(`by generating, uniformly at random, a partition with the desired conditions`):
		print(`N times (N should be large)`):
	elif procedure_name=EstStatAnal then
		print(`EstStatAnal(n,K,N,rv,m,C,DI)`):
		print(`estimates the expectation, variance, and the 3rd- through K-th scaled moments`):
		print(`of the random variable rv, of partitions of n `):
		print(`s.t. all parts are congruent to a member of C mod m`):
		print(`and the difference between two consecutive parts is NOT in DI`):
		print(`by generating, uniformly at random, a partition with the desired conditions`):
		print(`N times (N should be large)`):
	elif procedure_name=RandRestrictedPnk then
		print(`RandRestrictedPnk(n, k, m, C, DI)`):
		print(`Generates, uniformly at random, a partition of n with largest part k`):
		print(`s.t. all parts are congruent to a member of C mod m`):
		print(`and the difference between two consecutive parts is NOT in DI`):
	elif procedure_name=RandRestrictedPn then
		print(`RandRestrictedPn(n, m, C, DI)`):
		print(`Generates, uniformly at random, a partition of n`):
		print(`s.t. all parts are congruent to a member of C mod m`):
		print(`and the difference between two consecutive parts is NOT in DI`):
	elif procedure_name=numberParts then
		print(`numberParts(pi)`):
		print(`Inputs an integer partition (as a decreasing list of integers)`):
		print(`Outputs the number of parts`):
	elif procedure_name=largestPart then
		print(`largestPart(pi)`):
		print(`Inputs an integer partition (as a decreasing list of integers)`):
		print(`Outputs the largest part`):	
	elif procedure_name=smallestPart then
		print(`smallestPart(pi)`):
		print(`Inputs an integer partition (as a decreasing list of integers)`):
		print(`Outputs the smallest part`):	
	elif procedure_name=distinctParts then
		print(`distinctParts(pi)`):
		print(`Inputs an integer partition (as a decreasing list of integers)`):
		print(`Outputs the number of distinct parts`):	
	elif procedure_name=NuRestrictedPnk then
		print(`NuRestrictedPnk(n, k, m, C, DI)`):
		print(`The number of partitions of n with largest part k `):
		print(`s.t. all parts are congruent to a member of C mod m`):
		print(`and the difference between two consecutive parts is NOT in DI`):
	elif procedure_name=NuRestrictedPn then
		print(`NuRestrictedPn(n, m, C, DI)`):
		print(`The number of partitions of n`):
		print(`s.t. all parts are congruent to a member of C mod m`):
		print(`and the difference between two consecutive parts is NOT in DI`):
	elif procedure_name=Ave then
		print(`Ave(X)`):
		print(`inputs a discrete random variable X (given as  a list of numbers)`):
		print(`outputs the average (alias expectation, alias mean)`):
	elif procedure_name=Mom then		
		print(`Mom(X,k)`):
		print(`inputs a discrete random variable X (given as  a list of numbers) and a positive integer k`):
		print(`outputs the k-th moment about the mean`):
	elif procedure_name=SA then	
		print(`SA(X,K)`):
		print(`inputs a discrete random variable X (given as  a list of numbers) and a positive integer K`):
		print(`outputs a list of length K whose entries are`):
		print(`the expectation, variance, and the 3rd- through K-th scaled moments`):
	elif procedure_name=LD then
		print(`LD(L)`):
		print(`Inputs a list of positive integers L (of n:=nops(L) members)`):
		print(`outputs an integer i from 1 to n with the prob. of i being proportional to L[i]`):
	elif procedure_name=restrictedPnk then
		print(`restrictedPnk(n,k,m,C,DI)`):
		print(`The set of all partitions of n with largest part k `):
		print(`s.t. all parts are congruent to a member of C mod m`):
		print(`and the difference between two consecutive parts is NOT in DI`):
	elif procedure_name=restrictedPn then
		print(`restrictedPn(n,m,C,DI)`):
		print(`The set of all partitions of n `):
		print(`s.t. all parts are congruent to a member of C mod m`):
		print(`and the difference between two consecutive parts is NOT in DI`):
	elif procedure_name=partitionGFk then
		print(`partitionGFk(f, n, t, k, m, C, DI)`):
		print(`The generating function of the random variable f defined over`):
		print(`integer partitions of n with largest part k`):
		print(`s.t. all parts are congruent to a member of C mod m`):
		print(`and the difference between two consecutive parts is NOT in DI`):
		print(`i.e. sum(t^f(p),p in restrictedPnk(n,k,m,C,DI))`):
	elif procedure_name=partitionGF then
		print(`partitionGF(f, n, t, m, C, DI)`):
		print(`The generating function of the random variable f defined over`):
		print(`integer partitions of n`):
		print(`s.t. all parts are congruent to a member of C mod m`):
		print(`and the difference between two consecutive parts is NOT in DI`):
		print(`i.e. sum(t^f(p),p in restrictedPn(n,m,C,DI))`):
		print(`Try:`):
		print(`partitionGF(p -> 0, 10, t, 1, {0}, {})`):
		print(`partitionGF(largestPart, 10, t, 1, {0}, {})`):
		pritn(`partitionGF(smallestPart, 10, t, 1, {0}, {})`):
		print(`partitionGF(numberParts, 10, t, 1, {0}, {})`):
	elif procedure_name=partitionPGFk then
		print(`partitionPGFk(f, n, t, k, m, C, DI):`):
		print(`The probability generating function of the random variable f defined over`):
		print(`integer partitions of n with largest part k`):
		print(`s.t. all parts are congruent to a member of C mod m`):
		print(`and the difference between two consecutive parts is NOT in DI`):
	elif procedure_name=partitionPGF then
		print(`partitionPGF(f, n, t, m, C, DI):`):
		print(`The probability generating function of the random variable f defined over`):
		print(`integer partitions of n`):
		print(`s.t. all parts are congruent to a member of C mod m`):
		print(`and the difference between two consecutive parts is NOT in DI`):
	elif procedure_name=Zgf then
		print(`Zgf(k,t): The generating function in variable t of the sequence `):
		print(`a_n = total number of parts in all integer partitions of n with largest part at most k`):
		print(`i.e.  sum( Z(n,k)t^n, t=0..infinity), where`):
		print(`Z(n, k) = sum( sum(number of parts of p, p in P(n, k1)), k1=1..k)`):
	elif procedure_name=Sgf then
		print(`Sgf(k,t): The generating function in variable t of the sequence `):
		print(`a_n = total number of parts in all integer partitions of n with largest part k`):
		print(`i.e.  sum( S(n,k)t^n, t=0..infinity), where`):
		print(`S(n, k) = sum(number of parts of p, p in P(n, k))`):
	elif procedure_name=NiceSeq then
		print(`NiceSeq(k,N): Uses the first N terms of the Taylor series expansion of the`):
		print(`generating function, and tries to find a "nice" sequence for the total number`):
		print(`of parts in all integer partitions of n with largest part k by separating cases for`):
		print(`n mod lcm(1,...,k).`):
		print(`It outputs a list where the i-th entry is the total number of parts for partitions of`):
		print(`lcm(1,...,k)*m+i`):
	elif procedure_name=powerDiff then
		print(`powerDiff(f, q, m): (qD_q)^m f, where D_q is the differentiation operator with respect to q.`):
	elif procedure_name=Qgf then
		print(`Qgf(k, m, t): The generating function in t with respect to n of the sequence`):
		print(`Q(n, k, m) = sum(L(p)^m, p), where the sum extends over all partitions of n`):
		print(`with largest part AT MOST k.`):
	elif procedure_name=Q then
		print(`Q(n, k, m, t): sum(L(p)^m, p), where the sum extends over all partitions of n`):
		print(`with largest part AT MOST k.`):
		print(`Try:`):
		print(`Q(10, 2, 1, t)`):
	elif procedure_name=QExp then
		print(`QExp(n, k, m): The mth moment about 0 of the random variable "number of parts of a random partition of n with largest part at most k."`):
		print(`Try:`):
		print(`QExp(n, k, 1) should be roughly n * harmonic(k) / k.`):
		print(`evalf(seq(QExp(n, 2, 1) / n, n=1..20));`):
		print(`evalf(seq(QExp(n, 2, 2) / n^2, n=1..20));`):
	elif procedure_name=QMom then
		print(`QMom(n, k, m): the mth moment about the mean of the random variable "number of parts of a random partition of n with largest part at most k."`):
		print(`Try:`):
		print(`evalf(seq(QMom(n, 2, 1), n=1..20));`):
		print(`evalf(seq(QMom(n, 2, 2)/ n^2, n=1..20));`):
		print(`evalf(seq(QMom(n, 2, 3), n=1..20));`):
	elif procedure_name=QClosed then
		print(`QClosed(n, K, m): The closed form of Q(n,k,m)`):
	elif procedure_name=Child_PnkRest then
		print( `Child_PnkRest(k,m,C,DI,x): inputs positive integers k and m, sets C and DI, and variable x `):
		print(`Outputs (1) the equation relating Z[k] to its "children"`, `and (2) its "children"`):
		print(`where Z[k] is the set of partitions containing the empty partition and partitions where the largest part`):
		print(`is k, all parts mod m are congruent to an element of C, and the difference between any two terms is not in`):
		print(`the set DI`):
	elif procedure_name=f_PnkRest then
		print(`f_PnkRest(k,m,C,DI,x): inputs positive integers k and m, sets C and DI, and variable x`):
		print(`Outputs the generating function of the set of integer where the largest part  is k `):
		print(`all parts mod m are congruent to an element of C, and the difference between any two terms is not in`):
		print(`the set DI`):
	fi:
end:

with(gfun):

############################	General Partition Procedures	############################

#restrictedPnk(n,k,m,C,DI): The set of all partitions of n with largest part k
# s.t. all parts are congruent to a member of C mod m and the difference 
# between two consecutive parts is NOT in DI

restrictedPnk := proc(n, k, m, C, DI) local res, newLargest, prev, p:
    if n = 0 and k=0 then
        return {[]}:
    fi:

    if k > n or k < 0 then
        return {}:
    fi:

    if not((k mod m) in C) then
        return {}:
    fi:

    if k = n then
        return {[k]}:
    fi:

    res := {}:
    for newLargest from 1 to k do
        if k - newLargest in DI then
            next:
        fi:

        prev := restrictedPnk(n - k, newLargest, m, C, DI):
        res := {op(res), seq([k, op(p)], p in prev)}:
    od:

    res:
end:



#restrictedPn(n,m,C,DI): The set of all partitions of n
# s.t. all parts are congruent to a member of C mod m and the difference 
# between two consecutive parts is NOT in DI
restrictedPn := proc(n, m, C, DI) local k:
    {seq(op(restrictedPnk(n, k, m, C, DI)), k=1..n)}:
end:



#NuRestrictedPnk(n, k, m, C, DI): The number of partitions of n with largest part k
# s.t. all parts are congruent to a member of C mod m
# and the difference between two consecutive parts is NOT in DI

NuRestrictedPnk := proc(n, k, m, C, DI) local newLargest, res:
    if n = 0 and k=0 then
        return 1:
    fi:

    if k > n or k < 0 then
        return 0:
    fi:

    if not((k mod m) in C) then
        return 0:
    fi:

    if k = n then
        return 1:
    fi:

    res := 0:
    for newLargest from 1 to k do
        if k - newLargest in DI then
            next:
        fi:

        res := res + NuRestrictedPnk(n - k, newLargest, m, C, DI):
    od:

    res:
end:

#NuRestrictedPn(n, m, C, DI): The number of partitions of n 
# s.t. all parts are congruent to a member of C mod m
# and the difference between two consecutive parts is NOT in DI

NuRestrictedPn := proc(n, m, C, DI) local k:
    add(NuRestrictedPnk(n, k, m, C, DI), k=1..n):
end:


############################	Random Variable Procedures	############################
#numberParts(pi): Inputs an integer partition (as a decreasing list of integers)
# Outputs the number of parts
numberParts:=proc(pi):
	nops(pi):
end:

#largestPart(pi): Inputs an integer partition (as a decreasing list of integers)
# Outputs the largest part

largestPart:=proc(pi): 
	if pi = [] then
		RETURN(FAIL):
	fi:
	
	pi[1]:
end:

#smallestPart(pi): Inputs an integer partition (as a decreasing list of integers)
# Outputs the smallest part
smallestPart:=proc(pi):
	if pi = [] then
		RETURN(FAIL):
	fi:
	
	pi[-1]:
end:

#distinctParts(pi): Inputs an integer partition (as a decreasing list of integers)
# Outputs the number of distinct parts
	
distinctParts:= proc(pi)
    nops(convert(pi, set)):
end:



############################	Estimation Procedures	############################

#LD(L): Inputs a list of positive integers L (of n:=nops(L) members) 
# outputs an integer i from 1 to n with the prob. of i being proportional to L[i]

LD:=proc(L) local n,i,su,r:
	n:=nops(L):
	r:=rand(1..convert(L,`+`))():
	su:=0:

	for i from 1 to n do
		su:=su+L[i]:
		if r<=su then
			RETURN(i):
		fi:
	od:
end:


#RandRestrictedPnk(n, k, m, C, DI): Generates, uniformly at random, a partition of n with largest part k
# s.t. all parts are congruent to a member of C mod m
# and the difference between two consecutive parts is NOT in DI

RandRestrictedPnk := proc(n, k, m, C, DI) local L,L1,Die1,k1,S1,i:
    if not (type(n,integer) and type(k,integer) and n>=1 and k>=1) then
		RETURN([]):
	fi:
	
	if k>n then 
		RETURN([]):
	fi:
	
	S1:=C mod m: 
	
	if not ( member(k mod m,S1) ) then
		RETURN([]):
	fi:
	
	if k=n then
		if member(n mod m,S1) then 
			RETURN([n]):
		fi:
	fi:
	
	L:=[k]:
	
	L1:=convert({seq(1..min(n-k,k))} minus {seq(k-DI[i],i=1..nops(DI))},list):

	Die1:=[seq(NuRestrictedPnk(n-k, i, m, C, DI),i in L1)]:
	
	k1:=L1[LD(Die1)]:
	
	[k,op(RandRestrictedPnk(n-k, k1, m, C, DI))]:
end:


#RandRestrictedPn(n m, C, DI): Generates, uniformly at random, a partition of n 
# s.t. all parts are congruent to a member of C mod m
# and the difference between two consecutive parts is NOT in DI

RandRestrictedPn := proc(n, m, C, DI) local Die1,k,i:
	if not (type(n,integer) and n>=1) then
		RETURN([]):
	fi:

	if n=0 then RETURN([]): fi:
	
	Die1:=[seq(NuRestrictedPnk(n, i, m, C, DI),i=1..n)]:
	k:=LD(Die1): 
	
	RandRestrictedPnk(n, k, m, C, DI)
end:

#Ave(X): inputs a RV X given as  a list of numbers)
# outputs the average (alias expectation, alias mean)

Ave:=proc(X) local i:
	add(X[i],i=1..nops(X))/nops(X):
end:


#Mom(X,k): inputs a RV X (given as  a list of numbers)
# outputs the k-th moment about the mean

Mom:=proc(X,k) local i, mu:
	mu:=Ave(X):
	add((X[i]-mu)^k,i=1..nops(X))/nops(X):
end:


#SA(X,K): inputs a discrete random variable X (given as  a list of numbers) and a 
#positive integer K 
#outputs a list of length K whose 
#(i) first entry is the average
#(ii) second entry variance
#The remaining entries are the scaled third moment through K-th moment 

SA:=proc(X,K) local i,M,mu,sig2:
	if not (type(X,list) and type(K,integer) and K>=1) then
		print(`Bad input`):
		RETURN(FAIL):
	fi:
	
	mu:=Ave(X):
	
	if K=1 then
		RETURN([mu]):
	fi:
	
	sig2:=Mom(X,2):
	M:=[mu,sig2]:
	
	if sig2=0 then
		RETURN(M):
	fi:
	
	[op(M), seq(evalf(Mom(X,i)/sig2^(i/2)),i=3..K)]:
end:


#EstStatAnalk(n,k,K,N,rv,m,C,DI): estimates the expectation, variance, and 
# the 3rd- through K-th scaled moments of the random variable rv, of partitions 
# of n with largest part k s.t. 
# all parts are congruent to a member of C mod m
# and the difference between two consecutive parts is NOT in DI
# by generating, uniformly at random, a partition with the desired conditions
# N times (N should be large)

EstStatAnalk:=proc(n,k,K,N,rv,m,C,DI) local i,pi,X:
	X:=[]:
	for i from 1 to K do
		pi:=RandRestrictedPnk(n,k, m, C, DI):
		X:=[op(X),rv(pi)]:
	od:
	
	SA(X,N):
end:


#EstStatAnal(n,K,N,rv,m,C,DI): estimates the expectation, variance, and 
# the 3rd- through K-th scaled moments of the random variable rv, of partitions of n s.t. 
# all parts are congruent to a member of C mod m
# and the difference between two consecutive parts is NOT in DI
# by generating, uniformly at random, a partition with the desired conditions
# N times (N should be large)

EstStatAnal:=proc(n,K,N,rv,m,C,DI) local i,pi,X:
	X:=[]:
	for i from 1 to K do
		pi:=RandRestrictedPn(n, m, C, DI):
		X:=[op(X),rv(pi)]:
	od:
	
	SA(X,N):
end:



############################	Generating Function Procedures	############################

#partitionGFk(f, n, t, k, m, C, DI): The generating function of the random variable f defined over
#integer partitions of n with largest part k
#s.t. all parts are congruent to a member of C mod m
#and the difference between two consecutive parts is NOT in DI
#i.e. sum(t^f(p),p in restrictedPnk(n,k,m,C,DI))

partitionGFk := proc(f, n, t, k, m, C, DI) local p:
       add(t^(f(p)), p in restrictedPnk(n,k, m, C, DI)):
end:



#partitionGF(f, n, t, m, C, DI): The generating function of the random variable f defined over
#integer partitions of n s.t. all parts are congruent to a member of C mod m
#and the difference between two consecutive parts is NOT in DI
#i.e. sum(t^f(p),p in restrictedPn(n,m,C,DI))
#Try:
# partitionGF(p -> 0, 10, t, 1, {0}, {})
# partitionGF(largestPart, 10, t, 1, {0}, {})
# partitionGF(smallestPart, 10, t, 1, {0}, {})
# partitionGF(numberParts, 10, t, 1, {0}, {})

partitionGF := proc(f, n, t, m, C, DI) local p:
    add(t^(f(p)), p in restrictedPn(n, m, C, DI)):
end:

#partitionPGFk(f, n, t, k, m, C, DI): The probability generating function of the random variable f defined over
#integer partitions of n with largest part k
#s.t. all parts are congruent to a member of C mod m
#and the difference between two consecutive parts is NOT in DI

partitionPGFk := proc(f, n, t, k, m, C, DI)
    local N, total, p:
    N := 0:
    total := 0:
    for p in restrictedPnk(n, k, m, C, DI) do
        total := total + t^(f(p)):
        N := N + 1:
    od:

    total / N:
end:

#partitionPGF(f, n, t, m, C, DI): The probability generating function of the random variable f defined over
#integer partitions of n 
#s.t. all parts are congruent to a member of C mod m
#and the difference between two consecutive parts is NOT in DI

partitionPGF := proc(f, n, t, m, C, DI)
    local N, total, p:
    N := 0:
    total := 0:
    for p in restrictedPn(n, m, C, DI) do
        total := total + t^(f(p)):
        N := N + 1:
    od:

    total / N:
end:



#Zigf(k,i,t):

Zigf := proc(k, i, t) local j:
    mul(1 / (1 - t^j), j=1..i-1) *
    mul(1 / (1 - t^j), j=i+1..k) *
    t^i / (1 - t^i)^2:
end:



#Zgf(k,t): The generating function in variable t of the sequence 
# a_n = total number of parts in all integer partitions of n with largest part at most k
# i.e.  sum( Z(n,k)t^n, t=0..infinity), where
#  Z(n, k) = sum( sum(number of parts of p, p in P(n, k1)), k1=1..k)
Zgf := proc(k, t) local i:
    add(Zigf(k, i, t), i=1..k):
end:


#Sgf(k,t): The generating function in variable t of the sequence 
# a_n = total number of parts in all integer partitions of n with largest part k
# i.e.  sum( S(n,k)t^n, t=0..infinity), where
#  S(n, k) = sum(number of parts of p, p in P(n, k))

Sgf := proc(k, t)
    Zgf(k, t) - Zgf(k - 1, t):
end:

 
#NiceSeq(k,N): Uses the first N terms of the Taylor series expansion of the 
#generating function, and tries to find a "nice" sequence for the total number
#of parts in all integer partitions of n with largest part k by separating cases for 
#n mod lcm(1,...,k).
#It outputs a list where the i-th entry is the total number of parts for partitions of 
#lcm(1,...,k)*m+i 

NiceSeq:=proc(k,N) local f,l,m1,PS,i,s,j,guess:
	f:= taylor(Sgf(k,t),t=0,N):
	l:=lcm(seq(1..k)): 
	m1:=floor(N/l):
	PS:=[]:

	for i from 0 to l-1 do
		s:=[seq(coeff(f,t,l*j+i), j=0..m1-1)]: 
		guess:=convert(guessgf(s,t)[1],FormalPowerSeries,t,m):
		if guess=FAIL[1] then RETURN(FAIL): fi:
		PS:=[op(PS),factor(op(op(guess)[1])[1])]: 
	od:
	PS:
end:

# powerDiff(f, q, m): (qD_q)^m f, where D_q is the differentiation operator
# with respect to q.
powerDiff := proc(f, q, m)
    local fp, k:

    fp := f:
    for k from 1 to m do
        fp := q * diff(fp, q):
    od:

    fp:
end:

# Qgf(k, m, t): The generating function in t with respect to n of the sequence
# Q(n, k, m) = sum(L(p)^m, p), where the sum extends over all partitions of n
# with largest part AT MOST k.
# Try:
#   Qgf(1, 1, t)
#   Qgf(1, 2, t)
#   Qgf(2, 2, t)
Qgf := proc(k, m, t)
    local q, f,j:
    f := mul(1 / (1 - q * t^j), j=1..k):
    f := powerDiff(f, q, m):
    subs(q = 1, f):
end:

# Q(n, k, m, t): sum(L(p)^m, p), where the sum extends over all partitions of n
# with largest part AT MOST k.
# Try:
#   Q(10, 2, 1, t)
Q := proc(n, k, m)
    local f, t:
    f := Qgf(k, m, t):
    coeff(taylor(f, t, n + 1), t, n):
end:

# QExp(n, k, m): The mth moment about 0 of the random variable "number of parts
# of a random partition of n with largest part at most k."
# Try:
#   QExp(n, k, 1) should be roughly n * harmonic(k) / k.
#   evalf(seq(QExp(n, 2, 1) / n, n=1..20));
#   evalf(seq(QExp(n, 2, 2) / n^2, n=1..20));
QExp := proc(n, k, m)
    local pnk,rnk,j:
    rnk := add(NuRestrictedPnk(n, j, 1, {0}, {}), j=0..k):
    Q(n, k, m) / rnk:
end:

# QMom(n, k, m): The mth moment about the mean of the random variable "number
# of parts of a random partition of n with largest part at most k."
# Try:
#   evalf(seq(QMom(n, 2, 1), n=1..20));
#   evalf(seq(QMom(n, 2, 2)/ n^2, n=1..20));
#   evalf(seq(QMom(n, 2, 3), n=1..20));
QMom := proc(n, k, m)
    local mu, j:
    mu := QExp(n, k, 1):
    add(binomial(m, j) * QExp(n, k, j) * (-mu)^(m - j), j=0..m):
end:


#QClosed(n,k,m): The closed form of Q(n,k,m)

QClosed := proc(n, k, m)
    local f, t, formula:
    f := factor(Qgf(k, m, t)): 
    formula := op(1, convert(f, FormalPowerSeries,t,j)) / t^j:

    subs(j = n, factor(formula)):
end:


#Child_PnkRest(k,m,C,DI,x): inputs positive integers k and m, sets C and DI, and variable x.
#Outputs (1) the equation relating Z[k] to its "children", and (2) its "children",
# where Z[k] is the set of partitions containing the empty partition and partitions where the largest part 
#is k, all parts mod m are congruent to an element of C, and the difference between any two terms is not in
#the set DI

Child_PnkRest:=proc(k,m,C,DI,x) local L,i,r:
	if not member(k mod m,C) then 
		RETURN({1}): 
	fi:
	L:={seq(seq(i*m+r,i=0..floor((k-r)/m)), r in C)}:  
	L:=L minus {seq(k-i,i in DI)} minus {0}: 

	{-Z[k]+expand(x^k*add(Z[i],i in L))+x^k}, [seq(Z[i],i in L)]:
	end:

#f_PnkRest(k,m,C,DI,x): inputs positive integers k and m, sets C and DI, and variable x.
#Outputs the generating function of the set of integer partitions where the largest part  is k, 
#all parts mod m are congruent to an element of C, and the difference between any two terms is not in
#the set DI


f_PnkRest:=proc(k,m,C,DI,x) local r,ALREADYDONE,STILLDO,STATE,gu,eq,i,var,var1,lu,gu1,i1,L:
	if not member(k mod m,C) then 
		RETURN(0): 
	fi:   
	ALREADYDONE:={}:
	L:={seq(seq(i*m+r,i=0..floor((k-r)/m)), r in C)} minus {0}:  
	STILLDO:=convert(L,list):
	eq:={}: 
	
	while STILLDO<>[] do
		STATE:=STILLDO[1],m,C,DI: 
		ALREADYDONE=ALREADYDONE union {STATE}:
		gu:=Child_PnkRest(STATE,x):  
		eq:=eq union gu[1]:
		STILLDO:=STILLDO[2..-1]: 
	od: 

	var:={seq(Z[i],i in L)}:
	var1:=var minus {Z[k]}:  
	var:=[op(var1), Z[k]]: 
	gu:=subs(Z[k]=P,Groebner[Basis](eq,plex(op(var)))[1]):
	simplify(solve(gu,P)):
end: