Rework the *_roots functions and document them.

As I indicated in a [comment](https://github.com/scipy/scipy/pull/3795#issuecomment-48988505) to gh-3795. The *_roots routines need a good going over.  Here I've reworked them and documented them (although they don't appear in generated index). They are also comprehensively tested now.  This is also somewhat faster.  An example, before:
```
In [1]: from scipy import integrate, special

In [2]: f = lambda x: 1.0 / (1 + 16.0 * x**2)

In [3]: timeit integrate.quadrature(f, -1, 1, tol=1e-14, rtol=1e-14)
10 loops, best of 3: 34.5 ms per loop

In [4]: timeit special.t_roots(25)
1000 loops, best of 3: 397 µs per loop

In [5]: timeit special.p_roots(25)
1000 loops, best of 3: 507 µs per loop
```
After:
```
In [3]: timeit integrate.quadrature(f, -1, 1, tol=1e-14, rtol=1e-14)
100 loops, best of 3: 15.4 ms per loop

In [4]: timeit special.t_roots(25)
1000 loops, best of 3: 15.9 µs per loop

In [5]: timeit special.p_roots(25)
1000 loops, best of 3: 233 µs per loop
```
Maybe something like `verify_gauss_quadrature` or some such.