$p$ is an odd prime number,$S = \{x \mid 1 \leq x \leq 2p, x \in \mathbb{Z} \}$,

$A$ is a subset of $S$, satisfying

  1. $\operatorname{card}(A) = p$

  2. $\sum\limits_{x\in A} x \equiv 0 \pmod p$

For a given $p$, how many $A$ are there? How can I define a function in Mathematica, that given a prime $p$ will return the count of the $A$?

  • 2
    $\begingroup$ Phew. I don't think that this is an "easy" problem. While the straight-forward k = With[{sets = Subsets[Range[2 p], {p}]}, Length[sets] - Total[Unitize[Mod[sets.ConstantArray[1, p], p]]] ] will do it for tiny primes p, my computer's memory usage is maxed out already with p = 17. $\endgroup$ Mar 25, 2018 at 7:58
  • 5
    $\begingroup$ Have you seen this? $\endgroup$ Mar 25, 2018 at 8:05

1 Answer 1


The closed form solution in OEIS (see J.M.'s link) is obviously the simplest way to do this. If you want to play with it, an idea that occurred to me is to use a bit of generatingfunctionology.

Consider the product $$ G(p,t,x)=\prod_{k=1}^{2p}(1+tx^k). $$ This product contains a term for each subset $A\subseteq\{1,2,\ldots,2p\}$. The idea is that if $k\in A$, we pick the term $tx^k$, and when $k\notin A$, we picked the term $1$. So the subset $A=\{k_1,k_2,\ldots,k_n\}$, $k_1<k_2<k_3<\cdots<k_n$ yields the term $t^nx^{k_1+k_2+\cdots+k_n}$, and $$ G(p,t,x)=\sum_{A\subseteq \{1,2,\ldots,2p\}}t^{|A|}x^{\sum_{k\in A}k}. $$ In other words, the exponent of the variable $t$ contains the number of terms $k_i$, and the exponent of the other variable $x$ contains their sum.

  • You are interested only in those terms where the exponent of $t$ is exactly $p$. This is because you wanted $|A|=p$.
  • From those terms you want to further throw away the terms where the exponent of $x$ is not divisible by $p$. Your answer is then the sum of the coefficients of the remaining terms (each acceptable $A$ contributes $+1$ to one such term).
  • The largest acceptable exponent of $x$ is $\sum_{k=p+1}^{2p}k=p(1+3p)/2$.


In[1]:= G[p_]:=Coefficient[Expand[Product[(1 + t x^k), {k, 1, 2 p}]], t, p]
In[3]:=Sum[Coefficient[Out[2], x, 19*n], {n, 1, (1+3*19)/2}]

calculates the case $p=19$. You may want to code $p$ as a variable in that last formula also.

Checking this number with the OEIS formula adds to our confidence in the correctness of the method:

In[4]:=(Binomial[38, 19] + 2*18)/19

I tested this up to $p=47$. My Mathematica had no trouble whatsoever handling the generating function $G(47,t,x)$, and also gave the correct answer $34589385964790856394651396$. Took may be a second to compute (wasn't instantaneous). With $p=197$ calculating G[p] took in the ballpark of 2 minutes. I don't think I want to find the limit given that this method is inferior to the closed form. But, I was kind impressed that Mathematica handles the generating function easily enough - as long as you have the sense not to print it on the screen:-)

The difference to the method of iterating over the subsets of $\{1,2,\ldots,2p\}$ (see Henrik's comment under the question) is that organizing the calculation this way essentially recursively calculates the subset sums of subsets with largest element $\le \ell$, and only stores the relevant pieces of information in the partial generating function $$ G_\ell(p,t,x):=\prod_{k=1}^\ell(1+tx^k)\in \Bbb{Z}[t,x]. $$ Here the maximal $t$-degree of such a function is obviously $\ell\le2p$. The maximal $x$-degree is bounded from above by $\sum_{k=1}^{2p}k\approx 2p^2$. Therefore the number of terms in $G_\ell(p,t,x)$ is upper bounded by a cubic polynomial in $p$. It shouldn't be too difficult to see that the number of digits in any coefficient of $G_\ell(p,t,x)$ is upper bounded by a polynomial of $p$. The conclusion is that the calculation of $G(p,t,x)=G_{2p}(p,t,x)$ takes both polynomial time and polynomial space. This shows in my testing. A method requiring an exponential number of iterations will not unnaturally make any program to balk sooner than what happens with this method.

  • $\begingroup$ Very nice! ${}$ $\endgroup$ Mar 25, 2018 at 16:51
  • $\begingroup$ I can't test at the moment, but is it any faster if you use G[p_] := SeriesCoefficient[Product[(1 + t x^k), {k, 1, 2 p}], {t, 0, p}]? Avoiding Expand[] might help reduce the memory consumed. $\endgroup$ Mar 25, 2018 at 19:13
  • $\begingroup$ @J.M. I've been waiting for seven minutes now for G[197] to complete, so looks like that is a slower way. Don't know about memory consumption :-( I guess your idea was partly to automatically throw out $t^\ell$-terms when $\ell>p$. Should help, but may be the calculation is not organized optimally. $\endgroup$ Mar 25, 2018 at 19:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.