how to find cubic polynomial that an unknown subset of a set of integers satisfies I have a set, $S$, of positive integers and I have reason to believe that some infinite subset of them may be parametrized by a cubic polynomial with integer coefficients evaluated at integer arguments. I have computed a finite subset of $S$, call it $S_{1}$ (all elements of $S$ less than some fixed bound), and would like to use $S_{1}$ to try to find such a cubic polynomial.
Of course, there is a brute force method: for each 4-tuple in the set, $S_{1}$, I could apply Lagrange interpolation to fit the 4-tuple and then see if other integers in $S_{1}$ are also values of this cubic polynomial. However, this is not practical if the size of $S_{1}$ is say 1000 or so (in which case, there are over 41 billion such 4-tuples).
So my question is whether there exist more efficient techniques for finding such a polynomial.
Being able to do the same for polynomials of larger degree would also be useful for me, as here the above brute force method becomes even less practical.
Does anyone know of such techniques, have relevant references,...?
 A: In a comment to an earlier incorrect answer that I posted (and then deleted), Sum One made the following clarification.

Here is a concrete example of the general problem. Suppose I create 950 random numbers between 1 and 10,000,000 and then I add to this set 50 numbers of the form $P(n)$ for $n=1,2,\ldots$ all less than 10,000,000 for some cubic polynomial $P(n)$. I give you this set $S$ of 1000 numbers (without any further information), tell you that some of them (at least 10, say, to avoid trivial solutions) are $P(0),P(1),\ldots$ for some cubic polynomial and ask you to find this polynomial.

What follows is due to various colleagues of mine.  Here is an algorithm that is cubic in the size of $S$. Write $P(n) = p_3n^3+p_2n^2+p_1n+p_0$ for some unknown integers $p_3, p_2, p_1, p_0$.  For distinct $a,b,c$, we have the identity
$${(b-c)P(a) + (c-a)P(b) + (a-b)P(c) \over (b-c)(c-a)(a-b)} = -(p_3(a+b+c) + p_2).$$
Now set $a=1$ and exhaust over all $1000$ possible choices from $S$ for $P(a)=P(1)$.  For each such choice, exhaust over all $999\times 998$ choices of two other numbers from $S$, and form two separate lists $L_1$ and $L_2$.  For $L_1$, hypothesize that the two other chosen numbers are $P(12)$ and $P(45)$, and use the above identity to compute what the value of $-((1+12+45)p_3+p_2)$ = $-58p_3-p_2$ would be if the hypothesis happened to be correct.  For $L_2$, hypothesize that the two other chosen numbers are $P(23)$ and $P(34)$, and similarly compute what the value of $-((1+23+34)p_3+p_2)$ = $-58p_3-p_2$ would be if this hypothesis were correct.  Then combine these two lists and sort them to look for occurrences of the same number in both lists; random collisions should be rare, so there is now a small list of candidates for 5-tuples $(P(1),P(12),P(23),P(34),P(45))$ that can be interpolated and and tested over the whole dataset.
Actually I lied a bit.  For what I just described to work, I only used the fact that $12+45=23+34$ and I didn't need all five numbers to be congruent modulo $11$.  The reason that I chose them to be all congruent modulo $11$ is that $P(1)$, $P(12)$, $P(23)$, $P(34)$, $P(45)$ must all be congruent modulo $11$ since $P$ has integer coefficients.  Therefore I can partition $S$ into congruence classes modulo $11$ and perform the above computation separately on each congruence class.  The outer exhaust will then loop over only about $100$ choices for $P(1)$ and the lists $L_1$ and $L_2$ will be only about ten thousand long rather than a million long.  Of course we have to do $11$ such computations but it is still a win.
Congruences can be exploited in another way.  The following alternative algorithm is asymptotically slower but uses less memory and should be more effective for the stated parameter sizes.  We are going to have four nested loops, guessing the values of $P(1)$, $P(50)$, $P(27)$, and $P(12)$ respectively. Given a hypothesized value for $P(1)$, we know that $P(50)\equiv P(1) \pmod{49}$, so the possible values for $P(50)$ are cut down by a factor of roughly $49$.  At the next level, we must have $P(27)\equiv P(1) \pmod{26}$ and $P(27)\equiv P(50)\pmod{23}$, providing further cutdown, etc. The sequence $1,50,27,12$ should provide the greatest expected cutdown, generating on average fewer than 100000 interpolations instead of 41 billion.
A: I don't see how this can be done exhaustively for a large number of integers, but I can suggest an approach that may be useful in practice.  First sort your $N$ intergers as $x_1\lt x_2\lt\cdots\lt x_N$, and make an inverse index so that you can quickly test of a given integer is in the set.  Now choose a smaller integer $n$.  For each 4-tuple $i_1\lt i_2\lt i_3\lt i_4$ such that $i_4\le i_1+n$, try the cubic polynomial $p$ for which $p(1)=x_{i_1},\ldots,p(3)=x_{i_3}$. Trying it means computing $p$ for some other arguments to see if those values are in set too (that's what the inverse index is for).  This takes time $O(n^3 N)$.  The idea is to do it for as large an $n$ as you can manage, then you will find all cubic polynomials with the property that some sequence of four contiguous values are not interleaved by too many other values.
