BUG: special: fix overflow at low order in eval_hermite

Hermite polynomials are very large and using the relationship to
the laguerre polynomials overflow earlier than necessary (around
N = 67 or so.)  Using the Clenshaw recursion we can evaluate to
much higher orders.

From now we evaluate the hermite polynomial as a stretched and
scaled normalized hermite (He) polynomial since He is slightly
smaller at any given order.
BUG: special: improve h_roots accuracy at high order

The roots of the hermite polynomials are still computed by using the
Golub-Welsch algorithm.  However a single pass of newton's method is
applied to improve the roots.  The symmetry of the Hermite roots
about zero is used as well.  The weights are calculated and then
rescaled to ensure that the integral of 1 is correct.

These changes were largely copied from the hermgauss implementation
in numpy (by Charles Harris).

Fixes #1991
Merge pull request #2806 from ewmoore/hroot

BUG: special: eval_hermite overflows early and h_roots are inaccurate at high order

Fixes gh-1991