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