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