# To raise $e$ to the power of a matrix

Inspired by this 3Blue1Brown video

Given a square matrix $$\ M \$$, compute its matrix exponential $$\ \exp(M) \$$, which is defined, using an extension of the Maclaurin series for $$\ e^x \$$, as

$$\exp(M) = \sum_{r=0}^\infty \frac {M^r} {r!} = M^0 + M^1 + \frac 1 2 M^2 + \frac 1 6 M^3 + \cdots + \frac 1 {n!} M^n + \cdots$$

where $$\ n! \$$ represents the factorial of $$\ n \$$, and $$\ M^0 \$$ is the identity matrix for the dimensions of $$\ M \$$.

There are other ways to compute this, which you may use, as long as the result is sufficiently precise (see the rules below).

## Test cases

Input Output
0 00 0
1.0 0.00.0 1.0
1 23 4
51.968956198705044 74.73656456700328112.10484685050491 164.07380304920997
1 00 1
2.718281828459045 0.00.0 2.718281828459045
-10 -76 3
-0.17051293798604472 -0.220300006353908980.18882857687477908 0.23861564524264348
-2 160 7
0.1353352832366127 1949.32946336928330.0 1096.6331584284585
12 18 -5-13 13 1713 -6 2
951375.2972757841 1955306.8594829023 2179360.8077694285376625.60116007976 774976.2125979062 863826.1366984685773311.8986313189 1589134.8925863737 1771827.68268726
8 19 20 1919 -18 8 -117 -16 17 -1113 -15 -14 2
-809927951.1659397 682837927.821331 -2875315029.426385 166307199.77734298-114105964.84866604 423320553.28643256 -1532090815.2105286 379540651.37782615666012827.4455533 256623519.77362177 -454247177.983024 753881172.0779059-849659694.5821244 -147765559.4347415 -39332769.14778117 -766216945.8367432
15 -17 7 -1 18 5-20 1 -11 -11 -2 1614 -4 -6 -8 -4 5-18 2 -14 5 1 -11-16 18 19 -10 -17 1310 20 7 19 14 0
-84820410929.4261 -16367909783.470901 -68656483749.58916 3885773007.51203 -53912756108.37766 -68894112255.13809190335662933.039 -38645269722.440834 127312405236.2376 -13585633716.898304 90603945063.00284 75004079839.71536-68036952943.18438 -7733451697.302282 -53156358259.70866 3465229815.7224665 -41070570134.5761 -49564275538.347560712557398.76749 30529410698.827442 55820038060.925934 -1566782789.1900578 46171305388.15615 69179468777.9944123964494616.41298 -39882807512.560074 77695806070.41081 -9798106385.28041 53080430956.84853 33312855054.34455202240615797.98032 -49846425749.36303 132157848306.15779 -15002452609.223932 92731071983.4513 70419737049.6608
-3 3 9 -14 13 3 -19 11-3 16 -3 -2 -16 17 -7 14-16 -13 -19 -4 -19 -12 -19 4-19 2 -1 -13 -1 20 -18 20-15 -14 -17 4 -16 -7 -13 10-1 3 -2 -18 -13 -20 -18 8-6 5 17 4 -11 0 4 1-7 14 4 5 -10 1 11 -1
-961464430.42625 -3955535120.8927402 -458113493.1060377 1262316775.4449253 1876774239.173575 -1179776408.054209 710474104.2845823 -1223811014.558188728955217908.989292 119124631307.93314 13796523822.599554 -38015726498.96707 -56520887984.67961 35530121226.97329 -21396437283.72946 36856280546.42262-8410889774.023839 -34603239307.789085 -4007607155.9532456 11042781096.475042 16418151308.196218 -10320764772.97249 6215219812.505076 -10705984738.66510610215509474.424953 42027619363.9107 4867469315.8131275 -13412092189.39047 -19940786719.11994 12535160455.72014 -7548741937.235227 13003031639.209038-1859396787.0195892 -7649733581.4828005 -885954562.2162387 2441226246.193038 3629550445.402215 -2281610372.751828 1374002295.125188 -2366775855.5699253449955718.5164527 1851164998.6281173 214390574.08290553 -590752899.2082579 -878315768.622139 552129374.7322844 -332495739.50407004 572740581.36085164056736597.835622 16689783857.791903 1932941125.9578402 -5326143353.840331 -7918773134.746702 4977893918.896973 -2997723598.294145 5163693248.84186218572197375.577248 76407841992.77576 8849246673.162008 -24383706828.81331 -36253120255.06763 22789406560.399803 -13723910211.58447 23640014943.24763

## Rules

• Your outputs must be within $$\ \pm 1 \% \$$ of the outputs given in the test cases above
• If you use a boring builtin to do most of the computation, you should add it to the Community Wiki answer, or post another more interesting method as well
• You may assume $$\ M \$$ will be square and have side length in $$\ [2, 8] \$$
• The elements of $$\ M \$$ will all be integers in $$\ [-20, 20] \$$
• You may take $$\ M \$$ as a nested array, a built-in matrix type, a flat array which is a square-number in length, or any other sensible format
• Standard loopholes are forbidden
• Standard I/O rules apply
• This is , so the shortest code in bytes wins
• Sandbox Jun 30, 2021 at 15:27
• Comments are not for extended discussion; this conversation has been moved to chat. Jun 30, 2021 at 16:13
• An observation I don't see any answers using so far -- $e^M \approx (I+M/n)^n$ for large $n$.
– xnor
Jul 1, 2021 at 6:56
• @xnor fixed & thanks. Jul 1, 2021 at 9:17

If your solution uses a built-in function to perform the core of the challenge, please edit it in to this answer.

# Wolfram Language (Mathematica), 9 bytes

MatrixExp


# MATL, 2 bytes

Ye


Try it online!

# MATLAB/Octave, 5 bytes

@expm


Try it online!

# Python 3, 30 bytes

from scipy.linalg import*
expm


Try it online!

# R, 12 bytes

Matrix::expm


Try it online!

# Julia, 3 bytes

exp


Try it online!

# Jelly, 18 bytes

-©‘+³æ*‘ɼ¤÷®!¤¤$ÐL  Try it online! -5 bytes thanks to caird coinheringaahing This is probably the least Jelly-like Jelly answer ever. Goes for unbounded precision, even though 1/256! (which is what caird coinheringaahing's solution does) is definitely enough in practice. This only terminates due to limited precision. If we had infinite/unbounded precision, this would never stop running. It basically just keeps adding terms until the result stops changing, which means we've hit the precision limit. Or actually, this is more cursed: # Jelly, 34 bytes ŒṘ“çm(<ÇƤMʂƲṬVḅċñ!żŒȥu¢iİẠẈ“œ4»jŒV  Try it online! -3 bytes thanks to Unrelated String. • Do I count four ¤s? Jun 30, 2021 at 15:40 • @UnrelatedString yeah this is a really silly answer lol Jun 30, 2021 at 15:41 • @hyper-neutrino only way it could be sillier would be {compressed string for the sympy builtin to solve the challenge}ŒV Jun 30, 2021 at 15:42 • @pxeger done :D Jun 30, 2021 at 15:47 • @UnrelatedString nice lol Jun 30, 2021 at 15:59 # Jelly, 11 bytes æ*÷!}ɗⱮ⁹Ḷ¤S  Try it online! ## How it works Basically, we calculate $$\sum_{r=0}^{255} \frac {M^r} {r!}$$ which is large enough to not have to worry about precision. For one additional byte, we can instead calculate $$\sum_{r=0}^{10^9} \frac {M^r} {r!}$$, which is, for all intents and purposes, equal to $$\sum_{r=0}^\infty \frac {M^r} {r!}$$ æ*÷!}ɗⱮ⁹Ḷ¤S - Main link. Takes M on the left ¤ - Group the previous links into a constant: ⁹ - 256 Ḷ - Lowered range; [0, 1, ..., 255] ɗ - Group the previous 3 links into a dyad f(M, i): æ* - Raise M to the i'th power } - To i: ! - Factorial ÷ - Divide each element of Mⁱ by i! Ɱ - Over each i in [0, 1, ..., 255], calculate f(M, i) S - Sum  • Or, for a few extra bytes, we can do ȷ1000 which, at least according to Javascript, literally is infinity :P (and also will never terminate) Jun 30, 2021 at 16:17 • @hyper-neutrino Well, ȷ9 errors, so ȷ1000 would be awful :P Jun 30, 2021 at 16:18 # JavaScript (ES6), 144 bytes This iterates until $$\1/n!\$$ becomes $$\0\$$ because of IEEE 754 underflow. f=(m,M=m,S=m.map(r=>r.map(_=>j=0)),q=1)=>q?f(m,M.map((r,y)=>r.map((_,x)=>(S[y][x]+=(v=j?r.reduce((t,v,i)=>t+v*m[i][x],0):x==y)*q,v))),S,q/++j):S  Try it online! # R, 85816361 59 bytes Edit: -18 bytes & new approach thanks to xnor's comment, and then -2 more bytes thanks to ovs function(m){k=diag(nrow(m))+m/2^25;for(i in 1:25)k=k%*%k;k}  Try it online! Non-builtin solution (see here to save 69 47 bytes...). • – ovs Jul 1, 2021 at 9:27 • @ovs - Thankyou. Obviously I struggle with 3+3<1+1+6... Jul 1, 2021 at 9:30 # Pari/GP, 16 bytes Based on xnor's comment: An observation I don't see any answers using so far -- $$\e^M\approx(I+M/n)^n\$$ for large $$\n\$$. It seems that $$\n = 9!\$$ is enough for the $$\ \pm1 \%\$$ requirement of the testcases. Sorry. It turns out that $$\n = 9!\$$ is not good enough. It exceed the $$\ \pm1 \%\$$ requirement in exactly one value: the -39332769.14778117 in the 7th testcase, where $$\n = 9!\$$ gives -38882430.43811763. So I have to use $$\n = 99!\$$. m->(1+m/99!)^99!  Try it online! # Wolfram Language (Mathematica), 54 53 bytes I'm still not excellent at golfing, so additions are very welcome. Originally I tried to use only Sum[MatrixPower[#,n],{n,0,∞}]&, but that doesn't work for a matrix of only zeros as MatrixPower doesn't work on singular matrices. IdentityMatrix@Tr[1^#]+Sum[#~MatrixPower~n/n!,{n,∞}]&  It can respond with symbolic results for all inputs, so {{1,0},{0,1}} gives an answer of {{E,0},{0,E}}, and {{-2, 16}, {0, 7}} returns {{1+(1-E)^2/E^2, 16(E^9-1)/(9E^2)}, {0, E^7}}. If floating point answers are wanted N@ can be added in front of IdentityMatrix. -1 byte, thanks @att. Try it online! • -1 infixing #~MatrixPower~n – att May 13 at 18:41 # J, 31 bytes [:+/!@i.@99%~]+/ .*^:(<99)[:=#\  Try it online! Just implements the formula up to 99 terms, which seems to be good enough to be within 1% of the test cases. # Charcoal, 49 bytes ≔ＥθＥι⁼κμη≔ηζＦφ«≔Ｅη∕ＥκΣＥκ×ξ§§θπν⊕ιηＵＭζＥκ⁺μ§§ηλν»Ｉζ  Try it online! Link is to verbose version of code. Explanation: ≔ＥθＥι⁼κμη≔ηζ  Create an identity matrix and save it as the last power series term and also the running total. Ｆφ«  Loop 1,000 times. ≔Ｅη∕ＥκΣＥκ×ξ§§θπν⊕ιη  Multiply the last power series term by the input and divide by the 1-indexed loop index to create the next term. ＵＭζＥκ⁺μ§§ηλν  Add the term to the running total in-place, as the term is now a separate array and no longer shared with the running total. »Ｉζ  Print the final total. # Raku, 118 113 bytes {[»+«] (flat((1,0 xx$_)xx$_).rotor($_)[^$_],{cross(with=>(*Z* *).sum,@^x »/»++$,[Z] $_).rotor($_)}...*)[^99]}


Try it online!

• [»+«] reduces the following sequence of 99 matrices (each a list-of-lists) with elementwise addition. (99 seems like more than enough terms for any practical purpose.)
• flat((1, 0 xx $_) xx$_).rotor($_)[^$_] is the starting identity matrix.

I found a somewhat concise way to do matrix multiplication on lists-of-lists @m1 and @m2 in Raku:

cross(with => (* Z* *).sum, @m1, [Z] @m2).rotor(@m1)


cross here is the cross product operation. (* Z* *).sum is an anonymous function applied to each combination of rows from @m1 and columns from @m2, zipping them with multiplication and then summing. [Z] @m2 is a tricky way to get the transpose of the right-hand matrix @m2. That gives just a flat list of numbers, which .rotor(@m1) turns back into a square matrix.

@^x is the previous matrix in the series, which I divide elementwise by a steadily increasing anonymous state variable with »/» ++$ before multiplying by the input matrix $_ to get the same effect as dividing by the factorial in the original formula.