diff options
| author | Charles Harris <charlesr.harris@gmail.com> | 2014-10-04 12:37:11 -0600 |
|---|---|---|
| committer | Charles Harris <charlesr.harris@gmail.com> | 2014-10-04 12:37:11 -0600 |
| commit | 0650d04078666997c3dac3beec6fe872dca797dd (patch) | |
| tree | cce267dbba6b88798ad9125d9c7669b493271504 /numpy/polynomial/hermite.py | |
| parent | 58350f4608a22f4b4b66795f51eaefc206bd02b8 (diff) | |
| download | numpy-0650d04078666997c3dac3beec6fe872dca797dd.tar.gz | |
MAINT: Improve computation of scaled companion matrices.
The previous method used for hermite and hermite_e polynomials suffered
from double overflow for polynomials of large degree. Those numbers were
later scaled down by equally large numbers, but the result was NaN. The
wanted values are now computed in such a way that overflow in some
entries is replaced by underflow in others. The resulting zeros are a
negligible perturbation of the companion matrix.
Diffstat (limited to 'numpy/polynomial/hermite.py')
| -rw-r--r-- | numpy/polynomial/hermite.py | 6 |
1 files changed, 3 insertions, 3 deletions
diff --git a/numpy/polynomial/hermite.py b/numpy/polynomial/hermite.py index 1fd49d774..02b7fb53d 100644 --- a/numpy/polynomial/hermite.py +++ b/numpy/polynomial/hermite.py @@ -1577,13 +1577,13 @@ def hermcompanion(c): n = len(c) - 1 mat = np.zeros((n, n), dtype=c.dtype) - scl = np.hstack((1., np.sqrt(2.*np.arange(1, n)))) - scl = np.multiply.accumulate(scl) + scl = np.hstack((1., 1./np.sqrt(2.*np.arange(n - 1, 0, -1)))) + scl = np.multiply.accumulate(scl)[::-1] top = mat.reshape(-1)[1::n+1] bot = mat.reshape(-1)[n::n+1] top[...] = np.sqrt(.5*np.arange(1, n)) bot[...] = top - mat[:, -1] -= (c[:-1]/c[-1])*(scl/scl[-1])*.5 + mat[:, -1] -= scl*c[:-1]/(2.0 *c[-1]) return mat |