Distinguishable balls in distinguishable urns of a maximal size Let's assume we have $u\geqslant 1$ distinguishable urns, each of which can hold at most $m\geqslant 1$ balls, and $b\geqslant 1$ distinguishable balls.  In how many ways can we
(a) distribute all $b$ balls into the urns?
(b) distribute any number of balls up to $b$ into the urns?
The final configuration does matter. (E.g., for $u=2$, $m=3$, $b=4$,  urns U1,U2, and balls B1–B4, the configurations U1={B1,B2},U2={B3,B4} and U1={B3,B4},U2={B1,B2} and U1={B1,B2,B3},U2={B4} are three different configurations.)
The time order of putting the balls into urns doesn't matter. (For example, putting balls B1 and B2 into urn U1 and balls B3 and B4 into urn U2 gives the same result as putting B2 into U1, then B3 into U2, then B1 into U1, and then B4 into U2; both distribution oders lead to the same set {B1,B2} being the contents of U1 and {B3,B4} the contents of U2.)
In a boringly precise formula language, we count the number of tuples of the form $(U_1,\dots,U_u)$ where $\bigcup_{i=1}^u U_i=\{1,\dotsc,b\}$ (for (a); $\bigcup_{i=1}^u U_i \subseteq \{1,\dotsc,b\}$ for (b)), and $\lvert U_i\rvert \leqslant m$ ($1\leqslant i\leqslant u$), and the family of sets $\{U_i\mid 1\leqslant i\leqslant u\}$ is pairwise disjoint.
I need an exact closed-form expression, nonrecursive and involving only elementary operations, exponentiation, logarithm, factorials, big sums, big products, and multinomial coefficients. If unavailable, a polynomial-time algorithm (in the length of the input consisting of the binary representations of $b$, $m$, and $u$) would be great. If even that unavailable, a good lower bound (or an asymptotic formula with exact error margins) would be great. If even that is unavailable, a recursive formula would be great.
Remark. In my application, $u \leqslant b \leqslant mu$ and $m \leqslant b$; these assumptions might hypothetically help.
 A: So, if I understood what you mean to say

*

*you have $b$ dist. balls, i.e. labeled $1,2, \ldots , b$. You shuffle (permute) them , i.e. you relabel them,  -> no change wrt if they are labelled  orderly;

*you take the balls one by one, in sequence, and put it into any one of the free spaces available in the urn;

*since you do not take into account the order of the balls within each urn, that is the same as starting with the original labelling $1,2, \ldots ,b$,
deposit the balls sequentially and end up with a u-ple (or list) of
sub-sets(i.e. ordered) taken from the set $\{1,2, \ldots , b \}$ , and having $0,1, \ldots, max=m$ elements, that is to say something like
$$
\left[ {\left\{ {1,4,6} \right\},\left\{ \emptyset  \right\},\left\{ {2,5} \right\},\left\{ \emptyset  \right\},\left\{ 3 \right\}} \right]
$$
Under this acception, and in particular regarding the ball laying process, and with a given number $b$ of balls (case a) ) all the u-ples are equi-probable.
That premised, allow me to change your symbols
$$
{\rm N}{\rm .}\,{\rm balls}:\;b\; \Rightarrow \;s\;,\quad {\rm balls/box}\,
{\rm max}:\;m\; \Rightarrow \;r\;,\quad {\rm N}{\rm .}\,{\rm boxes}:\;\;u\; \Rightarrow \;m
$$
so that I can introduce a sketch I already prepared time ago for myself.

I hope that the sketch be quite self explanatory: we start from the Sets of non-empty sub-sets (counted by $ L_{\,b\,U\,D\,N} $) and from that we compute
the Lists of non-empty subsets which clearly are $m!$ times the precedent, and the cumulative, i.e.  possibly empty sub-sets case follow immediately.
Thus the key here are the Restrained Stirling N. 2nd kind.
You can refer for instance to this paper.
I will just remark here that they are related to the Bell Polynomials
and summarize a couple of the equivalent ways to define them
$$
\eqalign{
  & L_{\,b\,U\,D\,N} (s,r,m) = \left\{ \matrix{  s \cr  m \cr}  \right\}_{\,r}  =   \cr 
  &  = \sum\limits_{\left\{ {\matrix{
   {1\, \le \,p_{\,1}  \le p_{\,2}  \le  \cdots  \le p_{\,m} \, \le \,r}  \cr 
   {p_{\,1}  + p_{\,2}  + \, \cdots  + p_{\,m}  = s}  \cr 
 } } \right.} {{1 \over {t_{\,1} !t_{\,2} !\, \cdots \,t_{\,m} !}}
\left( \matrix{ s \cr   p_{\,1} ,\,p_{\,2} ,\, \cdots ,p_{\,m}  \cr}  \right)}  =   \cr 
  &  = \sum\limits_{\left\{ {\matrix{
   {0\, \le \,t_{\,j} \,\left( {\, \le \,m} \right)}  \cr 
   {t_{\,1}  + t_{\,2}  + \, \cdots  + t_{\,r}  = m}  \cr 
   {1\,t_{\,1}  + 2\,t_{\,2}  + \, \cdots  + r\,t_{\,r}  = s}  \cr 
 } } \right.} {\,{{s!} \over {\prod\limits_{1\, \le \,\,k\, \le \,r}
 {t_{\,k} !\left( {k!} \right)^{\,t_{\,k} } } }}}  \cr} 
$$
In the second line, the sum is taken over all the partitions $(p_1,p_2, \ldots , p_m)$ of $s$ into $m$ parts
of max size $r$.
$t_k$ counts the number of times that the part $p_k$ appears : $0 \le t_k \le m, \; (t_1+t_2+ \cdots + t_m)=m$ .
$$
\eqalign{
  & \left[ {p_{\,1} ,\,p_{\,2} ,\, \cdots ,p_{\,m} } \right]\;:\;
{\rm partition}\,{\rm of}\,s\,{\rm into}\,m\,{\rm parts}\, \le \,r  \cr 
  & \underbrace {p_{\,1} }_{t_{\,1} } < \,\underbrace {p_{\,2}  =
 p_{\,3}  =  \cdots  = p_{\,q} }_{t_{\,2} }\, <
 \underbrace {p_{\,q + 1} }_{t_{\,q + 1} } <  \cdots  <
 \underbrace {p_{\,m - u}  =  \cdots  = p_{\,m} }_{t_{\,m - u} } \cr} 
$$
In the third line, we rank the possible parts as $1 \le k \le r$ and assign to each the number of times $0 \le t_k \le m$
that the part $k$ appears in the partition.
Then, allowing some parts to be null, we get
$$
N_{\,b\,U\,D\,N} (s,r,m) = \sum\limits_{\left\{ {\matrix{
   {0\, \le \,p_{\,1}  \le p_{\,2}  \le  \cdots  \le p_{\,m} \, \le \,r}  \cr 
   {p_{\,1}  + p_{\,2}  + \, \cdots  + p_{\,m}  = s}  \cr 
 } } \right.} {{1 \over {t_{\,1} !t_{\,2} !\, \cdots \,t_{\,m} !}}
\left( \matrix{  s \cr   p_{\,1} ,\,p_{\,2} ,\, \cdots ,p_{\,m}  \cr}  \right)} 
$$
While, passing from partitions to Compositions
respectively "strong" and "weak", we get
$$
L_{\,b\,D\,D\,N} (s,r,m) = \sum\limits_{\left\{ {\matrix{
   {1\, \le \,c_{\,j} \, \le \,r}  \cr  {c_{\,1}  + c_{\,2}  + \, \cdots  + c_{\,m}  = s}  \cr  } } \right.}
 {\left( \matrix{ s \cr  c_{\,1} ,\,c_{\,2} ,\, \cdots ,c_{\,m}  \cr}  \right)} 
$$
$$
N_{\,b\,D\,D\,N} (s,r,m) = \sum\limits_{\left\{ {\matrix{
   {0\, \le \,c_{\,j} \, \le \,r}  \cr {c_{\,1}  + c_{\,2}  + \, \cdots  + c_{\,m}  = s}  \cr } } \right.}
 {\left( \matrix{  s \cr  c_{\,1} ,\,c_{\,2} ,\, \cdots ,c_{\,m}  \cr}  \right)} 
$$
You may verify that the above formulations in terms of multinomials satisfy the relations shown in the sketch.
Also, taking for example $r = 3$,  the above formulas check with
$ L_{\,b\,U\,D\,N} (s,3,m)$  : OEIS A111246
$ L_{\,b\,D\,D\,N} (s,3,m)$ : OEIS A189804
$ N_{\,b\,D\,D\,N} (s,3,m)$ : OEIS A234574
A: As far as I know, there is no non-recursive closed form expression for this. However, it can be effectively computed as follows.
If $N(b,m,u)$ is the number of placements of exactly $b$ balls, then
$$
N(b,m,u)=\sum_{i=0}^{\min(b,m)} \binom{b}{i}N(b-i,m,u-1)
$$
which follows by conditioning on all the ways to place at most $m$ balls in the first urn.
This lets you compute $N(b,m,u)$ by filling out the table indexed by $\{0,1,\dots,b\}\times \{1,\dots,u\}$ whose entry in the $r^{th}$ row and $c^{th}$ columns is $N(r,m,c)$. Each entry takes $O(m)$ arithmetic operations to fill out, so the complexity is $O(b\times m\times u)=O((mu)^2)$.
For a faster algorithm, you can use exponential generating functions to write
$$
N(b,m,u)=b![x^b]\left(1+\frac{x^1}{1!}+\frac{x^2}{2!}+\dots+\frac{x^{m}}{m!}\right)^u
$$
The notation $[x^n]f(x)$ refers to the coefficient of $x^n$ in the polynomial (or power series) $f(x)$.
This shows that $N(b,m,u)$ can be computed using $\log u$ polynomial multiplications between polynomials of degree $b$ (use exponentiation by squaring, and only keep track of terms up to degree $b$). If you use FFT multiplication, the complexity is $$O(b\log b\times \log u)=O(m\times u\times \log(mu)\times \log u).$$ Again, this complexity is measured in number of arithmetic operations.
Finally, both of these methods actually compute the entire list $\{N(a,m,u):0\le a\le b\}$ at no extra cost, so you can compute the number of ways to place at most $b$ balls by adding up all the values of $N(a,m,u)$ for $a\in \{0,1,\dots,b\}$.
