Truncated Exponential Series Modulo $p$: Deeper meaning for a Putnam Question. Apparently B6 of the Putnam this year asked:

Suppose $p$ is an odd prime. Prove that for $n\in \{0,1,2...p-1\}$, at least $\frac{p+1}{2}$ of the numbers $\sum^{p-1}_{k=0} k! n^{k}$ are not divisble by $p$.  

With some rearrangements, this is equivalent to showing that $$E_p(z):=\sum_{k=0}^{p-1} \frac{z^k}{k!}$$ has at most $\frac{p-1}{2}$ zeros.  A proof of this is at the end.
My question is: Can we improve the bound for the number of zeros?  Also is there a deeper connection here with other parts of mathematics motivating this problem?

Proof of problem:  Consider  $$Q(z)=z^{p}-z+\sum_{l=0}^{p-1}\frac{z^{l}}{l!}.$$ Then for each integer $Q(n)=E(n).$ However, $$Q^{'}(z)\equiv E^{'}(z)-1=E(z)-\frac{z^{p-1}}{(p-1)!}-1\equiv E(z)+z^{p-1}-1.$$ Then, if $Q(n)=0$  for $n\neq0$ , we must also have $Q^{'}(n)=0$  so that $n$ is a double root of $Q(n).$ Since $\deg Q(n)=p$, we see that at most half of the integers $n\in\{ 1,2,\dots,p-1\}$ satisfy $E(n)=0.$ Since $E(0)=1$, we conclude the desired result.

Remark:  This was asked in a slightly different form on math stack exchange.  I felt the answer I posted there was inadequate there, and  I personally became more curious while attempting to answer the question.
 A: Thinking about this problem again, I found the following simpler explanation which does not invoke orthogonal polynomials and uses the nonvanishing of only one Hankel determinant.  The determinant arises via Padé approximants but the exposition below is self-contained.
Suppose in general that $a_0 \neq 0$, $a_1, a_2, \ldots, a_{q-1}$ are elements of a finite field $F$ of $q$ elements for which the polynomial
$$
A(X) = \sum_{i=0}^{q-1} \phantom. a_i X^i
$$
vanishes at all but $t$ nonzero field elements, say $x_1,x_2,\ldots,x_t$, with $2t < q-1$.  Then
$$
A(x) = \frac{P(x)}{Q(x)} (1 - x^{q-1})
$$
for some polynomials $P,Q$ of degree $t$, where
$$
Q(X) = \prod_{m=1}^t (X-x_m) = \sum_{i=0}^t \phantom. q_i X^i
$$
for some field elements $q_i$.  Thus $A(X)$ is within $O(X^{q-1})$ of the power series about $X=0$ of the degree-$t$ rational function $P(X)/Q(X)$.  For any $t' \in [t, \phantom. (q-1)/2)$ and $n \in (0, \phantom. q-1-2t')$ it follows that the square Hankel matrix
$$
(a_{n+i+j})_{i,j=0}^{t'} = \left(
\begin{array}{ccccc}
a_n & a_{n+1} & a_{n+2} & \cdots & a_{n+t'} \\
a_{n+1} & a_{n+2} & a_{n+3} & \cdots & a_{n+1+t'} \\
\vdots & \vdots & \vdots & & \vdots \\
a_{n+t'} & a_{n+t'+1} & a_{n+t'+2} & \cdots & a_{n+2t'}
\end{array}
\right)
$$
of order $t'+1$ is singular, because the nonzero column vector
$(0,0,\ldots,0,q_d,q_{d-1},q_{d-2},\ldots,q_1,q_0)^{\rm T}$ [with $t'-t$ initial zeros] is in the kernel.  Thus a single invertible matrix of this form implies that there are more than $t' \geq t$ nonzero $x\in F$ at which $A(x) \neq 0$.
To apply this to the Putnam problem, take $q=p$ and  $a_i = i!$, and set $t' = t = \frac12(p-1) - 1$ and $n=1$.  The resulting Hankel matrix $((i+j+1)!)_{i,j=0}^t$ is invertible mod $p$ thanks to the formula $\prod_{k=0}^t k!(k+1)!$ for its determinant.  [This formula can be obtained from properties of the Laguerre orthogonal polynomials, but also has an elementary direct proof: see the second solution of this problem at the Putnam directory.]  Hence $A(x) \neq 0$ for at least $t+1 = (p-1)/2$ nonzero values of $x \bmod p$, QED.
It is also known that the Hankel matrix $(a_{i+j})_{i,j=0}^t$ is invertible, else we'd have $P/Q \equiv P_1/Q_1 \bmod X^{2t}$ for some $P_1,Q_1$ of degree less than $t$, whence $P/Q = P_1/Q_1$ identically and $P/Q$ is not in lowest terms.  This lets us connect the above solution with the orthogonal polynomials evaluated at $X^{-1}$ that arose in our previous solution.  Indeed, define a bilinear pairing $\langle\cdot,\cdot\rangle$ on $F[X]$ by $\langle P,Q \rangle = I(PQ)$, where $I: F[X] \rightarrow F$ is the linear form taking each $X^i$ to $a_i$.  The restriction of this pairing to the polynomials of degree at most $t$ is nondegenerate because its Gram matrix is the Hankel matrix we just proved invertible.  But $G(X) := X^t Q(1/X) = \sum_{i=0}^t q_i X^{d-i}$ is orthogonal to $X^j$ for each $j = 1,2,\ldots,t+1$.  Hence $XG$ is orthogonal to every polynomial of degree at most $t$, and is therefore a multiple of the orthogonal polynomial of degree $t+1$ for our inner product (which is unique up to scaling thanks to the nondegeneracy of the pairing).  It follows that in this case the orthogonal polynomial of degree $t+1$ must vanish at zero and at $t$ other elements of $F$.  This does not happen often for classical orthogonal polynomials, but (as noted in my previous answer) there is at least one infinite family of examples, the Čebyšev polynomials of the second kind $U_{(q-1)/2}$ when $q \equiv 3 \bmod 4$.
A: This is quite a coincidence.  Some 25+ years ago I observed that a very similar result follows from classical formulas and properties for Laguerre polynomials, i.e. the orthogonal polynomials for the measure $e^{-x}dx$ on $[0,\infty)$, whose moments $k! = \int_0^\infty x^k e^{-x} dx$ are the coefficients of the polynomial in problem B6.  I thought at the time that this was a curiosity of very little interest because one expects such a random polynomial to have no more than say $O(\log^2 p)$ roots mod $p$.  Now this problem appears on the Putnam exam.  I was able to reconstruct and modify my argument to produce this solution, but wondered how anybody would be expected to find that under contest conditions.  Your solution is much likelier to be the intended one.
I haven't answered either of your questions yet.  I see that while I was typing this F.Voloch did answer both.  It's somewhat mysterious that my very different technique yields exactly the same bound; it generalizes to the reduction mod $p$ of moment-generating polynomials of other distributions with tractable orthogonal polynomials.  In that setting the bound can actually be sharp.  If we use the moments
$$
1, \phantom. 0, \phantom.
\frac14, \phantom. 0, \phantom.
\frac18, \phantom. 0, \phantom.
\frac5{64}, \phantom. 0, \phantom.
\frac7{128}, \ldots
$$
of $(2/\pi) \sqrt{1-x^2} \phantom. dx$ on $(-1,1)$, and subtract $2$ from the $x^{p-1}$ coefficient of the resulting polynomial
$$
1 + \frac12 \frac{x^2}{2} + \frac12 \cdot \frac34 \frac{x^4}{3} + \frac12 \cdot \frac34 \cdot \frac56 \frac{x^6}{4} + \cdots \pm 2x^{p-1}
$$
mod $p$ (which is OK because the argument does not use this leading coefficient), we get a polynomial with $t$ or $t+1$ roots according as $t$ is even or odd, namely $x=\pm1$ with multiplicity $1$, and each $x$ for which $1-x^2$ is a quadratic residue with multiplicity $2$; and indeed the corresponding orthogonal polynomials, which are Čebyšev polynomials of the second kind, are "tractable" for our purpose, but the relevant one $U_t$ has simple roots at all other nonzero $x \bmod p$ — and also at $x=0$ when $t$ is odd, which did not happen in Putnam B-6 and explains how an extra root can appear here.
For the B-6 polynomials, G.Myerson reports that the paper Voloch cited actually gets a bound $O(p^{2/3})$ on the number of roots mod $p$, which is much lower than $p/2$ but still well above what we expect to be true.  Here's some numerical evidence: the gp code
B6(p) = poldegree(gcd(Mod(1,p)*(x^p-x), Mod(1,p)*sum(k=0,p-1,x^k*k!)))
forprime(p=3,200,print([p,B6(p)]))

finds only one $p<200$ for which the polynomial has as many as $5$ roots mod $p$, namely $p=151$, and only three each for $4$ roots ($p=37, 97, 167$) or $3$ (these being $53,191,199$).
P.S. This would make another good example for this MO question (#69737: Contest problems with connections to deeper mathematics).
A: I've seen this trick in a paper of Mit'kin, Math Zametki 1992. There he improves the bound. This is related to Stepanov's method to bound the number of solutions of equations over finite fields.
