Constrained sum of product binomial coefficients

The classic example of the sum of product binomial coefficients is Vandermonde's Identity:

$$\sum_{a_{1,1}+a_{1,2}+a_{1,3}=x}\binom{n}{a_{1,1}}\binom{n}{a_{1,2}}\binom{n}{a_{1,3}}=\binom{3n}{x}$$

However, I'm interested in the following twist. To motivate the question, consider the simple case in which I'd like to compute the following sum. Assume that $0\leq c_1\leq n$, $0\leq c_2\leq 2c_1$, $n$ and $x\leq 3c_1$ are given.

$$\sum_{\bigcap\left(\substack{a_{1,1}\leq c_1,a_{1,2}\leq c_1,a_{1,3}\leq c_1\\ a_{1,1}+a_{1,2}=a_{2,1}\leq c_2,a_{1,2}+a_{1,3}=a_{2,2}\leq c_2\\a_{1,1}+a_{1,2}+a_{1,3}=x}\right)}\binom{n}{a_{1,1}}\binom{n}{a_{1,2}}\binom{n}{a_{1,3}}$$ After some algebra, this can be simplified to:

$$\sum_{\bigcap\left(\substack{x-c_2\leq a_{1,1}\leq c_1\\x-c_1\leq a_{2,1}\leq c_2\\a_{2,1}-a_{1,1}\leq c_1}\right)}\binom{n}{a_{1,1}}\binom{n}{a_{2,1}-a_{1,1}}\binom{n}{x-a_{2,1}}$$

which only involves two terms $a_{1,1}$ and $a_{2,1}$ that we can sum over.

This is fairly simple to code and corresponds to the case $L=3$. However, I will eventually run into the complicated case where I have the following constrained sum for general $L$:

$$\sum_{\bigcap\left(\substack{\bigcup_{k=1}^{L-1}\bigcup_{r=1}^{L-k+1}\{a_{k,r}\leq c_k \}\\ a_{L,1}=x}\right)}\binom{n}{a_{1,1}}\binom{n}{a_{1,2}}\binom{n}{a_{1,3}}\cdots\binom{n}{a_{1,L}}$$

or alternatively, like above,

$$\sum_{\bigcap\left(\substack{ \bigcup_{k=1}^{L-1}\{x-c_{L-k}\leq a_{k,1}\leq c_k \} \\ \bigcup_{k=1}^{L-2}\bigcup_{r=1}^{L-k-1}\{a_{k+r,1}-a_{r,1}\leq c_k \}}\right)}\binom{n}{a_{1,1}}\binom{n}{a_{2,1}-a_{1,1}}\binom{n}{a_{3,1}-a_{2,1}}\cdots\binom{n}{a_{L-1,1}-a_{L-2,1}}\binom{n}{x-a_{L-1,1}}$$

where $a_{k,r}=\sum_{s=r}^{r+k-1} a_{1,s}$

I'm interested in exact solutions to the problem. Is there an efficient way to calculate this sum for general $L$ using R (or similar languages as MATLAB and Python)?

I believe I have a working solution that works for general $L>1$:

binom.coeff <- function(x,l,size,cvec){
sum(apply(valid.counts(l=l,cvec=cvec,x=x),1,function(j) prod(choose(size,j))))
}
valid.counts <- function(l,cvec,x){
vec = seq(from=0,to=cvec[1])
lst = lapply(numeric(l), function(i) vec)
mat = as.matrix(expand.grid(lst))

if(l==2){
return(mat[which(rowSums(mat)==x),])
}
else{
#k=1 is redundant since all must be less than or equal to cvec[1]
for (k in seq(from=2, to=l-1)){
partial = NULL
for (r in seq(l-k+1)){
partial = which(rowSums(mat[,r:(r+k-1),drop=FALSE])<=cvec[k])
mat = mat[unique(partial),]
}
}
return(mat[which(rowSums(mat)==x),])
}
}