(DM42) Matrix exponential

[Albert Chan]

(08-23-2023 09:39 PM)Albert Chan Wrote:  exmp1(k*Jn) = expm1(k*n) * (Jn/n)

My old code epsilon has enough safety margin for small matrices.
But, in general, perhaps we need to cover the edge case.

If n×n elements sized about the same, say k, FNORM(X) ≈ k*n

Code:

< 40 P=MAX(0,IROUND(LOG2(MAXAB(X)))+10) @ MAT X=(1/2^P)*X

> 40 P=MAX(0,IROUND(LOG2(FNORM(X)))+9) @ MAT X=(1/2^P)*X

Example, expm1(J20) = expm1(20)/20 * J20 = 24258259.7205 * J20

Old Patch: P=10 ---> 24258259.1465 * J20
New Patch: P=13 --> 24258259.7183 * J20