This program solves a version of the subset-sum problem: given a set of non-negative integers \$S=\{a_1,a_2,\cdots,a_n\}\$ (they may not be distinct) and a non-negative integer \$s\$, decide whether there exists a subset of \$S\$ (let's call it \$T\$) such that \$\sum\limits_{b\in T}b=s\$.

Let the sum of all elements in \$S\$ be \$t\$ (\$t=\sum\limits_{j=1}^n a_j\$). Consider the polynomial $$f(x)=\prod\limits_{j=1}^{n}\left(1+x^{a_j}\right)$$ which can be expanded to the form $$f(x)=c_0+c_1x^1+c_2x^2+\cdots+c_tx^t\quad\left(t=\sum\limits_{i=1}^n a_i\right)$$ It is easy to see that the problem is equivalent to deciding whether \$c_s=0\$. Let \$m=\max\{s,t-s\}+1\$ and \$\zeta=e^{2\pi i/m}\$, then we have $$\sum\limits_{k=0}^{m-1}{\left(\zeta^p\right)}^k=\begin{cases} 0,&\text{if }m\text{ divides }k\\ m,&\text{if }m\text{ does not divide }k \end{cases}$$ where \$p\in\mathbb{Z}\$. Therefore we can isolate the coefficient \$c_s\$ in \$f(x)\$. Consider $$g(x)=x^{-s}f(x)=c_0x^{-s}+c_1x^{1-s}+\cdots+c_s+\cdots+c_tx^{t-s}$$We have $$\sum\limits_{k=0}^{m-1}g\left(\zeta^k\right)=\sum\limits_{j=0}^{t}c_j\sum\limits_{k=0}^{m-1}{\left(\zeta^k\right)}^{j-s}=mc_s$$ Therefore $$c_s=\sum\limits_{k=0}^{m-1}\zeta^{-ks}f\left(\zeta^k\right)$$ So we only need to compute \$c_s\$ and see if \$|c_s|>\frac{1}{2}\$ (we do not directly compare \$c_s\$ to zero since there are rounding errors during floating-point computation). When we compute \$f\left(\zeta^k\right)\$, we do not directly pass \$\zeta^k\$ as the argument; instead, we pass its argument \$2\pi ik/m\$, and then we can compute \${\left(\zeta^k\right)}^{a_j}\$ with the Euler formula $$e^{2\pi ia_j k/m}=\cos(2\pi a_j k/m)+i\sin(2\pi a_j k/m)$$ without directly computing the exponent of a complex number.

Here is the Python code:

import math
from typing import *

class SubsetSumSolver:
    def __init__(self, data: List[int], desired_sum: int):
        self.data = data # a
        self.desired_sum = desired_sum # s
        self.len_data = len(data) # n
        self.sum_data = sum(data) # t
    def evaluate_polynomial(self, theta: float) -> complex:
        result = 1.
        for element in self.data: # element: a_j
            arg = element * theta
            result *= 1 + math.cos(arg) + 1j * math.sin(arg)
        return result
    def solve(self) -> bool:
        if self.desired_sum > self.sum_data:
            return False
        if self.desired_sum == self.sum_data:
            return True
        if self.desired_sum == 0:
            return True
        order = max(self.desired_sum,
            self.sum_data - self.desired_sum) + 1 # m
        result = 0.
        double_pi = 2 * math.pi
        for k in range(order):
            theta = double_pi * k / order
            f_result = self.evaluate_polynomial(theta)
            theta_s = -double_pi * k * self.desired_sum / order
            term_s = math.cos(theta_s) + 1j * math.sin(theta_s)
            result += term_s * f_result
        result /= order
        return abs(result) > 0.5

a = [2, 3, 5, 6, 8, 13, 27, 36]
s = 67 # test data
solver = SubsetSumSolver(a, s)
print(solver.solve()) # True

I have a question here: If \$n\$ is very large, will it be possible that this code doesn't work because of huge rounding errors?