Why expm(2*A) != expm(A) @ expm(A)

Question

According to the rule exp(A+B) = exp(A)exp(B), which holds for commuting matrices A and B, i.e. when AB = BA, we have that exp(2A) = exp(A)exp(A). However when I run the following in Python:

import numpy as np
from scipy.linalg import expm

A = np.arange(1,17).reshape(4,4)

print(expm(2*A))
[[ 306.63168024  344.81465009  380.01335176  432.47730444]
 [ 172.59336774  195.36562731  214.19453937  243.76985501]
 [ -35.40485583  -39.87705598  -42.94545895  -50.01324379]
 [-168.44316833 -190.32607875 -209.76427134 -237.72069322]]

print(expm(A) @ expm(A))
[[1.87271814e+30 2.12068332e+30 2.36864850e+30 2.61661368e+30]
 [4.32685652e+30 4.89977229e+30 5.47268806e+30 6.04560383e+30]
 [6.78099490e+30 7.67886126e+30 8.57672762e+30 9.47459398e+30]
 [9.23513328e+30 1.04579502e+31 1.16807672e+31 1.29035841e+31]]

I get two very different results. Note that @ is just the matrix product.

I also tried it in Matlab and the two results are the same as expected. What am I missing here?

Edit: I have NumPy 1.15.3, SciPy 1.1.0, Python 3.6.4, Windows 7 64-bit

As suggested in comments by Warren Weckesser, using A = A.astype(np.float64) solves the problem.

Answer 1

In brief: there's a bug in scipy 1.1.0 which seems to have been fixed in 1.2.0.

Installing the latest scipy (1.2.1), the following passes just fine:

import numpy as np
from scipy.linalg import expm

A = np.arange(1,17).reshape(4,4)

assert (expm(A) @ expm(A) == expm(2 * A)).all()