A symmetric matrix $A$ has an eigenvalue in $(-1,1)$ iff $A^2$ has an eigenvalue in $[0,1]$. Equivalently this means that
$$
\min_{\Vert v\Vert=1} \frac{\Vert Av\Vert^2}{\Vert v\Vert^2} <1.
$$
Hence you need to decide if the minimum of a nonnegative quadratic function on the unit sphere is $<1$.



I can think of two ways of doing this  but I cannot comment on the computational complexity.

The first is a Monte Carlo type method. Sample $v$ uniformly on the sphere and compute $\Vert A v\Vert$. If the sample is large  and there exists an eigenvalue  of $A^2$ that is $<1$ you can pick detect this with high confidence. Technically, if the size of of the matrix is $n\times n$, then choose the coordinates  $v_1,\dotsc, v_n$  of the random vector $v$ to  be independent  standard normal random variables  and test if
$$
\Vert A v\Vert^2 <\Vert v\Vert^2.
$$

Here's a justifications. One can choose an orthonormal basis  of your ambient space consisting of eigenvectors  $e_1,\dotsc, e_n$ of $A^2$ corresponding to eigenvalues $\lambda_1\leq \cdots \leq \lambda_n$. If $\lambda_1<1$ then $\Vert Av\Vert^2<1$ in a neighborhood $U$ of $e_1$ on the unit sphere $\{\Vert v\Vert=1\}$.  The probability that a random vector  $v$ lands in $U$ is proportional to the size of $U$. One can  provide lower bounds on the volume of $U$  (I won't do it here) that will give give you a lower bound on the probability of a random $v$ landing in $U$. That will give you and idea of the size of of the sample. I cannot comment on the computational complexity of this approach.

Your samples yield only numbers $\Vert Av\Vert$ substantially bigger than $1$ then you can say with high confidence that $A$ has no small eigenvalues. If the numbers $\Vert A v\Vert$ are bigger than $1$  "many" are close to $1$ there is a bit of ambiguity.  

The second approach is  by gradient descent. Consider the   function
$$
f:\{\Vert v\Vert=1\}\to [0,\infty),\;\;f(v)=\Vert Av\Vert^2.
$$
A  flow line  $v(t)$ of the negative gradient flow
$$
\frac{dv}{dt}=-\nabla f(v)
$$ 
will converge exponentially to an eigenvector of $A^2$.  If the initial condition $v(0)$ is uniformly random on the unit sphere, then, with probability $1$, this flow line will converge exponentially to an eigenvector corresponding to a minimal eigenvalue.

Solve numerically this equation with random initial condition $v(0)$. I speculate that this discretization will  lead you very fast to a decision concerning small e-values.  The speed of convergence depends on how "packed" are the eigenvalues of $A^2$: if they are all packed in a small interval, the convergence will be slower.

Again, I cannot comment on the computational complexity of this approach.