knitr::opts_chunk$set(echo = TRUE)
![](`r system.file("help/figures/partitions.png", package = "partitions")`){width=10%}
A generating function for restricted partitions---originally due, as
far as I can tell, to Wilf (2000)---is presented and R idiom using the
spray
package given. The generating function approach is shown to
be not particularly efficient compared to the direct enumeration used
in
the partitions
package.
The partitions
package gives functionality for various integer
partition enumeration problems including that of restricted
partitions, function restrictedparts()
:
library("partitions") library("disordR") (jj <- restrictedparts(7,3)) ncol(jj)
Here I will consider function R()
, which calculates the size of
the matrix required:
R(3,7,include.zero=TRUE)
Function R()
is very basic; all it does is to go through all the
restricted partitions, counting them one by one until the recursion
bottoms out:
unsigned int numbrestrictedparts(int *x, const int m){ unsigned int count=1; while(c_nextrestrictedpart(x, &m)==0){ count++; } return count;
To implement a potentially more efficient method, we can use generating functions. Here we follow Wilf and, using his terminology, define an infinite polynomial $P(x,y)$ as follows:
\begin{equation}\label{GF} P(x,y) = \prod_{r=0}^\infty\frac{1}{1-x^ry} \end{equation}
Or, expanding:
\begin{equation}\label{GF2} P(x,y) = \left(1+y+y^2+y^3 + \cdots\right) \left(1 + xy + x^2y^2 + x^3y^3+\cdots\right) %\left(1 + x^2y + x^4y^2+x^6y^3+\cdots\right) \cdots \left(1+x^ry+x^{2r}y^{2}+x^{3r}y^{3}+\cdots\right)\cdots \end{equation}
The power of $x$ counts the total of the chosen integers (the size of
the partition), and the power of $y$ counts the number of integers
chosen (the length of the partition). Thus the number of partitions
of $k$ into at most $n$ parts is the coefficient of $x^ky^n$ in
$P(x,y)$. In numerical work it is convenient and efficient to ignore
terms with a power of $x$ higher than $n$ (sum of integers chosen
exceeds $n$), or with power of $y$ higher than $k$ (number of integers
chosen exceeds $k$). Taking R(3,7,include.zero=TRUE)
as an example
we would truncate equation \ref{GF2} as follows:
\begin{eqnarray}\label{GF3} P(x,y) = &{}& \left(1+y+y^2+y^3\right) \left(1 + xy + x^2y^2 + x^3y^3\right) \left(1 + x^2y + x^4y^2+x^6y^3\right)\times\nonumber\ &{}&\left(1 + x^3y + x^6y^2\right) \left(1 + x^4y\right) \left(1 + x^5y\right) \left(1 + x^6y\right) \left(1 + x^7y\right) \end{eqnarray}
and the coefficients of $P(x,y)$ up to $x^7y^3$ would correctly count the restricted partitions. Note that we need consider only at most four terms in each bracket (powers of $y$ above three being irrelevant) and we may stop the continued product at the $x^7$ term as further brackets contain only one and powers of $x$ above the eighth.
The R implementation uses the spray
package, in particular
function ooom(x)
which returns $\frac{1}{1-x}$.
library("spray") R_gf <- function(k,n){ # version 1 x <- spray(cbind(1,0)) y <- spray(cbind(0,1)) P <- ooom(y,k) # term x^0; number of zeros chosen for(i in seq_len(k)){ # starts at 1 P <- P*ooom(x^i*y,n) } return(drop(coeffs(P[k,n]))) }
Thus
R_gf(7,3)
We can do slightly better in terms of efficiency by ruthlessly cutting out powers higher than needed:
strip <- function(P,k,n){ # strips out powers higher than needed ind <- index(P) val <- elements(coeffs(P)) wanted <- (ind[,1] <= k) & (ind[,2] <= n) spray(ind[wanted,],val[wanted]) }
which is used here:
R_gf2 <- function(k,n,give_poly=FALSE){ x <- spray(cbind(x=1,y=0)) y <- spray(cbind(x=0,y=1)) P <- ooom(y,k) # term x^0 for(i in seq_len(k)){ # starts at 1 P <- strip(P*ooom(spray(cbind(i,0))*y, min(n,ceiling(k/i))),k,n) } if(give_poly){ return(P) } else { return(drop(coeffs(P[k,n]))) } }
then
R_gf2(7,3)
We can test the computational efficiency of the generating function approach using larger values of $k$ and $n$:
k <- 140 n <- 4 system.time(jj1 <- R(n,k,include.zero=TRUE)) system.time(jj2 <- R_gf2(k,n)) jj1==jj2
So the generating function approach is not particularly efficient, at
least not in this sort of use-case with the spray
package. It
might be better with the mvp
package; I don't know.
Of course, R_gf2()
calculates the generating polynomial
which gives very much more information than is returned. Perhaps this
is why it is so slow compared to function R()
, although it is
surprising to see direct enumeration so heavily outperforming a
generating function.
H. S. Wilf, 2000. "Lectures on Integer Partitions"
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.