Aquí tienes una solución. Editaré más tarde con la explicación.
Entrada: un conjunto $S = \{s_1, s_2, \ldots, s_N\}$, y un entero $k$.
Define matrices $P,Q$ de tamaño $(N+1) \times k$ inicializando la primera fila
$P_{0,i} = i$, $Q_{0,i} = 1
y rellenando filas sucesivas mediante
$$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}$$
mientras se reduce todo (incluidos los índices) módulo $k$.
Salida $Q_{N,0}$.
Rellenar cada entrada de la matriz es solo multiplicación, adición y reducción módulo $k$, cada una de las cuales creo que se puede hacer en tiempo $\log k$, por lo que en total el algoritmo es del orden $Nk\log k$.
Implementé esto en 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]
Puedes comprobar que por ejemplo, aur$([1,2,3], 20) = 0$ y aur$([1,2,3], 3000) = 2160$ como en los ejemplos del autor original.
Editar: Aquí está la explicación.
Define un polinomio $$P_S(x) = \prod_{T \subseteq S} (x + \sum_{i \in T} i).$$ Si establecemos $x=0$, casi obtenemos lo que queremos, excepto que el producto también incluye el conjunto vacío que (asumo) no debería estar allí. Sin embargo, será conveniente incluirlo. A partir de $P_S(x) todavía podemos recuperar nuestra respuesta: primero dividimos por $x$ para eliminar el conjunto vacío y luego establecemos $x=0$. De forma equivalente, tomamos $P'(0)$.
La razón por la que definimos este polinomio es porque obedece una bonita relación de recurrencia. Si establecemos $S_n = \{s_1, \ldots, s_n\}$, entonces
$$P_{S_n}(x) = P_{S_{n-1}}(x) P_{S_{n-1}}(x+s_n).$$
Esto se deduce del hecho de que un subconjunto $T$ de $S_n$ es o bien un subconjunto $T$ de $S_{n-1}$, o una unión $T \cup \{s_n\}$ para $T \subseteq S_{n-1}$. Además, podemos empezar con la condición inicial $P_{\emptyset}(x) = x$.
Calcular los coeficientes de $P_{S}(x)$ de esta manera es suficiente para calcular nuestra salida, pero este polinomio tendrá un grado $2^N$, por lo que no es factible. En su lugar, usamos la regla del producto, estableciendo $Q_S(x) = P_S'(x).$ Esto da
$$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).$$
Luego podemos obtener $Q_{S_N}(0)$ y reducir módulo $k$. De esta manera, no necesitamos la lista completa de coeficientes exponencialmente grandes de estos polinomios, sino que solo necesitamos evaluarlos para $x = 0, 1, \ldots, k$.