Next Edit:
Let's now construct (or, rather, prove the existence of) a non-zero polynomial $P$ with coefficients $0,\pm 1$ of degree at most $n$ that has at least $c\sqrt n$ roots on $(0,1)$. This construction will, finally, justify my casual remark in the beginning of this long discussion (which I still hope to continue) and show that my memory, albeit failing, can still be occasionally trusted somewhat. If you ask me when, where, and by whom this argument was first invented, I have no idea.
We'll proceed as before but do everything way more carefully. First, we shall show that if $a>0$ is a small number, then every polynomial $P(z)=\sum_{k=0}^n a_kz^k$ with coefficients $a_k\in [-1,1]$ and $a_0=\pm 1$ satisfies $\|P\|_{L^0(I)}\ge \exp(-C/a)$ where $I=[1-2a,1-a]$ and for a non-negative function $f$ on an interval $J$, $\|f\|_{L^0(J)}=\exp\left[\frac{1}{|J|}\int_J\log| f|\right]$ is the geometric mean of $f$ on $J$.
This is achieved by considering the domain $\Omega$ that is a disk centered at some point $z_0\in(0,\frac 13)$ (say, $z_0=\frac 16$) of radius $1-|z_0|$ with a slit $[1-3a,1]$ (the blue circle on the figure below).
Applying Jensen's inequality, we see that
$$
0 =\log|P(0)|\le\int_{\partial\Omega}\log|P|d\omega
$$
where $\omega$ is the harmonic measure on $\partial\Omega$ with respect to $0$. Now we first estimate the integral of $\log_+|P|$ using the trivial bound $|P(z)|\le\frac 1{1-|z|}$. We split $\partial\Omega$ into the slit part $S$ and the circle part $C$. Note that the conformal mapping of $\Omega$ to the unit circle is given by an explicit formula, so the density of $\omega$ can be found exactly, but I still prefer a back of the envelope computation that isn't algebraicly heavy even if it is a bit lengthier.
First, we consider the harmonic function $u(z)=\log\frac 1{|1-z|}$ in $\Omega$. we have
$$
0=u(0)=\int_S\log\frac 1{|1-z|}\,d\omega(z)+\int_C\log\frac 1{|1-z|}\,d\omega(z)\,.
$$
Since $\log\frac 1{|1-z|}\ge -\log 2$ on $C$, we conclude that
$$
\int_S\log_+|P|\,d\omega\le \int_S\log\frac 1{1-|z|}\,d\omega(z)
\\
=\int_S\log\frac 1{|1-z|}\,d\omega(z)\le\log 2\,.
$$
Now note that on $C$, the harmonic measure $\omega$ is dominated by the harmonic measure for the disk without a slit, which has bounded density with respect to the Lebesgue measure on the circumference, so the integral of $\log\frac{1}{1-|z|}$ with respect to $\omega$ is uniformly bounded by some constant depending on (fixed) $z_0$ but not on the size of the slit. Thus, $\int_{\partial\Omega}\log_+|P|\,d\omega\le C$ independently of $a$ and, therefore, $\int_I \log_-|P|\,d\omega\le \int_{\partial\Omega}\log_-|P|\,d\omega\le C$ as well.
Now we need a more clear idea of what $\omega$ is on $I$. First, map the disk conformally to the upper half-plane so that $0$ is mapped to $i$ and $1$ to $0$, say. Then the slit $S$ will be mapped to $[0,hi]$ for $h$ comparable to $3a$ and the mapping will be bi-Lipshitz on the slit and $\omega$ will be mapped to the harmonic measure $\omega'$ in the half-plane with the new slit.
Now, we use the usual $w\mapsto\sqrt{w^2+h^2}$. Then we'll get $d\omega'(w)\approx \frac{|w|}{\sqrt{|w^2+h^2|}}|dw|$. When $z$ is in $I$, the corresponding point $w$ is in the "middle part" of $[0,hi]$, so the factor in front of $|dw|$ is comparable to $1$. Thus, on $I$, we have $d\omega(z)\approx |dz|$ and we conclude that $\int_{I}\log_-|P(x)|\,dx\ge -C$, i.e., $\|P\|_{L^0(I)}\ge \exp[-C/|I|]=\exp[-C/a]$.
If $a_0=0$ and $a_k\in\{0,\pm 1\}$, consider the least $m\in[0,n]$ for which $a_m=\pm 1$. The factor $x^m\ge(1-2a)^n\ge \exp[-3an]$ on $I$, so in this case
$$
\|P\|_{L^0(I)}\ge\exp\left[-\tfrac Ca-3an\right]\,.
$$
We shall now fix $a=\frac 1{\sqrt n}$, so
$$
\|P\|_{L^0(I)}\ge\exp\left[-C\sqrt n\right]\,.
$$
This part has been definitely known to Borwein and Erdelyi. Now we play the same game as before but only on the interval $I$. Formally, for a polynomial $P$, we define
$$
\widetilde P(t)=P(1-a-at), \qquad t\in[0,1]\,.
$$
The above result states now that for every non-zero $P$ with coefficients $0,\pm 1$, we have
$$
\|\widetilde P\|_{L^0([0,1])}\ge \exp[-C\sqrt n]\,.
$$
Now, noting that $\|\widetilde p\|_{L^\infty([0,1])}\le n$ for any $p$ with coefficients $0,1$ and using the pigeonhole as before, we can find a non-zero polynomial $P$ with coefficients $0,\pm 1$ such that
$$
\left|\int_0^1 \widetilde P(t) t^\ell\,dt\right|\le n2^{-n/m},\qquad \ell=0,1,\dots,m-1\,.
$$
Assuming that $\widetilde P$ has only $m-1$ (or fewer) sign changes on $[0,1]$ (i.e., that $P$ has fewer than $m$ zeroes on $I$), we again construct $Q(t)=\prod_j(t-t_j)$ as before and note that $\widetilde PQ$ preserves sign on $[0,1]$. Also, $\|Q\|_{L^0([0,1])}\ge\exp[-Cm]$ (each factor $|t-t_j|$ has geometric mean uniformly bounded from below. So, we get the chain of inequalities
$$
\exp[-C\sqrt n-Cm]\le \|\widetilde P Q\|_{L^0([0,1])}\le\|\widetilde P Q\|_{L^1([0,1])}
\\
=\left|\int_0^1 \widetilde P Q\right|\le 2^m n2^{-n/m}\,,
$$
and when $m=c\sqrt n$ with small enough $c>0$, we get a contradiction.