(DM42) Matrix exponential

[Albert Chan]

(08-16-2023 09:48 PM)Albert Chan Wrote:  Perhaps we could do z = expm1(x/2^n), so that overflow is not an issue.

On decimal machine, we like z = expm1(x/10^n), to remove scaling errors.

\( (e^x-1)(e^y-1) = e^x e^y - e^x - e^y + 1 = (e^x e^y-1) - (e^x-1) - (e^y-1) \)

expm1(x+y) = expm1(x) * expm1(y) + expm1(x) + expm1(y)

Using this identity, below code does expm1(log1p(r)*n) = (1+r)^n - 1

This is an important term used in TVM calculations.
(Formula arranged to produce Plus42 menus order: N, I, PV, PMT, FV, BEGIN)

TVM:0=(1/EXPM1(N*LNP1(I))*(PV+0*PMT+FV)+PV+BEGIN*PMT)*I+PMT

Code:

function f(r, n)                -- expm1(log1p(r)*n), integer n>=0
    if n<=1 then return r*n end
    local R = f(r, floor(n/2))
    R = R * (R+2)               -- expm1(log1p(r)*(n-n%2))
    return n%2==0 and R or r*R + r + R
end

lua> r = 0.01
lua> f(r,0), f(r,1), f(r,2), f(r,3)
0      0.01      0.0201      0.030301

expm1(log1p(r)) = r, we can build from expm1(r) to expm1(r * n)

lua> n = 100
lua> f(expm1(r), n)
1.7182818284590453
lua> expm1(r * n)
1.7182818284590453