Consider the following constraint satisfaction problem: Let $\alpha_1 , \ldots, \alpha_k \in \mathbb{R}$ be given as well as an error parameter $\epsilon$. Find $p_1, \ldots, p_n$ such that  

(i) $0 \le p_1, \ldots, p_n \le 1$

(ii) For $1\le j \le k$, $|(\sum_{i=1}^n p_i^j) - \alpha_j|\le \epsilon$.

I am interested in the time complexity of this problem: In particular, based on some symmetry considerations, one can show that it suffices to restrict one's search to the case where $p_1, \ldots, p_n$ all come from a set of size at most k. Based on this observation, one can consider all possible ways of partitioning $p_1, \ldots, p_n$ into k sets (which takes time $n^k$) and subsequently solve a problem which requires one to solve polynomial equations over the reals (in k variables and of degree bounded by k). This takes time $k^k$. The dependence of n^k is prohibitive for me and I was wondering if there is a way to solve this in time $O(n^{O(1)} \cdot k^{k})$ or at least better than $n^k$.