TVM solve for interest rate, revisited

[Albert Chan]

f = ((pv+fv)/(K-1) + pv)*x + pmt      , where K = (1+x)^n

To improve accuracy, we noticed ULP(f) ≥ ULP(pmt)
If pv=0, we have ULP(f) = ULP(pmt)

We like to make |pv| smaller, if possible, to increase f precision.
We can do this with time-symmetry (formula does not care time direction):

NPV = NFV / K  = NPMT / C

K = (1+x)^n
C = (x*n)/(1-1/K) = (x*n)/(K-1) * K = (x*n)/(K-1) + n*x

NPV(n,i,pv,pmt,fv) = NFV(-n,i,fv,-pmt,pv)
NFV(n,i,pv,pmt,fv) = NPV(-n,i,fv,-pmt,pv)
NPMT(n,i,pv,pmt,fv) = NPMT(-n,i,fv,-pmt,pv)

f = NPMT/n      → f(n,i,pv,pmt,fv) = -f(-n,i,fv,-pmt,pv) = f(-n,i,-fv,pmt,-pv)

Code:

function _find_rate(iter_i, found, n,pv,pmt,fv,i0, verbose)
    if abs(pv) > abs(fv) then n,pv,fv = -n,-fv,-pv end
    local c, f0 = 2
    for i,eps,f in iter_i(n,pv,pmt,fv,i0) do
        if verbose then print(i,f) end
        if c>0 then c=c-1; f0=f; continue end
        if found(f,f0) then return i-eps/2, true end
        if abs(f) < abs(f0) then f0=f; continue end
        return i-eps/2, (i == i-eps*0.001)
    end
end

For convenience, I made function fn, that act like Python's lambda.
Note: Halley's method for f=0, convergence may not be 1-sided.

Code:

function find_rate(...) return _find_rate(iter_i, fn'a,b: a*b<=0', ...) end
function find_rate2(...) return _find_rate(iter_i2, fn'a,b: a==0', ...) end

Update, Jul 11,2022:

Rate iterators actually returned (i+eps), eps, f
In other words, (i,f) is the point, (i+eps) extrapolated value.

Final f of only a few ULP's, (i+eps) may not be that good.
A reasonable guess is true rate between i and (i+eps).
Patched code returned (i+eps/2) = (i+eps) - eps/2
This is patched result:

lua> n,pv,pmt,fv = 10, -100, 10, 1e10
lua> find_rate2(n,pv,pmt,fv,nil,true)
0.3749999875677088      999999995.5
0.7951945637769868      161944286.30682382
1.3009220024991999      22940767.56093494
1.9230635724160559      3128383.8915700433
2.694052015186147       422117.19292708224
3.6428250853054216      56675.26724831162
4.70010922758184        7473.075259553417
5.291135425226106       837.9935225675524
5.32155201095865        25.713350899574266
5.321554259156788       0.0018612216055089448
5.321554259156787       -1.2505552149377763e-012
5.321554259156789       2.3874235921539366e-012
5.3215542591567875      true

Above verbose output, f value should read as if whole column get shifted up.

(i1, f1) = (0.375, 162.e6)
(i2, f2) = (0.795, 22.9e6)
...

It showed a sign flip when i ends in 788, flipped again when i ends in 787
With f down to few ULP's, correction is worse, with final i ends in 789, outside bracketed rate.
Patched code assumed true i ends in middle of (787, 789)

Update, Jul 18,2022:

fuzz factor for rate found widened, with less false negatives.

< return i-eps/2, (i == i-eps*0.01)
> return i-eps/2, (i == i-eps*0.001)