# The Prime Polynomial : Generating Prime Numbers

First of all, i'll confess i'm no math geek. I'm from Stackoverflow, but this question seemed more apt here, so i decided to ask you guys :)

Now, i know noone has discovered (or ever will) a Polynomial that generates Prime Numbers. But i've read about Curve Fitting (or Polynomial Fitting) so i was wondering if there was a way, we could have a simple n-degree Polynomial that could generate the first 1000 (or X) primes accurately.

I don't need it to generate all the primes, maybe just upto 1 million, since we already have the data, can we deduce that polynomial?

How big will be the polynomial for it to be accurate? Could you give an example for the first 100 primes? Am i just plain naive?

-
By Lagrange interpolation one can fit a polynomial through any finite set of values. Doing it for 1000 primes would not be helpful as it would have degree 999 and huge coefficents. – Robin Chapman Oct 14 '10 at 6:39
@Robin, so it can't be done with degree less than 999? :( – st0le Oct 14 '10 at 6:48
Of course there may be some other polynomial that does the trick, but it's not clear how to find it. Here's a link that contains more examples (not for n=100 though): mathworld.wolfram.com/Prime-GeneratingPolynomial.html – Hans Lundmark Oct 14 '10 at 6:53
@Robin, clearly the degree will be at most 999, but why not less? – Hans Oct 14 '10 at 7:42
@Hans: If we need to generate the first 1000 primes in order (i.e P(n) is the nth prime), there is a unique polynomial of degree <= 999, and that could be found by guessing the degree, finding it using interpolation and verifying for the remaining values... – Aryabhata Oct 14 '10 at 14:27

While Robin was commenting, I cooked up a Mathematica one-liner which gives you the Lagrange interpolation polynomial for the first n primes:

f[n_] := Sum[Prime[k] Product[x-m,{m,Join[Range[1,k-1], Range[k+1,n]]}] / Product[k-m,{m,Join[Range[1,k-1], Range[k+1,n]]}] ,{k,1,n}] // Expand


And I don't think anyone would call the results "simple". For n=100 it's too big to paste here, but n=15 gives

4748 - (1752758563*x)/120120 + (71043957851*x^2)/3783780 - (4320411427*x^3)/316800 + (378496362427*x^4)/59875200 - (7239131749*x^5)/3628800 + (11528263*x^6)/25920 - (2075348983*x^7)/29030400 + (5090997277*x^8)/609638400 - (977071*x^9)/1382400 + (743507*x^10)/17418240 - (568871*x^11)/319334400 + (2207*x^12)/45619200 - (3151*x^13)/4151347200 + (89*x^14)/17435658240

-
Thanks Hans! Kind of what i was going for...guess storing the numbers directly is much more efficient. I'll wait till the end of the day for more answers and accept yours if there's none better. – st0le Oct 14 '10 at 7:02
On the other hand, InterpolatingPolynomial[Prime[Range[n]],x] will give the Newton form of the interpolating polynomial. Still, as mentioned by numerous people in this thread, the relatively complicated structure of the primes makes this construct impractical. – J. M. Oct 17 '10 at 15:50
@J. M.: Nice! I didn't know there was a builtin function for that. Of course, wrapping your expression in Expand[ ... ] gives the same result as mine, but in a much cleaner way. – Hans Lundmark Oct 17 '10 at 18:39

From Mathworld

However, there exists a polynomial in 10 variables with integer coefficients such that the set of primes equals the set of positive values of this polynomial obtained as the variables run through all nonnegative integers, although it is really a set of Diophantine equations in disguise (Ribenboim 1991). Jones, Sato, Wada, and Wiens have also found a polynomial of degree 25 in 26 variables whose positive values are exactly the prime numbers (Flannery and Flannery 2000, p. 51).

Unfortunately, the primes do not come out in order, so this will not help for what you want. But it is interesting.

-