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.