Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fix coefficients in expm helper function #9705

Merged
merged 1 commit into from Apr 19, 2019
Merged

Fix coefficients in expm helper function #9705

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
8 changes: 3 additions & 5 deletions scipy/sparse/linalg/matfuncs.py
View file
Open in desktop
Expand Up @@ -844,19 +844,17 @@ def _ell(A, m):
if len(A.shape) != 2 or A.shape[0] != A.shape[1]:
raise ValueError('expected A to be like a square matrix')

p = 2*m + 1

# The c_i are explained in (2.2) and (2.6) of the 2005 expm paper.
# They are coefficients of terms of a generating function series expansion.
choose_2p_p = scipy.special.comb(2*p, p, exact=True)
abs_c_recip = float(choose_2p_p * math.factorial(2*p + 1))
choose_2m_m = scipy.special.comb(2*m, m, exact=True)
abs_c_recip = float(choose_2m_m * math.factorial(2*m + 1))

# This is explained after Eq. (1.2) of the 2009 expm paper.
# It is the "unit roundoff" of IEEE double precision arithmetic.
u = 2**-53

# Compute the one-norm of matrix power p of abs(A).
A_abs_onenorm = _onenorm_matrix_power_nnm(abs(A), p)
A_abs_onenorm = _onenorm_matrix_power_nnm(abs(A), 2*m + 1)

# Treat zero norm as a special case.
if not A_abs_onenorm:
Expand Down