(DM42) Matrix exponential

[Gjermund Skailand]

updated 12/8

A matrix exponential function for square matrix for the DM42, Free42

This is just a simple implementation based argument reduction and series expansion.
Exp(M) = I + M + 1/(2!)M^2 + 1/(3!)M^3... 1(n!)M^n
I is identity matrix

The argument reduction reduces the fnrm to below 1.0 prior to the series expansion.

Note that running times will increase for arrays with large inputs.
Program runs in 4 stack mode.
Flag 24 and 25 shall be cleared.
Uses registers R00-R02
R00: counter number of convergence iterations
R01: stopping criteria is difference between row norm for two iterations to be less than value set in line 30. Test examples use 1e-32, but a more sensible value is probably 1e-16.
R02: scratch, used for argument reduction and squaring

Code:


@ Matrix exponential by series expansion. 
@ scaling by argument reduction  
@ uses R0-02, flag 24,25 cleared
00 { 104-Byte Prgm }
01▸LBL "MEXP"
02 FUNC 11
03 L4STK
04 2
05 STO 00
06 RCL ST Y
07 FNRM
08 LN
09 RCL 00
10 LN
11 ÷
12 IP
13 1
14 STO 01
15 +
16 STO 02
17 Y↑X
18 ÷
19 ENTER
20 ENTER
21 EDIT
22 ↑
23▸LBL 01
24 ↓
25 RCL+ 01
26 →
27 FC? 77
28 GTO 01
29 EXITALL
30 1ᴇ-32   @ 1ᴇ-16 convergence test 
31 STO 01
32▸LBL 02
33 R↓
34 X<>Y
35 RCL÷ 00
36 ISG 00
37 NOP
38 RCL× ST Z
39 X<>Y
40 RCL+ ST Y
41 STO- ST L
42 LASTX
43 RNRM
44 RCL- 01
45 X>0?
46 GTO 02
47 R↓
48 0=? 02
49 GTO 04
50▸LBL 03
51 ENTER
52 ×
53 DSE 02
54 GTO 03
55▸LBL 04
56 FNRM
57 LASTX
58 END

example 1
M =
[[ 1 4 ]
[ 1 1 ]]
MEXP=
[[ 10.226708, 19.717657]
[ 4.929414, 10.226708]]

exact (Matrix exponential - Wikipedia):
[[(exp(4)+1)/(2e), (exp(4)-1)/e ]
[(exp(4)-1)/(4e), (exp(4)+1)/(2e)]]

difference between MEXP(M)- exact =
[[ 3.0E-32, 6.0E-32]
[ 1.9E-32, 3.0E-32 ]]

example 2, calculated in about 1.6 s when running on batteries
M =
[[1 2 3 ]
[4 5 6]
[7 8 9]]
MEXP(M) =
[[ 1.118907e6, 1.374815e6, 1.630724e6]
[ 2.533881e6, 3.113415e6, 3.692947e6]
[ 3.948856e6, 4.852013e6, 5.755171e6]]

ps It might be more sensible to use R1=1E-16, as it cuts down calculation time almost in half, and if in doubt about the result, do a second run with R1=1E-32 to compare the first 16 digits.

Best regards
Gjermund