(DM42) Matrix exponential

[Albert Chan]

EXPM1(Matrix X) for HP71B

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 S=10^MAX(0,EXPONENT(MAXAB(X))+4) @ MAT X=(1/S)*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
80 MAT R=R*X @ MAT X=R
90 DISP FNF(S) @ DISP @ MAT DISP R @ END

200 DEF FNF(K) @ IF K<=1 THEN 250
210 DISP FNF(K DIV 2);
220 MAT T=R*R @ MAT R=R+R @ MAT R=R+T
230 IF MOD(K,2)=0 THEN 250
240 MAT T=R*X @ MAT T=T+X @ MAT R=R+T
250 FNF=K @ END DEF

>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  4  9  19  39  78  156  312  625  1250  2500  5000  10000

1118905.69939      1374815.06290      1630724.42641
2533881.04184      3113414.03128      3692947.02074
3948856.38427      4852012.99966      5755169.61507