Sum of $k$ elements of $\{1,...,p\}$ divisible by $p$ Let $k,p \in \mathbb{N}$, $A=\{1,2,\ldots p\}$. In how many ways can we pick $k$ (not necessarily distinct) numbers from $A$ such that their sum is $0\;{\rm mod}\; p$?
 A: Partial Answer (The Case $\gcd(k,p)=1$)
If $\gcd(k,p)=1$, then the answer is $n_{p,k}:=\frac{1}{p}\,N_{p,k}$, where $N_{p,k}$ is the number of multisets of size $k$ consisting of elements of $\{1,2,\ldots,p\}$.  The number $N_{p,k}$ is simply $\binom{p+k-1}{k}$, so
$$n_{p,k}=\frac{1}{p}\,\binom{p+k-1}{k}\,.$$ It can be easily seen that the operation $s$ on a multiset of size $k$ consisting of elements of $A$ increases the sum (in modulo $p$) by the amount $k$, where $s$ increases each element of $A$ by $1$ (in modulo $p$).  
For example, if $k=4$ and $p=5$, then $s(T)=\{1,2,4,4\}$ for $T=\{1,3,3,5\}$.  Thus, the cyclic group $C_p=\langle \gamma\rangle$ of order $p$ acts on the set $\mathcal{M}$ of multisets of size $k$ consisting of elements of $A$ via $\gamma\cdot T=s(T)$ for each $T\in\mathcal{M}$.  This action is without a fixed point (this is where the assumption $\gcd(k,p)=1$ is needed).  That is,
$$\big|\mathcal{M}/\!/C_p\big|=\frac{|\mathcal{M}|}{\left|C_p\right|}=\frac{1}{p}\,|\mathcal{M}|\,.$$
In fact, this argument also shows that the number of multisets in $\mathcal{M}$ whose element sum is congruent modulo $p$ to a fixed $t\in\{0,1,2,\ldots,p-1\}$ is $n_{p,k}$.)  
A: We ask  about the probability that  a multiset of size  $k$ drawn from
$[n]$ has  sum divisible by  $n$ where we  use the notation  from this
MSE link.
The exponential formula tells us that  the cycle index $Z(S_k)$ of the
unlabeled multiset operator $\mathfrak{M}$ on $k$ slots has OGF
$$Z(S_k) = [w^k]
\exp\left(\sum_{l\ge 1} a_l \frac{w^l}{l}\right).$$
The number of these multisets is thus given by
$$\left. Z(S_k)
\left(\sum_{q=1}^n z^q\right)\right|_{z=1}
= \left. [w^k]
\exp\left(\sum_{l\ge 1} 
\left(\sum_{q=1}^n z^{ql}\right) \frac{w^l}{l} \right)
\right|_{z=1}
\\ = [w^k]
\exp\left(n \sum_{l\ge 1} \frac{w^l}{l}\right)
= [w^k] \frac{1}{(1-z)^n} = {n+k-1\choose k}.$$
The desired probability is then given by
$${n+k-1\choose k}^{-1} \frac{1}{n}
\sum_{p=0}^{n-1} 
\left. Z(S_k)
\left(\sum_{q=1}^n z^q\right)\right|_{z=\exp(2\pi ip/n)}.$$
This is
$${n+k-1\choose k}^{-1} \frac{1}{n} 
\sum_{p=0}^{n-1} \left. [w^k]
\exp\left(\sum_{l\ge 1} \left(\sum_{q=1}^n z^{ql}\right)
\frac{w^l}{l}\right)\right|_{z=\exp(2\pi ip/n)}.$$
Evaluating the contribution for $p=0$ first we get
$${n+k-1\choose k}^{-1} \frac{1}{n} 
[w^k] \exp\left(\sum_{l\ge 1} n
\frac{w^l}{l}\right) = \frac{1}{n}.$$
This was  the same  as when  we counted.  It  remains to  evaluate the
contribution  from $1\le  p\le n-1.$  Now for  these $p$  if $l$  is a
multiple of $m = n/\gcd(p, n)$ we have
$$\sum_{q=1}^n \exp(2\pi ip/n)^{ql} = n.$$
We get zero otherwise. This yields for the remaining terms without the
scalar in front
$$\sum_{p=1}^{n-1} [w^k]
\exp\left(\sum_{l\ge 1} n
\frac{w^{ml}}{ml} \right)
= \sum_{p=1}^{n-1} [w^k]
\exp\left(\frac{n}{m} \sum_{l\ge 1} 
\frac{w^{ml}}{l} \right)
\\ = \sum_{p=1}^{n-1} [w^k]
\exp\left(\frac{n}{m}\log\frac{1}{1-w^m}\right)
= \sum_{p=1}^{n-1} [w^k] 
\frac{1}{(1-w^{n/\gcd(p,n)})^{\gcd(p,n)}}.$$
Using an Iverson bracket this becomes
$$\sum_{p=1}^{n-1} [[n/\gcd(p,n)|k]] \times
{k\gcd(p,n)/n + \gcd(p, n) - 1\choose k\gcd(p,n)/n}.$$
Putting it all together we thus obtain
$$\bbox[5px,border:2px solid #00A000]{
\frac{1}{n} +
{n+k-1\choose k}^{-1} \frac{1}{n} 
\sum_{p=1}^{n-1} [[n/\gcd(p,n)|k]] \times
{k\gcd(p,n)/n + \gcd(p, n) - 1\choose k\gcd(p,n)/n}.}$$
Note that $n/\gcd(p,n)$ is  a divisor of $n$ that is  at least two (we
would need $p=n$  to get $n/\gcd(p,n) = 1$ but  $p\lt n$).  This means
when $\gcd(k, n) =  1$ the Iverson bracket fails for  all $p$ and only
the first term remains. 
We  also  have a  Maple  routine  to  compute  and verify  these  data
including two routines for enumeration, the second one optimized.

with(numtheory);
with(combinat);

CF := (n,k) ->
add(coeftayl(1/(1-w^(n/gcd(p,n)))^gcd(p,n), w=0, k), p=1..n-1);

P := (n,k) -> 1/n+1/n/binomial(n+k-1,k)*CF(n,k);

# enumeration routine A, translates problem definition
ENUMA := proc(n, k)
local gf, terms, res, term;
option remember;

    if n=1 then return 1 fi;

    gf := mul(1/(1-u*z^q), q=1..n);
    terms := coeftayl(gf, u=0, k);

    res := 0;

    for term in terms do
        if degree(term) mod n = 0 then
            res := res + lcoeff(term);
        fi;
    od;

    res/binomial(n+k-1,k);
end;

# enumeration routine B, memory allocation optimized
ENUMB := proc(n, k)
local iter, res;
option remember;

    if n=1 then return 1 fi;

    res := 0;

    iter :=
    proc(count, sval, pos)
    local nxt;

        if count = k then
            if sval mod n = 0 then
                res := res + 1;
            fi;
            return;
        fi;

        if pos > n then return fi;

        for nxt from 0 to k - count do
            iter(count+nxt, sval+pos*nxt, pos+1);
        od;
    end;

    iter(0, 0, 1);
    res/binomial(n+k-1,k);
end;

P2 := (n, k) -> 1/n + 1/binomial(n+k-1,k)/n
*add(`if`(k mod n/gcd(p,n)  = 0,
          binomial(k*gcd(p,n)/n+gcd(p,n)-1, k*gcd(p,n)/n),
          0),
     p=1..n-1);

