Product of sums of all subsets mod $k$? Input: set of $N$ integers and $k$ value (modulo)
Output: product of sum of all subsets of input set (modulo $k$)
e.g. for input {1,2,3} and $k = 20$, we have 1 * 2 * 3 * (1+2) * (1+3) * (2+3) * (1+2+3) = 2160.
and:
$2160 \equiv 0 \pmod{20}$
Output is 0.
Could anyone suggest if you can implement such an algorithm optimally? (in polynomial time - if this is possible).
 A: Here's a solution.  I'll edit later with the explanation.

Input: a set $S = \{s_1, s_2, \ldots, s_N\}$, and an integer $k$.
Define matrices $P,Q$ of size $(N+1) \times k$ initializing the first row
$P_{0,i} = i$, $Q_{0,i} = 1$
and filling in successive rows by settting
$$P_{n,i} = P_{n-1,i} \cdot P_{n-1, i+s_n}$$
$$Q_{n,i} = P_{n-1,i} \cdot Q_{n-1, i+s_n} + P_{n-1,i} \cdot Q_{n-1, i+s_n}$$
while reducing everything (including indices) modulo $k$.  
Output $Q_{N,0}$.

Filling in each entry of the matrix is just multiplication, addition, and reducing mod $k$, each of which I believe can be done in time $\log k$, so in total the algorithm is of order $Nk\log k$.
I implemented this in Python:
def aur(S, k):
    P = [list(range(k))]
    Q = [[1]*k]
    for n in range(1,len(S)+1):
        P.append([])
        Q.append([])
        for i in range(k):
            P[n].append( (P[n-1][i] * P[n-1][(i+S[n-1]) % k]) % k   )
            Q[n].append(  ( P[n-1][i] * Q[n-1][i] + Q[n-1][i]* P[n-1][(i+S[n-1]) % k]  ) % k )
    return Q[n][0]

You can check that e.g., aur$([1,2,3], 20) = 0$ and aur$([1,2,3], 3000) = 2160$ as in the OP's examples.  
Edit:  Here is the explanation.  
Define a polynomial $$P_S(x) = \prod_{T \subseteq S} (x + \sum_{i \in T} i).$$  If we set $x=0$, we almost get what we want, except that the product also includes the empty set which (I'm assuming) should not be there.  Nevertheless, it will be convenient to include it.  From $P_S(x)$ we can still recover our answer: we first divide by $x$ to eliminate the empty set, and then set $x=0$.  Equivalently, we take $P'(0)$.
The reason we define this polynomial is because it obeys a nice recurrence relation.  If we set $S_n = \{s_1, \ldots, s_n\}$, then 
$$P_{S_n}(x) = P_{S_{n-1}}(x) P_{S_{n-1}}(x+s_n).$$  
This follows from the fact that a subset $T$ of $S_n$ is either a subset $T$ of $S_{n-1}$, or a union $T \cup \{s_n\}$ for $T \subseteq S_{n-1}$.  Furthermore, we can start with the initial condition $P_{\emptyset}(x) = x$.
Computing the coefficients of $P_{S}(x)$ this way is enough to compute our output, but this polynomial will have degree $2^N$, so this is not feasible.  Instead, we use the product rule, setting $Q_S(x) = P_S'(x).$  This gives  
$$Q_{S_n}(x) = P_{S_{n-1}}(x) Q_{S_{n-1}}(x+s_n) + Q_{S_{n-1}}(x) P_{S_{n-1}}(x+s_n).$$  
We can then output $Q_{S_N}(0)$ and reduce mod $k$.  In this way we don't need the full list of exponentially many coefficients of these polynomials, but we only need to evaluate them for $x = 0, 1, \ldots, k$.
A: Here is an attempt at a polynomial time solution: We can of course assume that the given integers $a_i$ are in the range $\{0,\ldots,k-1\}.$ Denote by $v_m$ the number of times $m$ appears as a sum modulo $k$, i.e.
$$ v_m = |\{ S \subset [N] : \sum_{j \in S} a_j \equiv_k m\}| \leq 2^N,$$
for $m =0,\ldots,k-1$. The required sum is zero if $v_0 > 1$, otherwise it is
$$ \prod_{m=1}^{k-1} m^{v_m} \pmod k,$$
which can be computed with $\sum_{m=1}^{k-1} \log(v_m) \leq k N$ multiplications modulo $k$, so it is enough to compute $\{v_m\}$. Now consider the polynomial
$$ P(x) = \prod_{i=1}^N (1+x^{a_i}) = \sum_{S \subset [N]} x^{\sum_{j \in S} a_j}.$$
Then by taking $w_k = e^{-2 \pi i/ k}$, we can compute
$$ X_\ell = P(w_k^\ell) = \sum_{m=0}^{k-1} v_m w_k^{\ell m}.$$
$\{X_\ell\}$ is the discrete Fourier transform of $\{v_m\}$ so we just need to compute the inverse.
