Sum of $k$ elements of $\{1,...,p\}$ divisible by $p$

Question

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$?

Answer 1

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}$.)

Answer 2

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);