Product of sums of all subsets mod $k$?
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$.
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.