(DM42) Matrix exponential

[Albert Chan]

I might have over-complicated HP71B code with 10^P scalings.

2^P may be good enough. (5^17<1E12, 1/2^17 to 12 decimals is exact)
With below code, P>17 only if MAXAB(X) ≥ 2^7.5 ≈ 181

Gain in 10^P scaling down is offset by more computation to then scale back up.
With 2^P, we just do P squarings, FNF(K) code not needed.

10 DESTROY ALL @ OPTION BASE 1
20 INPUT "N=";N @ DIM X(N,N),T(N,N),R(N,N)
30 MAT INPUT X
40 P=MAX(0,IROUND(LOG2(MAXAB(X)))+10) @ MAT X=(1/2^P)*X
50 MAT R=(1/4)*X @ MAT T=IDN @ MAT R=R+T
60 MAT R=(1/3)*R @ MAT R=R*X @ MAT R=R+T
70 MAT R=(1/2)*R @ MAT R=R*X @ MAT R=R+T @ MAT R=R*X
80 FOR I=1 TO P @ MAT T=R*R @ MAT R=R+R @ MAT R=R+T @ I; @ NEXT I
90 @@ MAT DISP R @ END

>run
N=3
X(1,1)? 1
X(1,2)? 2
X(1,3)? 3
X(2,1)? 4
X(2,2)? 5
X(2,3)? 6
X(3,1)? 7
X(3,2)? 8
X(3,3)? 9

1  2  3  4  5  6  7  8  9  10  11  12  13

1118905.69945      1374815.06297      1630724.42651
2533881.04197      3113414.03144      3692947.02095
3948856.38452      4852012.99997      5755169.61548

(08-16-2023 04:38 PM)Gjermund Skailand Wrote:  there is a nice article
"Nineteen dubious ways to calcate the exponential of a matrix, tventyfive years later."

A nasty example from the article, less squaring computations might win over perfect scaling.

XCas> approx(EXPM1([[-49,24], [-64, 31]])) /* reference */

-1.73575875814      0.551819099658
-1.47151759909      0.103638240716

HP71B, SCALE BY 10^5 (≈ 2^(5*3.322) = 2^16.6)

-1.73575875912      0.551819100440
-1.47151760119      0.103638242380

HP71B, SCALE BY 2^16 or 2^17

-1.73575875814      0.551819099660
-1.47151759909      0.103638240710