Approximating a magic function

I got to wondering about a function to generate the primes: $\Pi(n)=$ n'th prime. Obviously that would have truly weird behavior: an infinite number of primes are only 2 apart, and the gaps between primes can also be exceedingly large too. But just for laughs, suppose one existed.

I'd bet that you couldn't use a Taylor series to calculate it – maybe locally, but nothing like the famous $1 + x/1 + x^2/2! + x^3/3! ...$.

Suppose you approximated it with longer and longer polynomials. Name the polynomial that fits the first $N$ primes ${}_{N}\Pi(x)$ (with ${}_{N}\Pi(k) = p_k$; the k'th prime, with k less or equal to $N$), and the coefficient of the $x^j$ term call ${}_{N}c_j$.<\p>

As $N$ increases, and the new polynomial fits more and more primes, do the coefficients converge? The first (the constant term of the polynomial) ${}_Nc_0$ is always $1$. How about the second ${}_Nc_1$ (coefficient of $x$) and third ${}_Nc_2$ (coefficient of $x^2$)? (By the construction of these, the polynomial to fit $N$ primes will only have $N+1$ coefficients.)

It won't come as a great surprise to see that they don't seem to converge. The polynomials resulting from fitting the first 30 primes gives this for the behavior of those two coefficients. They look like they're about to blow up.

But its not that simple. Expand the graph to include the first 50 polynomial second and third coefficients, and they switch directions and start to blow up the other way. You see that the deviation that looked so large in the plot above is invisibly small in the one below.

Not a big surprise – we didn't expect that the magic $\Pi(x)$ function to find all the primes was going to be simple to approximate. After all, the magic function has to be extremely "jumpy" and polynomials are nice and smooth. But the variation is certainly dramatic.

Of course this isn't entirely fair – trying to fit polynomials to points is famously ugly and unstable. But this is pretty dramatic.

UPDATE: If you were wondering if I was plotting round-off error, the answer is no. I did the calculations using pari/gp 2.13.3, and only turned the integer rational numbers into floating point at the printing step. If you are curious, I include the script below:

Top = 50
coeffs=matrix(Top, Top)
for(N=2, Top, \
 target=primes(N) - vector(N, k, 1); \
 arr = matrix(N, N, i, j, i^j); \
 co = (1/arr)*target~; \
 for(i=1, N, coeffs[N,i] = co[i];););
\\
for(i=1,Top,\
 for(j=1, Top, print1(1.*coeffs[j,i],","););print(" ");)