17
\$\begingroup\$

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 0
0 0
1.0 0.0
0.0 1.0
1 2
3 4
51.968956198705044 74.73656456700328
112.10484685050491 164.07380304920997
1 0
0 1
2.718281828459045 0.0
0.0 2.718281828459045
-10 -7
6 3
-0.17051293798604472 -0.22030000635390898
0.18882857687477908 0.23861564524264348
-2 16
0 7
0.1353352832366127 1949.3294633692833
0.0 1096.6331584284585
12 18 -5
-13 13 17
13 -6 2
951375.2972757841 1955306.8594829023 2179360.8077694285
376625.60116007976 774976.2125979062 863826.1366984685
773311.8986313189 1589134.8925863737 1771827.68268726
8 19 20 19
19 -18 8 -11
7 -16 17 -11
13 -15 -14 2
-809927951.1659397 682837927.821331 -2875315029.426385 166307199.77734298
-114105964.84866604 423320553.28643256 -1532090815.2105286 379540651.37782615
666012827.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 16
14 -4 -6 -8 -4 5
-18 2 -14 5 1 -11
-16 18 19 -10 -17 13
10 20 7 19 14 0
-84820410929.4261 -16367909783.470901 -68656483749.58916 3885773007.51203 -53912756108.37766 -68894112255.13809
190335662933.039 -38645269722.440834 127312405236.2376 -13585633716.898304 90603945063.00284 75004079839.71536
-68036952943.18438 -7733451697.302282 -53156358259.70866 3465229815.7224665 -41070570134.5761 -49564275538.3475
60712557398.76749 30529410698.827442 55820038060.925934 -1566782789.1900578 46171305388.15615 69179468777.9944
123964494616.41298 -39882807512.560074 77695806070.41081 -9798106385.28041 53080430956.84853 33312855054.34455
202240615797.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.5581887
28955217908.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.665106
10215509474.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.5699253
449955718.5164527 1851164998.6281173 214390574.08290553 -590752899.2082579 -878315768.622139 552129374.7322844 -332495739.50407004 572740581.3608516
4056736597.835622 16689783857.791903 1932941125.9578402 -5326143353.840331 -7918773134.746702 4977893918.896973 -2997723598.294145 5163693248.841862
18572197375.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
\$\endgroup\$
4
  • \$\begingroup\$ Sandbox \$\endgroup\$
    – pxeger
    Jun 30, 2021 at 15:27
  • \$\begingroup\$ Comments are not for extended discussion; this conversation has been moved to chat. \$\endgroup\$
    – hyper-neutrino
    Jun 30, 2021 at 16:13
  • 5
    \$\begingroup\$ An observation I don't see any answers using so far -- \$e^M \approx (I+M/n)^n\$ for large \$n\$. \$\endgroup\$
    – xnor
    Jul 1, 2021 at 6:56
  • 1
    \$\begingroup\$ @xnor fixed & thanks. \$\endgroup\$ Jul 1, 2021 at 9:17

10 Answers 10

6
\$\begingroup\$

Trivial built-in answers

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!

\$\endgroup\$
4
\$\begingroup\$

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.

\$\endgroup\$
8
  • \$\begingroup\$ Do I count four ¤s? \$\endgroup\$ Jun 30, 2021 at 15:40
  • 1
    \$\begingroup\$ @UnrelatedString yeah this is a really silly answer lol \$\endgroup\$
    – hyper-neutrino
    Jun 30, 2021 at 15:41
  • 1
    \$\begingroup\$ @hyper-neutrino only way it could be sillier would be {compressed string for the sympy builtin to solve the challenge}ŒV \$\endgroup\$
    – pxeger
    Jun 30, 2021 at 15:42
  • \$\begingroup\$ @pxeger done :D \$\endgroup\$
    – hyper-neutrino
    Jun 30, 2021 at 15:47
  • 1
    \$\begingroup\$ @UnrelatedString nice lol \$\endgroup\$
    – hyper-neutrino
    Jun 30, 2021 at 15:59
3
\$\begingroup\$

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
\$\endgroup\$
2
  • \$\begingroup\$ 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) \$\endgroup\$
    – hyper-neutrino
    Jun 30, 2021 at 16:17
  • 1
    \$\begingroup\$ @hyper-neutrino Well, ȷ9 errors, so ȷ1000 would be awful :P \$\endgroup\$ Jun 30, 2021 at 16:18
3
\$\begingroup\$

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!

\$\endgroup\$
3
\$\begingroup\$

R, 85 81 63 61 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...).

\$\endgroup\$
2
  • \$\begingroup\$ 61 bytes? \$\endgroup\$
    – ovs
    Jul 1, 2021 at 9:27
  • \$\begingroup\$ @ovs - Thankyou. Obviously I struggle with 3+3<1+1+6... \$\endgroup\$ Jul 1, 2021 at 9:30
3
\$\begingroup\$

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!

\$\endgroup\$
3
\$\begingroup\$

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!

\$\endgroup\$
1
  • 1
    \$\begingroup\$ -1 infixing #~MatrixPower~n \$\endgroup\$
    – att
    May 13 at 18:41
2
\$\begingroup\$

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.

\$\endgroup\$
1
\$\begingroup\$

Charcoal, 49 bytes

≔EθEι⁼κμη≔ηζFφ«≔Eη∕EκΣEκ×θπν⊕ιηUMζEκ⁺짧ηλν»Iζ

Try it online! Link is to verbose version of code. Explanation:

≔EθEι⁼κμη≔ηζ

Create an identity matrix and save it as the last power series term and also the running total.

Fφ«

Loop 1,000 times.

≔Eη∕EκΣEκ×θπν⊕ιη

Multiply the last power series term by the input and divide by the 1-indexed loop index to create the next term.

UMζEκ⁺짧ηλν

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.

»Iζ

Print the final total.

\$\endgroup\$
1
\$\begingroup\$

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.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.