Quadrature speed-up: scipy.special.orthogonal.p_roots with cache

The p_roots now caches its results. This is because scipy.integrate.quadrature, which makes use of this function, calls that function for basically the same inputs (integer values), over and over again. In my case, that single function was way more expensive than the function being integrated.

With the modification, I easily get 100-fold of speed up for scipy.integrate.quadrature. Here is the benchmark:

```python
import numpy as np
from scipy.integrate import quadrature

def link(x):
    return 1.0 / (1.0 + np.exp(-x))

def gauss(x, mu, sig2):
    return np.exp(- (x - mu)**2 / sig2)

class Cases(object):
    def __init__(self):
        np.random.seed(1)
        
        self._n_small = 100
        self._mu_small = np.random.rand(self._n_small) - 0.5
        self._sig2_small = np.random.gamma(1, 0.5, self._n_small)

        self._n_big = 1000
        self._mu_big = np.random.rand(self._n_big) - 0.5
        self._sig2_big = np.random.gamma(1, 4.0, self._n_big)

    def run_small(self):
        for i in xrange(self._n_small):
            std = np.sqrt(self._sig2_small[i])
            a = self._mu_small[i] - 6 * std
            b = self._mu_small[i] + 6 * std
            quadrature(lambda x: x * link(x) * gauss(x, self._mu_small[i], self._sig2_small[i]), a, b)

    def run_big(self):
        for i in xrange(self._n_big):
            std = np.sqrt(self._sig2_big[i])
            a = self._mu_big[i] - 6 * std
            b = self._mu_big[i] + 6 * std
            quadrature(lambda x: x * link(x) * gauss(x, self._mu_big[i], self._sig2_big[i]), a, b,
                       miniter=100, maxiter=101)
```
And here are the results comparing scipy.integrate.quadrature using the previous and new p_roots:

| Case  | previous          | new           | speedup  |
|---------|--------------|---------------|--------------|
| run_small | 1.28 sec  | 137 msec  | x9           |
| run_big     | 21.6 sec  | 107 msec  | x201       |


The run_big case is the kind of setting I'm working with: I start with a miniter that is high enough for the quadrature to converge in one iteration.