Counting in how many ways rocks can be put in boxes How can I figure out the following questions?
How many possible combinations can be done by having 26 boxes and 15 red rocks, and 15 black rocks? Each box can have up to 15 rocks in it. We can have none or all 15 rocks in one box, or distributed in the boxes. All the rocks do not have to be used. But we can not mix the red and black rocks together. Meaning if one of the boxes has at least 1 black rock in it, we can not add any red rocks to that box. The rocks have no value, other than that they are red or black.
 A: Suppose:


*

*$i$ boxes contain red rocks.

*$j$ boxes contain black rocks.

*$s$ red rocks are in boxes.

*$t$ black rocks are in boxes.


(No, I'm not trying to write a Dr Seuss book.)
With these restrictions, there are ${26 \choose {i,j,26-i-j}}$ ways to choose the boxes for the rocks (note: this is the multinomial coefficient).  Once we've decided on which boxes red rocks belong to, we need to partition them into boxes.  Let $N_{n,k}$ be the number of ordered partitions of $n$ into $k$ non-zero parts.  Then the number of ways of putting the red rocks into their chosen boxes is $N_{s,i}$.  Similarly, the number of ways of putting the black rocks into their chosen boxes is $N_{t,j}$.
Hence, if we sum this over all valid $i,j,s,t$ we get the result \[\sum_{i=0}^{15} \sum_{j=0}^{15} \sum_{s=0}^{15} \sum_{t=0}^{15} {26 \choose {i,j,26-i-j}} N_{s,i} N_{t,j}.\]
(Technically, we need the upper bound on $j$ to be $\mathrm{min}(15,26-j)$, but if we assume that all multinomial coefficients with a negative parameter are $0$, then it's fine.)
This can be computed in GAP by:
Multinomial:=function(n,L)
  local x,i;
  if(ForAny(L,i->i<0)) then return 0; fi;
  if(Sum(L)<>n) then return 0; fi;
  if(n=0) then return 1; fi;
  x:=Factorial(n);
  for i in L do x:=x/Factorial(i); od;
  return x;
end;;

num:=0;
for i in [0..15] do
  for j in [0..15] do
    for s in [0..15] do
      for t in [0..15] do
        num:=num+Multinomial(26,[i,j,26-i-j])*NrOrderedPartitions(s,i)*NrOrderedPartitions(t,j);
      od;
    od;
  od;
od;
Print(num,"\n");

which gave the result 38996875553354247211.
As a quick check of this number, we can act on the set of arrangements (by permuting box indices) with the group $\langle (1,2,\ldots,p)\rangle$ for prime $p$ in the range $15<p \leq 26$.  Orbits have size $p$ unless the first $p$ boxes are empty (since $p>15$).
The number of arrangements with the first $p$ boxes empty is listed below (computed using an edited version of the above code), and the congruences match up.
p    nr arrangements with boxes 1..p empty
17   25995200573833414041
19   3229548304818580291
23   1017952468781

As another check, we can act on the set of arrangements with the group $\langle (1,2,\ldots,13)(14,15,\ldots,26)\rangle$.  Non-identity permutations in this group have no fixed points and have order $13$.  Thus, orbits have size $13$, unless they are one of the ones listed below:
............. .............
............. rrrrrrrrrrrrr
............. bbbbbbbbbbbbb
rrrrrrrrrrrrr .............
rrrrrrrrrrrrr bbbbbbbbbbbbb
bbbbbbbbbbbbb .............
bbbbbbbbbbbbb rrrrrrrrrrrrr

(note: we must have only one copy of a symbol in each box, since $2 \times 13>15$).  This is consistent with $38996875553354247211 \equiv 7 \pmod {13}$.
