Looking for TVM contributions

[robve]

(06-13-2024 09:28 PM)Albert Chan Wrote:  

(06-12-2024 08:48 PM)robve Wrote:  | 11 | robve | 60x24x365 | 1/6% → ? | 0 | -0.01 | ?→ FV | =N | yes |
Sharp PC-1421 FV=5258.756614 => I%=1.073741824E+09 !!!

(06-12-2024 11:56 PM)dm319 Wrote:  Very interesting! Weird how they can be similar and so different, and what is going on with 11?!

Initial guess (hard coded?) may be too high for the problem.
[...]
What if we supply a guess, say 0.1% ?
[...]
Why does it fail? Because we use Newton's method from the wrong edge!
The other edge = pmt/-pv = ∞, but 1e-3 is just as bad. (technically, ∞ is not an edge)

Great explanation with details and plots, Albert. This machine may look very bad from this result, but it is a weakness of root finding methods in practice when implemented and it gets worse with limited numerical precision on vintage calculators. I suspected that the slope would be close to zero for it to jump that far from the root. Corrections to Newton can be made, see e.g. R.W. Hamming Numerical Methods for Scientists and Engineers "modification is that if |f(z_new)| > |f(z)| then halve step, adjust z_new and try again". But adding more logic to implement TVM may not have been possible for the PC-1421 which has only 40K ROM and is otherwise the same as the PC-1401. It is remarkable SHARP managed to squeeze all common financial functions into a tiny bit of remaining bytes in unused PC-1401 ROM.

Based on comments on this and other HP forum threads, I was curious to test Newton versus Secant for TVM root finding, especially on SHARP and similar machines that have only 10 to 12 digits BCD to work with.

After some experimentation in C code with double precision and in SHARP BASIC, I am not sure if others noticed this as well, but I noticed that Newton suffers from instability when TVM rates converge to zero. It is particularly bad with 10 to 12 BCD calculators like the SHARP where the root finder does not get closer than 1E-5 to 1E-7 when the root I%=0 or very close to it. Perhaps Newton's convergence rate is too high, which could explain why minor rounding errors cause significant overshooting?

Lines 31, 33 and 34 are replacing the Secant method lines 31, 33 to 35 and the rest of the program is the same as in the Secant version. The parts updated are (added context for completeness):

Code:

' check if i%=0 (eval order of P+F+N*M is important)
29 C=0 : IF P+F+N*M=0 LET I=0 : RETURN
' solve Y=NPMT=0 for rate I% with the newton method, pick initial nonzero guess
30 Y=(P+F)/N,W=F-P-Y,V=(M+Y)/W,W=(N*N-1)*Y*V/W,I=100*V*(W-3)/(W-1.5)
31 IF I=0 LET I=1E-5
' newton loop to find root of Y=NPMT
33 GOSUB 40 : H=(K*M-(P+F*S)*J/R)*100/(B*M-(P+F)*(1+N*S/(R+R/J))/R-F)
34 I=I-H : If I+H*H>I GOTO 33
' round very small I% to zero
36 I=(ABS I>1E-9)*I : RETURN
' check edge case I%=0 to set J,K,S,R and compute N when requested
39 C=0 : IF I=0 LET J=1,K=1,S=1,R=-N,N=N-(N=0)*(P+F)/(M+(M=0)) : RETURN
' compute L=LN(1+J) and S=(1+J)^-N and R=S-1 using LNP1 and EXPM1 formulas
40 J=.01*I,K=1+J*B,L=1+J,L=LN L-(L-1-J)/L,S=-N*L,R=EXP S,R=R-1-(LN R-S)*R,S=R+1 : RETURN

I wrote the long expression for H to make good use SHARP's 12 BCD internal precision. I also tested the slightly longer formulation from Albert's previous posts on Newton with an HP-71B example:

Code:

H=(-(P+F)*J/R+(B*M-F)*J+M)*100/(-(P+F)/R*(1+N*S/(R+R/J))+(B*M-F))

Note that there are some minor differences to Albert's code: R is S but not negated and S=R+1, K=1+J*B with J the rate (I%/100) see line 40. This way I can reuse K, L, R and S to compute N, PV, PMT and FV as before. Note that the convergence test I+H*H>I is also used, but it has problems because H*H assumes ideal convergence by the next step without rounding issues. In reality H may not be quadratically converging to zero. This can be observed in the following cases using 10 to 12 digits BCD arithmetic (I suspect IEEE754 binary may suffer similar issues though much more subtly and less often):

4) n=480 pv=100000 pmt=-208.333 fv=0 does not terminate because I% gets stuck around 1E-07 and H around 1E-08

5) n=10 pv=50 pmt=-30 fv=400 is slightly worse than before with 14.43587131 and does not converge to 14.435871133 after forced Newton iterations

7) n=10 pv=-100 pmt=10 fv=1e-10 does not converge after removing the P+F+N*M=0 check for zero rates (a simple "cheat") because I% gets stuck around 1E-05

12) n=40 pv=900 pmt=-400 fv=-1000 BGN mode gives 80.00000003% instead of 80% or 80.00000001%. Additional iterations will flip/flop around 80.00000001%. It is a minor error, but an interesting difference nonetheless.

I suspect that the convergence rate of Newton is too high, which causes rounding errors and cancellations to accumulate and cause the method to behave more erratically for small rates to solve.

Another practical downside of Newton on these vintage machines is that the step computations are more costly to perform than Secant. This negates the benefit of faster theoretical quadratic convergence with Newton versus Secant that converges theoretically with the power of the golden ratio 1.618... which is still very good. Secant also appears to correct itself when it overshoots the root and in the presence of rounding errors. But that's still anecdotal given the smallish set of test cases (about 25 that I use, including the 12 of this thread).

Furthermore, an initial good rate guess (almost) eliminates the advantage of using Newton for speed compared to Secant. So perhaps Secant is the better method? At least for 10 to 14 digit BCD calculators and pocket computers? I don't see how Newton can be improved to improve over the Secant results. It's a bit disappointing because the formulas posted here and here are clever and elegant.

- Rob

[Albert Chan]

(06-14-2024 03:23 PM)robve Wrote:  I suspected that the slope would be close to zero for it to jump that far from the root.

NPV iterating to infinity is not a rounding issue ... go figure!

NPV = (fv-pmt/x) * ((1+x)^-n-1) + (pv+fv)

NPV(n>1, x→∞) ≈ (fv-0) * (0-1) + fv = -fv + (pv+fv) = pv      // =0, for previous post example

Update: it may be better to use another NPV version.
Time-reversal to get NFV, then scale it back down for NPV

NFV = (pv+pmt/x) * ((1+x)^n-1) + (pv+fv)
NPV = NFV / (1+x)^n = (pv+pmt/x) * (1-(1+x)^-n) + (pv+fv)/(1+x)^n

NPV(n>1, x→∞) = (pv+pmt/x) --> pv

Quote:After some experimentation in C code with double precision and in SHARP BASIC, I am not sure if others noticed this as well, but I noticed that Newton suffers from instability when TVM rates converge to zero.

Secant's method have this issue too. Both methods need accurate f(ε) calculation.
It is not noticed because solved rate is not small enough, not hitting f(ε) issue yet.
Also, with 2 very close points, f'(ε) estimated from secant slope is likely worse.

After I use {f(0), f'(0), f''(0)} to estimate {f(ε), f'(ε)}, this issue mostly goes away.

BTW: small rate branch is primarily designed in case edge rate has wrong sign.
Next iteration(s) may pass through f(ε), before it reaches true rate on the other side.
Getting true tiny rate solution is really a bonus.

If you like to try Newton's method, use this instead (convergence test using f(i), not i)
Convergence depends on ULP, not quadratic convergence, thus will never get stuck.
https://www.hpmuseum.org/forum/thread-16...#pid187993

Note that K=1+B*I is removed. It just seems to introduce unnecessary error.
Assuming (pv+pmt) and (fv-pmt) is exact, this may be better. (B=1 → B=0)

Code:

function tvm_begin(n,i,pv,pmt,fv, verbose)
    if not fv  then return tvm(n,i,pv+pmt,pmt,fv)+pmt end
    if not pv  then return tvm(n,i,pv,pmt,fv-pmt)-pmt end
    if not pmt then return tvm(n,i,pv,pmt,fv) / (1+i) end
    return tvm(n,i,pv+pmt,pmt,fv-pmt, verbose)
end

Quote:Another practical downside of Newton on these vintage machines is that the step computations are more costly to perform than Secant. This negates the benefit of faster theoretical quadratic convergence with Newton versus Secant that converges theoretically with the power of the golden ratio 1.618... which is still very good.

Per function calls, Secant's method may beat Newtons.
However, f(i) and f'(i) shared expensive expm1(log1p(i)*n) factor.
My guess is f'(i) cost very little extra ... say 1/4 of a function call.
Secant's method needs to calculate slope too ... it is not free.

Besides getting f(i) and slope, there are bookkeeping cost too.
Per timing, with less loops, my guess is Newtons beat Secant.

Of course, this is all guessing. We know once the code is done.

[robve]

(06-14-2024 06:40 PM)Albert Chan Wrote:  

(06-14-2024 03:23 PM)robve Wrote:  I suspected that the slope would be close to zero for it to jump that far from the root.

NPV iterating to infinity is not a rounding issue ... go figure!

Right! Initial guesses matter. Can happen with any TVM calc. In the TVM code that I had before I halved I% from an initial guess I%=0.01 until |Y-W| (unscaled "slope" sort of) is larger than the I% rate itself, which prevented this effectively. But that is very ad-hoc.

(06-14-2024 06:40 PM)Albert Chan Wrote:  

Quote:After some experimentation in C code with double precision and in SHARP BASIC, I am not sure if others noticed this as well, but I noticed that Newton suffers from instability when TVM rates converge to zero.

Secant's method have this issue too.

Agree. But for Newton running on this 10 BCD machine the longer f(i)/f'(i) expression seems more likely to introduce floating point inaccuracies than the simpler Secant. Whereas the IEEE754 double precision version I have is more well behaved. It's one possible explanation.

(06-14-2024 06:40 PM)Albert Chan Wrote:  Both methods need accurate f(ε) calculation.
It is not noticed because solved rate is not small enough, not hitting f(ε) issue yet.
Also, with 2 very close points, f'(ε) estimated from secant slope is likely worse.

The Secant slope estimation is very simple and it is not the true slope f'(ε), so surely Secant will dance around the root as well. That's what we know. I'm just seeing that it converges more nicely than Newton for small rates ε for the example cases we have and the convergence test to use is reliable too. A more expansive set of test examples could tell if this is an artifact of these examples or if this is general.

Here is one example. With Newton when pushing on to get to y=0 it hovers for a moment:

Code:

n=480 pv=100000 pmt=-208.333 fv=0
 1/12  i%=             1e-05  y=  0.00501045666647   9.99992016789144771e-06
 2/11  i%= 7.98321085531e-11  y=   3.999917908e-08   7.98320053540206081e-11
 3/10  i%= 1.03199085624e-16  y= 2.84217094304e-14   4.22703454717260951e-17
 4/9  i%= 6.09287401526e-17  y= 2.84217094304e-14   7.48692358995212396e-17
 5/8  i%=-1.39404957469e-17  y=                 0                         0
i%=0%

With Secant this converges right away:

Code:

n=480 pv=100000 pmt=-208.333 fv=0
i%=1e-05    y=0.00501046
 3/12  i%= 3.99161666854e-11  y= 1.99996748051e-08   4.99996008383331496e-06
 4/12  i%= 1.35138613467e-16  y= 2.84217094304e-14   3.99160315468350466e-11
 5/12  i%= 7.84135180081e-17  y=                 0   5.67250954586066735e-17
i%=0%

Now, this is just an example, so don't take my word for it.

What I see is that convergence issues with the Newton version just happen more in 10 BCD precision when there is a loss of 4 to 5 digits. Also when to stop is tricky. Even in IEEE754 double precision, such as for these:

n=40 pv=900 pmt=1000 fv=-1000 BGN
n=360 pv=0 pmt=-1000 fv=1e+06
n=5 pv=-369494 pmt=17500 fv=540000
n=30 pv=-3200 pmt=-717.44 fv=60000
n=60 pv=243400 pmt=-1363.29 fv=-222976
n=6 pv=-32000 pmt=0 fv=28347
n=60 pv=8000 pmt=-150.97 fv=0

(06-14-2024 06:40 PM)Albert Chan Wrote:  After I use {f(0), f'(0), f''(0)} to estimate {f(ε), f'(ε)}, this issue mostly goes away.

Makes sense, but this code isn't using {f(0), f'(0), f''(0)}. Also, won't that increase the cost to execute?

(06-14-2024 06:40 PM)Albert Chan Wrote:  If you like to try Newton's method, use this instead (convergence test using f(i), not i)
Convergence depends on ULP, not quadratic convergence, thus will never get stuck.
https://www.hpmuseum.org/forum/thread-16...#pid187993

The rate gets stuck to about 1E-5 in one case and 1E-7 in another which is too far from the root at zero to be detected as a possible zero. One explanation is the use of 10 to 12 BCD (2 guard digits) with lower precision EXPM1 and LOG1P. So losing as much as 4 to 5 digits which is not great. A reliable convergence test should take this into account.

(06-14-2024 06:40 PM)Albert Chan Wrote:  Per function calls, Secant's method may beat Newtons.
However, f(i) and f'(i) shared expensive expm1(log1p(i)*n) factor.
My guess is f'(i) cost very little extra ... say 1/4 of a function call.
Secant's method needs to calculate slope too ... it is not free.

But it's almost free, right? Just iterate \( i_{\rm new} := i-(i-i_{\rm old})f/(f-f_{\rm old}) \) with one evaluation of npmt to get f(i) to update i.

(06-14-2024 06:40 PM)Albert Chan Wrote:  Besides getting f(i) and slope, there are bookkeeping cost too.
Per timing, with less loops, my guess is Newtons beat Secant.

This is true for Newton that runs fewer iterations. But this benefit diminishes rapidly with more precise initial rate guesses, especially when we only use single precision or only 10 digits then only a few iterations suffice in both cases.

- Rob

[robve]

I've compiled a comparison as an addendum to my previous post.

The Secant code is as before, with the main parts being:

Code:

  i = iguess();
  if (i == 0.0)
    i = 1e-5;
  y = npmt();
  v = i;
  i /= 2.0;
  w = y;
  y = npmt();
  do {
    do {
      double t = i;
      h = (i - v)*y/(y - w);
      i -= h;
      v = t;
      if (i != 0.0) {
        w = y;
        y = npmt();
      }
      if (isnan(i))
        return;
    } while (fsign(y) == fsign(w) && fabs(y) < fabs(w));
  } while (y != 0.0 && y != w && g-- > 0);

The Newton code main parts, where I've made sure the convergence test doesn't quit early:

Code:

  i = iguess();
  if (i == 0.0)
    i = 1e-5;
  do {
    w = y;
    y = npmt();
    d = b*pmt - (pv + fv)*(1.0 + n*s/(r + r/j))/r - fv;
    h = 100*y/d;
    i -= h;
    if (isnan(i))
      return;
  } while (y != 0 && y != w && g-- > 0);

The npmt function:

Code:

double npmt()
{
  j = 0.01*i;
  k = 1.0 + j*b;
  l = log1p(j);
  r = expm1(-n*l);
  s = r + 1.0;
  return k*pmt - (pv + fv*s)*j/r;
}

To compare the root finding error of the rate between these methods, I've computed the new FV given the rate found and compared it to the TVM input FV. The absolute error of the FV is larger when the rate is less accurate.

Code:

printf(" %.2g\n", fabs(fv - (k*pmt*r/j-pv)/s))

The absolute error of FV is zero when the rate found is accurate, perhaps not to the last digit, but at least so it can be used to compute accurate TVM solutions for PV, PMT and FV.

Why FV? The FV appears to be more sensitive to rate errors than PV and PMT, i.e. those can have a zero error when FV has a nonzero error.

I will post updated results.


- Rob

[Albert Chan]

(06-14-2024 10:11 PM)robve Wrote:  Agree. But for Newton running on this 10 BCD machine the longer f(i)/f'(i) expression seems more likely to introduce floating point inaccuracies than the simpler Secant. Whereas the IEEE754 double precision version I have is more well behaved. It's one possible explanation.

Perhaps you missed my post regarding small rate branch. It is required to allow edge rate to work.
It is just simple quadratic f(i) ≈ f(0) + f'(0)*i + f''(0)/2*i²      (see LINE 20, FNZ)

f(0) = (pv+fv)/n + pmt
f'(0) = (pv-fv + (pv+fv)/n)/2
f''(0) = (pv+fv)/n * (n*n-1)/6

(06-07-2024 01:45 AM)Albert Chan Wrote:  10 INPUT "B,N,P,M,F? ";B,N,P,M,F
12 IF B THEN P=P+M @ F=F-M ! end mode
14 IF ABS(P)>ABS(F) THEN X=F @ F=-P @ P=-X @ N=-N ! time reversed
20 DEF FNZ(I,A,B,C) @ D=(A+B)/2 @ B=A*C @ FNZ=A+M+(D+B)*I @ D=D+2*B @ END DEF
30 DEF FNF(I) @ IF 1+N*I*I=1 THEN FNF=FNZ(I,(P+F)/N,P-F,(N*N-1)/12*I) @ END
40 S=EXPM1(LOGP1(I)*N) @ D=(P+F)/S*(1-N*(S+1)/(S+S/I))+P
50 FNF=((P+F)/S+P)*I+M @ END DEF
60 I=M/F @ Y=FNF(I) ! "smallest" edge guess
70 DISP 100*I,Y @ H=-Y/D @ I=I+H @ Y0=Y @ Y=FNF(I)
80 IF SGN(Y0)=SGN(Y) AND ABS(Y0)>ABS(Y) THEN 70
90 DISP 100*I,Y @ J=I-Y/D/2 @ DISP 100*J ! half-correction

BTW, I think your Newton's code termination test is wrong.
It does not look like HP71B LINE 80 at all!

do { ... } while (y != 0 && y != w && g-- > 0)

Quote:But it's almost free, right? Just iterate \( i_{\rm new} := i-(i-i_{\rm old})f/(f-f_{\rm old}) \) with one evaluation of npmt to get f(i) to update i.

Newton's method is not that bad. Also, (P+F)/S factor may also be shared.

FNF=((P+F)/S+P)*I+M
D=(P+F)/S*(1-N*(S+1)/(S+S/I))+P

Quote:This is true for Newton that runs fewer iterations. But this benefit diminishes rapidly with more precise initial rate guesses, especially when we only use single precision or only 10 digits then only a few iterations suffice in both cases.

Yes, I played with many such rate guess formulas.
The problem is they are not guaranteed to always work, unlike edge rate.

[robve]

(06-15-2024 03:39 AM)Albert Chan Wrote:  BTW, I think your Newton's code termination test is wrong.
It does not look like HP71B LINE 80 at all!

Your wish is my command! I changed the termination test to yours. The results are much worse:


  tvmperfcomp.pdf (Size: 36.85 KB / Downloads: 7)

EDIT: I should mention that the Newton results in the PDF above use i+h*h>i as was used in Albert's HP-71B code before. Actually, I also tried using fsign(y) == fsign(w) && fabs(y) < fabs(w) as shown n Albert's post above but then Newton only iterates once so the results are much worse. I don't understand why that works on the HP-71B. I'm probably doing something stupid.

- Rob

[robve]

I already had considered Newton termination tests and tried them before I posted. I don't post about considerations that don't work or work as well.

Disclaimer: it's experimental, so things can and probably will change when better alternative implementations or combinations emerge.

fabs(y) >= 1E-11 1E-12 does not terminate, but 1E-11 is too permissive with larger rate errors and does not offer termination guarantee

i + h*h > i this tests i% (scaled 100x), but still permits rate errors (h squared is too optimistic as a practical Newton convergence expectation) results in the PDF

j + h*h/10000 > j Albert's version with fractional rate (h is scaled 100x), but this is even worse than the results in the PDF, because of quicker quadratic downscaling of h

The three I'm more comfortable with are (edited to add context and the ones I've posted before):

  } while (y != 0.0 && y != w && g-- > 0);

    } while (fsign(y) == fsign(w) && fabs(y) < fabs(w));
  } while (y != 0.0 && y != w && g-- > 0);

    } while (y != 0.0 && fabs(y) < fabs(w));
  } while (y != 0.0 && y != w && g-- > 0);

where w is the previous y=npmt(). The second above repeats until sign change or step increases rather than decreases, then retries again a maximum of g times. The third is similar, but does stop when the sign changes, which I slightly prefer. We pick g >= 3 to be competitive to Secant (for Secant g does not matter too much, g=2 works just fine as in the code shown in previous posts and used in this comparison). Pick g too large and some test cases run too long without improving the rate with Newton.

With g=2, 3, 4 and 8 the Newton versus Secant results are:

  tvmperfcomp-g2348.pdf (Size: 145.98 KB / Downloads: 1)

Perhaps a better rate finder might be a hybrid:
1. make a good initial rate guess
2. do one or two Newton steps, if either finds the root then we're done
3. use the two points from 1) 2) to repeat Secant until convergence

- Rob

[Albert Chan]

Hi, robve

I think Newton's vs Secant's is just like integrate using Double-Exp vs Romberg.
Double-Exp integrate amazingly fast, using very few points, but required good f

Newton's method, with a single point, f and f' need to be accurate.
We only need accurate f/f', but I don't know how to get that without good f, f'.

Code:

npmt = k*pmt - (pv + fv*s)*j/r;      // k=1+j*b, s=(1+j)^-n, r=s-1

pv and fv normally are big, likely opposite sign, (pv+fv*s) may have cancellation issue.
But, (pv+fv), (pv-fv), (pv+pmt), (fv-pmt) probably OK

k = 1+j*b, for b=1, also may lose precision.

npmt = (1+j*b)*pmt - (pv + fv*(r+1))*j/r = pmt + j*b*pmt - j*fv - (pv+fv)*j/r

npmt = (-(pv+fv)/r - (fv-b*pmt))*j + pmt
npmt' = -(pv+fv)/r*(1+n*(r+1)/(r+r/j)) - (fv-b*pmt)

Or, if you prefer, its time-symmetric version, where R = subst(r, n=-n)
This setup coded for HP71B (with test samples), and seems to work well.

npmt = ((pv+fv)/R + (pv+b*pmt))*j + pmt
npmt' = (pv+fv)/R*(1-n*(R+1)/(R+R/j)) + (pv+b*pmt)

BTW, I would not do {j = 0.01*i; ...; i=100*j} for binary code.
I would do everything in j. If % is needed, simply multply by 100.

Sorry for the broken record, code need a branch to handle tiny rate!
Without that, edge guess is not safe to use.

A hard coded guess rate is not ideal. I remember Plus42 had hard coded i = 0.0001
Looking at the code, it is easy to make a situation where this would fail TVM solver.
This is why edge rate were later used, and user cannot supply their own guess.
Perhaps this hand holding is overkill, but it does make rate search robust.

[robve]

(06-15-2024 02:36 PM)Albert Chan Wrote:  Newton's method, with a single point, f and f' need to be accurate.
We only need accurate f/f', but I don't know how to get that without good f, f'.

Code:

npmt = k*pmt - (pv + fv*s)*j/r;      // k=1+j*b, s=(1+j)^-n, r=s-1

pv and fv normally are big, likely opposite sign, (pv+fv*s) may have cancellation issue.
But, (pv+fv), (pv-fv), (pv+pmt), (fv-pmt) probably OK

k = 1+j*b, for b=1, also may lose precision.

npmt = (1+j*b)*pmt - (pv + fv*(r+1))*j/r = pmt + j*b*pmt - j*fv - (pv+fv)*j/r

npmt = (-(pv+fv)/r - (fv-b*pmt))*j + pmt
npmt' = -(pv+fv)/r*(1+n*(r+1)/(r+r/j)) - (fv-b*pmt)

Or, if you prefer, its time-symmetric version, where R = subst(r, n=-n)
This setup coded for HP71B (with test samples), and seems to work well.

npmt = ((pv+fv)/R + (pv+b*pmt))*j + pmt
npmt' = (pv+fv)/R*(1-n*(R+1)/(R+R/j)) + (pv+b*pmt)

Thanks. I have already used, tested and compared (offline) the "original formula" as had mentioned in a previous post and didn't find it that different in terms of results. I did that to make sure I'm not injecting unnecessary calculation errors. I'm aware of such numerical problems and the differences pertaining cancellations and roundoff.

Sorry to disappoint, but Secant still wins with smallest FVerror (the losing values are marked red):

  tvmperfcomp-formulas.pdf (Size: 41.84 KB / Downloads: 3)

These results are with Newton g=4 but g=5,6 are the same as it bottoms out so I stuck to g=4 for now. The first column is Newton with the "simpler" formulas of npmt and npmt' I've used before, the next column is Newton with the above formulas in your post. The third column is Secant with the simpler npmt formula.

The edge case of a tiny rate to branch is not the culprit, but it does badly affect this single test case for Newton when Secant is OK:

n=480.000000 pv=100000.000000 pmt=-208.333333 fv=0.000000
=> I%=0% FVerror=540000

I also implemented a prototype hybrid algorithm. This uses g=3 (g=2 is not as good). Here compared to Secant (with g=2) it performs well:

  tvmperfcomp-hybrid.pdf (Size: 37.14 KB / Downloads: 2)

The number of evals is between Newton and Secant, which is to be expected. The FVerror looks alright, a bit better than Newton and a little less good than Secant. But this is preliminary. There is some room for improvement.

- Rob

[robve]

(06-15-2024 05:48 PM)robve Wrote:  I also implemented a prototype hybrid algorithm. This uses g=3 (g=2 is not as good). Here compared to Secant (with g=2) it performs well:

The number of evals is between Newton and Secant, which is to be expected. The FVerror looks alright, a bit better than Newton and a little less good than Secant. But this is preliminary. There is some room for improvement.

One nice aspect of Secant is that it does not have to deal with the problematic tiny rates where Newton goes haywire (see above and Albert's has extensively posted about this).

This updated hybrid method runs about as many npmt evaluations as Newton does, but has no issue with tiny rates and offers mostly higher accuracy (for these tests):

  tvmperfcomp-hybrid2.pdf (Size: 41.42 KB / Downloads: 8) (updated with a more robust version)

Newton takes 103 evaluations total over all tests, Secant takes 162 evaluations (though they are a bit cheaper), Hybrid takes 104 evaluations and these are cheaper than Newton.

The Hybrid method solver is pretty straightforward to compute rate i% (updated to a more robust version):

Code:

  g = 0; // no retries
  i = iguess();
  if (i == 0.0)
    i = 1e-5;
  v = i;
  w = y;
  y = npmt();
  d = b*pmt - (pv + fv)*(1.0 + n*s/(r + r/j))/r - fv;
  h = 100*y/d;
  i -= h;
  if (y == 0 || y == w || isnan(y))
    return;
  w = y;
  y = npmt();
  if (y == 0 || y == w || isnan(y))
    return;
  do {
    do {
      double t = i;
      h = (i - v)*y/(y - w);
      i -= h;
      v = t;
      if (i != 0.0) {
        rr = r; // only to compute FV error
        w = y;
        y = npmt();
      }
      if (isnan(i))
        return;
    } while (y != 0.0 && fabs(y) < fabs(w));
  } while (y != 0.0 && y != w && g-- > 0);
  if (y != 0) {
    i = v;
    r = rr;     // only to compute FV error
    s = r + 1;  // only to compute FV error
  }
  j = 0.01*i; // only to compute FV error
  printf("i%%=%.18g%%  FVerror=%.2g\n", fabs(fv - (k*pmt*r/j-pv)/s);

The Secant code is practically the same, but lacks the Newton step and instead uses i/2 for the next point, and with g=2 to repeat and retry if necessary.

EDIT: minor edit to correct testing nan(y) instead of nan(i) in the initialization, this can save one unnecessary npmt evaluation when y=nan.

- Rob

[Albert Chan]

(06-15-2024 09:52 PM)robve Wrote:  Newton takes 103 evaluations total over all tests, Secant takes 151 evaluations (though they are a bit cheaper), Hybrid takes 104 evaluations and these are cheaper than Newton.

From your C code, I am guessing this is done with IEEE double?
What is g? Why is it needed?

With lua tvm(), it finished it all with 132+12 = 144 evaluations (last 3 was duplicated tests)

24 test samples, f evaluations + final rate included (mode BEGIN → END, for first 3 cases)

Code:

rate_test = {
{40,i,900-40,-40,-1000+40       ,4, 0.04753367188571021},
{40,i,900+1000,1000,-1000-1000  ,2, -0.49999999999997724},
{40,i,900-400,-400,-1000+400    ,3, 0.8000000000098451},    -- x2
{31536000,i,0,-.01,3331667.0067 ,8, 1.1694311188017701e-07},
{525600,i,0,-.01,5282.367769    ,6, 1.902587523612461e-08},
{525600,i,0,-.01,5260.382426    ,6, 3.170979198140072e-09}, -- x2

{32,i,1-1e6,0,1e6               ,3, 3.125001611329216e-08}, -- x2
{360,i,0,-1e3,1e6               ,7, 0.004980365840090327},
{328,i,35e3,-324,0              ,5, 0.008720569148340681},
{58,i,-775,-49.56,4e3           ,6, 0.0026058787514755544},
{32,i,-6e3,0,1e4                ,6, 0.016091394922276913},
{48,i,19198,-450,0              ,6, 0.004918000645880568},

{5,i,-369494,17500,540000       ,6, 0.1200000561803591},
{348,i,243400,-1363.29,0        ,5, 0.00437501742139102},
{30,i,-3200,-717.44,60000       ,6, 0.048750016694561145},
{60,i,243400,-1363.29,-222975.98,4, 0.00437499988799056},
{24,i,0,-50.26,1281.34          ,7, 0.005209314607602606},
{6,i,-32000,0,28346.96          ,4, -0.019999978033474884},

{60,i,8000,-150.97,0            ,7, 0.004166696417110639},
{456,i,270000,-1215.333333,0    ,6, 0.0036443486419549583},
{480,i,100000,-208.333333,0     ,6, -6.652806639543341e-12},
{10,i,50,-30,400                ,7, 0.1443587132807996},
{10,i,50,-30,80                 ,5, -0.36893369874177734},
{10,i,-100,10,0                 ,7, 0},
}

Update: tested again, but with i=false instead of i=nil, to use big sized edge.
Again, it finished in 132+12 = 144 evaluations (last 3 was duplicated tests)

Basically the same solved rate, except for #22 and #23, with 2 sign changes. (0 or 2 roots)

Code:

#22 n=10 pv=50 pmt=-30 fv=400
small edge = -30/400 = -0.075 --> i =  0.1443587132807996
  big edge = -30/-50 =  0.60  --> i =  0.5317221326838473

#23 n=10 pv=50 pmt=-30 fv=80
small edge = -30/80 = -0.375  --> i = -0.36893369874177734
  big edge = -30/-50 =  0.60  --> i =  0.5846195526622099

[robve]

I'm not sure how to put this without offending anyone intellectually, but I think this needs to be said to explain some very basic concepts and my thoughts on the subject of solving TVM rates.

Solving is never as accurate as the floating point precision MachEps (or ULP) suggests. No matter how we rewrite the algebraic formulations of the function and its derivative. Inaccuracy increases with increasing number of terms in the expressions of the function to solve and its algebraic or numerical derivative. Sure, some forms are more accurate than others in preventing cancelations and unnecessary rounding, but with TVM the input parameter space is large and unconstrained (except for i% and n).

When solving roots we cannot expect nor should we visualize function \( f(x) \) to cut a perfect skinny line through the x-axis at which there is a single point \( x \) that is the root \( f(x)=0 \) to find. No! Instead, it is a fat line with some thickness \( \varepsilon \). This thickness depends on the inaccuracy of evaluating \( y=f(x) \). This fat line cuts the x-axis with an interval \( [x-\varepsilon,x+\varepsilon] \) on which all \( x \) values return \( f(x)=0 \). If \( f(x) \) becomes shallow close to the root, then the fat line may project an even wider interval (any deviation in the height \( y \) moves \( x \) faster along the axis). Root finders like Newton and Secant jump to the next point \( x \) closer to the root based on the behavior (slope) of the function and its derivative at the previous point \( x \). But we only see a fat line at the previous point \( x \), not a skinny line. Getting accurate values of \( f(x) \) and \( f'(x) \) to stay on course to converge to the root is problematic with a fat line. What is the value of \( f(x) \) exactly? And what is the slope \( f'(x) \) at this point? Since we cannot get highly accurate answers, we may end up jumping away from the root and far enough to cause an overcorrection in the next step to jump to the other side of the root and so on, possibly ad infinitum.

Let me explain my thoughts on using a hybrid approach, that is, iterate with Secant after a few (one or two) Newton TVM rate solver steps. Theoretically Newton is a better method, but suffers more than Secant with "fat lines" as we know for TVM: tiny rates cause Newton to go crazy and lower floating point precision causes non-termination such as on 10 digit BCD calculators. (Disclaimer: unless more advanced termination conditions are used that strike a good balance between execution cost and accuracy.) Furthermore, claiming that Newton is always quadratically or even better to converge is nonsense in practice. Initially the first steps are converging faster for y=npmt(i%) toward the root. But this is not guaranteed to continue to the root and in practice there is a significant convergence slow down closer to the root or even overshooting and subsequent correcting.

Using Secant after Newton is conceptually similar to solving higher-order polynomial equations. To solve these, Euclidean division can be used to remove factor \( x-r \) after finding one root \( r \). This repeats until all roots were found. But this is not stable. Accuracy suffers greatly with higher polynomial orders and less precise arithmetic. Therefore, Weierstrass aka Durand–Kerner must be used as a final step to improve the roots iteratively (see e.g. NR and other sources).

Finally, the i% rate solution we find may slightly vary in practice and still be correct. There is no single value that is the only correct one. There is a range of correct i% values (the fat line, remember). It depends on the floating point precisions and on the given TVM parameters how wide this range of i% values is that are all roots of the TVM equation. Therefore, be very skeptical when comparing decimal places in the rate solution i%. Because TVM is a model, any i% that maintains TVM equilibrium is correct. Checking the TVM model just requires computing N, PV, PMT or FV given the solution we found for the rate i%. So after checking over two dozen TVM cases, I found that FV is the most sensitive to i%, i.e. \( {\rm error} = | FV - FV(i) | \) is nonzero for errors in i% when the errors in \( PV(i) \) and \( PMT(i) \) are sometimes zero. The "FV error" is a good metric to evaluate the accuracy of the root i%. Let's not use \( N(i) \) because the direct formula of N has an expression involving logarithms and terms that may introduce inaccuracies that have nothing to do with how accurate the root i% is.

Hope this helps.

- Rob

[robve]

I've updated the Hybrid solver code in my post #50 above. This version is a bit more robust. I realized I could use the previous rate stored in v when it overshoots. Results in the PDF are also updated.

- Rob

[dm319]

Fascinating discussion! Must confess I don't understand everything (understatement), but I think Rob you are saying that the 'correct' solve for i may not be a precise number given the impreciseness of representing the other values (PV/FV etc) using limited precision?

So what I'm about to say is going to sound like I'm ignoring all that...

I'm blaming bxparks for sending me on a rabbit hole of trying to programmatically measure accuracy of these solvers. It has involved all sorts of trouble, including use of AWK (did you know it has arbitrary precision BTW?), discovering that tidyverse isn't compatible with the gnu multiple precision library, but finding that data.table is. And also learning a bit about what 'accuracy' is.

Anyway, I think I have finally arrived at some results which look somewhat correct. I have taken absolute accuracy for interest rate and relative accuracy for all other values.

abs accuracy = -log10(abs(ref - result))
rel accuracy = -log10(abs(ref - result) / abs(ref))

I've given a +3 bonus for precise result (Inf accuracy) and a -3 penalty for failing to come to a solution. I think these results are truer than my 'eyeball, and round in my head the number of digits of precision', which was prone to errors that I kept finding as I added more data.

Code:

|calculator                 |    1|   1b|    2|    3|    4|   5A|   5B|   6A|   6B|    7|  8|    9|   10|   11|   12| acc_mean|
|:--------------------------|----:|----:|----:|----:|----:|----:|----:|----:|----:|----:|--:|----:|----:|----:|----:|--------:|
|R Rmpfr                    | 16.8| 36.1| 36.3| 36.7|   NA| 35.3| 35.0| 35.5| 35.3| 35.5| 42| 36.4| 36.4|   NA|   NA|     34.8|
|NSTK TVM V8                | 16.8| 33.1| 33.2| 32.7| 21.4| 32.3| 32.0| 32.5| 32.3| 32.5| 39| 33.4| 33.4| 28.4| 36.1|     31.3|
|C47 eadc4507f              | 16.8| 32.7| 33.3| 30.6| 18.4| 32.3| 35.0| 32.5| 32.3| 25.0| 37| 31.1| 30.9| 27.7| 33.1|     29.9|
|Plus42 1.1.9               | 16.8| 33.1| 32.4| 33.7| 21.4| 31.8|   NA|   NA| 32.3| 31.4| 36| 29.8| 30.3| 31.0| 23.7|     29.5|
|NSTK TVM V5                | 16.8| 31.0| 33.1| 25.2|  9.3| 31.3| 31.7| 31.4| 32.3| 19.8| 32| 33.4| 25.7|   NA|   NA|     27.2|
|NSTK TVM V6                |  1.3|  3.0| 33.1| 25.2| 21.4| 32.3| 32.0| 32.5| 32.3| 32.5| 39| 29.8| 30.3|   NA|   NA|     26.5|
|TI-83 Plus                 | 13.8| 12.9| 12.8| 13.3| 21.4| 11.7|   NA|   NA| 11.9| 20.4| 19|  5.6|  5.2| 10.8|  1.0|     12.3|
|RPN83P v0.9.0              | 13.1| 12.9| 13.4| 13.4| 21.4| 12.0| 11.3| 11.4| 11.5| 11.4| 13|  5.6|  5.2| 11.7|  1.9|     11.3|
|Casio FC-200               |  9.6|  9.3|   NA| 10.6| 12.1|  5.7|  8.3|  5.4|  5.4| 19.7| 12|   NA|   NA|   NA|   NA|      9.8|
|HP-17BII                   | 12.4| 10.7| 11.7| 12.2| 11.1|  5.7|  8.3|  5.4|  5.4| 11.7| 12|   NA|   NA|   NA|   NA|      9.7|
|HP-17B                     | 12.4| 10.7| 11.7| 12.2| 11.1|  5.7|  8.3|  5.4|  5.4| 11.9| 12|  6.6|  6.2|   NA|   NA|      9.2|
|TI BA II Plus Professional |  9.6|  9.5|  8.7| 10.6| 13.2|  8.7|   NA|   NA|  8.4| 16.7| 12|  4.6|  4.5|  9.5|  3.3|      9.2|
|HP-12cp                    | 12.4| 11.0| 11.7| 12.2| 11.1|  5.7|  8.3|  5.4|  5.4| 21.4| 12|  3.2|  3.5| 11.0|  2.8|      9.2|
|HP-20b                     | 12.4| 10.7| 11.7| 12.2| 11.1|  5.7|  8.3|  5.4|  5.4| 11.7| 12|  6.6|  6.2| 10.9|  2.8|      8.9|
|HP-30b                     | 12.4| 10.7| 11.7| 12.2| 11.1|  5.7|  8.3|  5.4|  5.4| 11.7| 12|  6.6|  6.2| 10.9|  2.8|      8.9|
|Sharp PC-1403H HPF 188177  |  9.6|  8.5|  8.8| 10.6| 21.4|  8.7|   NA|  5.4|  5.4|  9.7| 15|  1.7|  1.2|  6.1|  1.0|      8.1|
|HP-12c comma               |  9.6|  9.5| 11.7|  7.2|  6.0|  8.7|   NA|  8.4|   NA|  6.7| 12|  3.2|  3.2|   NA|   NA|      7.8|
|HP-92                      |  9.6|  8.5|  8.8| 10.6| 21.4|  5.7|  8.3|  5.4|  5.4|  9.7| 12|  1.7|  1.6|  4.3|  1.2|      7.6|
|Sharp EL-735               |  9.6|  8.3|  8.7|  3.5| 21.4|  5.7|  8.3|  5.4|  5.4|  9.7|  9| 12.3|  2.7|  1.2|  1.6|      7.5|
|HP-12c                     |  9.6|  8.6|  8.8| 10.6|  9.1|  5.7|  8.3|  5.4|  5.4| 19.7| 12|  1.7|  1.6|  4.6|  1.0|      7.5|
|HP-37E                     |  9.6|  8.5|  8.8| 10.6|  9.0|  5.7|  8.3|  5.4|  5.4|  9.7| 12|  1.7|  1.6|  6.0|  1.0|      6.9|
|HP-38E                     |  9.6|  8.5|  8.8| 10.6|  9.0|  5.7|  8.3|  5.4|  5.4|  9.7| 12|  1.7|  1.6|  6.0|  1.0|      6.9|
|HP-27                      |  9.6|  9.5|  8.8|  9.0|  9.0|  5.7|  8.3|  5.4|  5.4|  6.7|  6|  1.5| -1.8|  4.3| -2.0|      5.7|
|HP-22                      |  9.6|  8.6|  8.8|  9.0|  9.0|  5.7|  8.3|  5.4|  5.4|  6.7|  6|  1.7| -1.8|  4.7| -2.0|      5.7|
|Sharp PC-1421              |  9.6|  8.2|  8.7|  3.5|  6.0|  8.7|   NA|  8.4|   NA|  9.7|  9|  1.7|  1.6| -9.0|  1.6|      5.2|
|HP-70                      |  9.3|  6.6|  8.8|  1.2|  9.0|  5.7|  8.3|  5.4|  5.4|  6.7|  6|  1.7| -1.8| -1.1| -2.0|      4.6|
|HP-80                      |  9.6|  5.1|  9.4|  1.2|  6.0|  5.7|  8.3|  5.4|  5.4|  6.7|  6|  1.3| -1.8| -1.0| -2.0|      4.3|

I do have pretty graphs but they will have to wait until I'm happy these results are robust.

More importantly: does anyone have an HP-41 with financial pack? Also I'm missing the Prime. Wouldn't mind a bit more Casio too while I'm at it.

[Albert Chan]

(06-16-2024 05:18 PM)robve Wrote:  When solving roots we cannot expect nor should we visualize function \( f(x) \) to cut a perfect skinny line through the x-axis at which there is a single point \( x \) that is the root \( f(x)=0 \) to find. No! Instead, it is a fat line with some thickness \( \varepsilon \). This thickness depends on the inaccuracy of evaluating \( y=f(x) \). This fat line cuts the x-axis with an interval \( [x-\varepsilon,x+\varepsilon] \) on which all \( x \) values return \( f(x)=0 \). If \( f(x) \) becomes shallow close to the root, then the fat line may project an even wider interval (any deviation in the height \( y \) moves \( x \) faster along the axis). Root finders like Newton and Secant jump to the next point \( x \) closer to the root based on the behavior (slope) of the function and its derivative at the previous point \( x \). But we only see a fat line at the previous point \( x \), not a skinny line. Getting accurate values of \( f(x) \) and \( f'(x) \) to stay on course to converge to the root is problematic with a fat line. What is the value of \( f(x) \) exactly? And what is the slope \( f'(x) \) at this point? Since we cannot get highly accurate answers, we may end up jumping away from the root and far enough to cause an overcorrection in the next step to jump to the other side of the root and so on, possibly ad infinitum.

Nicely put. BTW, fat line may not be smooth, due to rounding errors.

We should stop Newton's iteration when f(x) cannot be trusted.
This is why Plus42 does final half-correction, to reduce error radius.
The idea is: f is likely no good, but perhaps we can trust its sign?

The best we could hope for is solved x actually get f(x) = 0
Since f=npmt, should we use npmt(final i) as accuracy metric?
I would think round-trip for PMT is what we wanted, not FV.

My setup tried to make line skinny, with time-symmetry.

Code:

local s = EFF(i,n)
local k = (pv+fv)/s
local d = k*(1-n*(s+1)/(s+s/i)) + pv
f = (k + pv)*i + pmt
eps = -f / d

Note the + pv in f and d

If pv is big relatively to fv, we can apply time-symmetry --> n, pv, fv = -n, -fv, -pv
For real life financial situations, many times pv or fv is 0, thus make the line very skinny!

lua> n,i,pv,pmt,fv = 48,nil,19198,-450,0 -- test sample #12

For reference, Plus42 to get I%YR: 4.918000645880570070062687539208835e-1

find_rate() has time-symmetry flip disabled (tvm() does the flipping)

lua> find_rate(n,i,pv,pmt,fv, true)
0.023439941660589644    220.49758413166887
0.006730726899179227    19.372332869813135
0.004940893361022922    0.24155367517261084
0.004918004408642906    3.96964566675706e-05
0.004918000645880669    1.0231815394945443e-12
0.004918000645880572    -5.684341886080802e-14
0.004918000645880574

By design, f(i) would not overshoot, if calculated accurately.
The fact that it does, suggested f(i) is not very sharp.

For this example, even sign of final f is wrong!
Final half-correction make it worse, errors = 5 ULP.
But at least it is better than full correction, error = 8 ULP

If we apply Secant's method from last 2 points, error = 8 ULP too.
Based on this example alone, Secant's method is not better.
Secant's slope is no where near its true value.

secant slope  = 11117.714285714286 -- from last 2 points, almost touching
newton slope = 10549.817399409634 -- from lua code last d value
correct slope  = 10549.817399409620 -- mpmath with mp.dps=34

lua> find_rate(-n,i,-fv,pmt,-pv, true)
0.023439941660589644    220.49758413166887
0.006730726899179217    19.37233286981302
0.004940893361022924    0.24155367517272452
0.004918004408642897    3.969645661072718e-05
0.004918000645880665    1.0231815394945443e-12
0.004918000645880568    0
0.004918000645880568

With time-symmetry, we have PV=0, which make f(i) sharper. No overshoot, errors = -2 ULP

Comment 6/17/24:
After some thought, this does not show time-symmetry version is better.
The few ULP's gain in accuracy may be simply lucky

Overshoot can happen by accumulated rounding errors.
Even sharp f may overshoot, even though technically it shouldn't.

[Albert Chan]

I just updated tmv() code. Verbose option now gives i, f, fp

(06-17-2024 02:33 PM)Albert Chan Wrote:  lua> tvm(10,nil,1000,-287,2000,true)

Code:

-0.1435                 116.13771101577868      -1087.5631610761484
-0.036712904618008635   29.240213109966078      -536.0217925151969
0.017837513887852915    7.537006560375005       -262.9585107578723
0.04649985510150145     1.9377897954382775      -129.06227911365863
0.06151423272404261     0.5045985944708491      -62.25029703788141
0.0696201955624815      0.14253250628991054     -27.205477349325747
0.07485930621875686     0.05843692701773762     -4.949693275931963
0.08666547758333192     0.29079749268407795     44.04186343829406
nan

When I try this no solution case with TI BA-35 Solar, it hanged

[Werner]

(06-16-2024 05:18 PM)robve Wrote:  I'm not sure how to put this without offending anyone intellectually [..]

Absolutely no need to be defensive, very well explained, Rob!
Werner

[SlideRule]

from an unknown but astute source;

EXCEL uses a slightly different TVM equation in solving for interest rate … making it most useful

PV(1+RATE)^NPER + PMT (1+RATE*TYPE)[{(1+RATE)^NPER} - 1] / RATE + FV = 0

This equation holds true i.e. it equals ZERO when at least one or at most two of the three variables of PV, FV, and PMT are negative.

BEST!
SlideRule

[robve]

(06-17-2024 06:37 PM)SlideRule Wrote:  from an unknown but astute source;

EXCEL uses a slightly different TVM equation in solving for interest rate … making it most useful

PV(1+RATE)^NPER + PMT (1+RATE*TYPE)[{(1+RATE)^NPER} - 1] / RATE + FV = 0

This equation holds true i.e. it equals ZERO when at least one or at most two of the three variables of PV, FV, and PMT are negative.

BEST!
SlideRule

Assuming TYPE is 0 for END and 1 for BGN, then we recognize PV and FV are switched and their signs inverted and N is inverted:
PV*S + PMT*K*R/I + FV = 0
with K=1+B*I and R=EXPM1(-N*LNP1(I)) and S=R+1

Compare this to

FV = (K*PMT*R/I-PV)/S
K*PMT*R/I - PV - FV*S = 0

This gives the same TVM problem.

EDIT: I just checked with the Hybrid solver:

./tvm 328 '?' 35000 -324 0
n=328.000000 pv=35000.000000 pmt=-324.000000 fv=0.000000
evals=6, I%=0.872056914834067998%, FVerror=1.3e-10

./tvm -328 '?' 0 -324 -35000
n=-328.000000 pv=0.000000 pmt=-324.000000 fv=-35000.000000
evals=5,I%= 0.872056914834067887%, FVerror=0

So this aligns with the observation Albert made.

- Rob

[robve]

(06-17-2024 09:03 PM)robve Wrote:  Assuming TYPE is 0 for END and 1 for BGN, then we recognize PV and FV are switched and their signs inverted and N is inverted:
PV*S + PMT*K*R/I + FV = 0
with K=1+B*I and R=EXPM1(-N*LNP1(I)) and S=R+1

Compare this to

FV = (K*PMT*R/I-PV)/S
K*PMT*R/I - PV - FV*S = 0

This gives the same TVM problem.

This way we can control the error in npmt() and optimize the result, which Albert mentioned before i.e. "time reversal".

Lets call this the forward npmt():
j = 0.01*i;
k = 1.0 + j*b;
l = log1p(j);
r = expm1(-n*l);
s = r + 1.0;
return k*pmt - (pv + fv*s)*j/r;

The backward dual npmt() is:
j = 0.01*i;
k = 1.0 + j*b;
l = log1p(j);
r = expm1(n*l);
s = r + 1.0;
return k*pmt + (fv + pv*s)*j/r;

The npmt() version to pick to solve the rate depends on the values of s and r that we get. My thought on this is that we should pick the version in which |s| and |r| are relatively close.

Let me illustrate with the "time reversed" TVM problem, which is essentially the same. The forward npmt() is best for this TVM example (the backward npmt() line is fatter):
./tvm 5 '?' -369494.09 17500 540000
3/0 i%= 11.9999997132 y=-2.05654941965e-08
4/0 i%= 11.9999997132 y= 0
4, 11.9999997132285614%, 0
s=0.567426862982966185 r=-0.43257313701703376

./tvm -5 '?' -540000 17500 369494.09
3/0 i%= 11.9999997132 y=-2.05764081329e-08
4/0 i%= 11.9999997132 y=-7.27595761418e-12
5/0 i%= 11.9999997132 y= 2.91038304567e-11
5, 11.9999997132285596%, 0
s=1.76234166063797959 r=0.762341660637979479

The |s| and |r| values are closer in the forward than in the backward npmt().

Another TVM example:
./tvm 328 '?' 35000 -324 0
3/0 i%= 0.87205693754 y= 6.9640698257e-06
4/0 i%= 0.872056914834 y= 1.54045665113e-11
5/0 i%= 0.872056914834 y= 5.68434188608e-14
6/0 i%= 0.872056914834 y=-5.68434188608e-14
6, 0.872056914834067998%, 1.3e-10
s=0.0579632092841858571 r=-0.942036790715814143

./tvm -328 '?' 0 -324 -35000
3/0 i%= 0.87205693754 y= 6.96406976886e-06
4/0 i%= 0.872056914834 y= 1.54614099301e-11
5/0 i%= 0.872056914834 y= 0
5, 0.872056914834067887%, 0
s=17.2523228501226313 r=16.2523228501226313

Because we don't know the accurate rate yet to solve for, we can't decide ahead to pick the forward or backward npmt() based on the rate and n. But a good guess (which we have) will help!

It might also be interesting to see how this affects N, PV, PMT and FV computations, wben we aren't using the rate solver. We can pick the forward or backward npmt() to get parameters j, k, l, r, s and use the corresponding N, PV, PMT, and FV formulas to compute those as follows:

For the forward npmt():
n = log((k*pmt-fv*j)/(k*pmt+pv*j))/l);
pv = k*pmt*r/j-fv*s;
pmt.= (fv*s+pv)*j/k/r;
fv = (k*pmt*r/j-pv)/s;

For the backward npmt():
n = log((k*pmt-fv*j)/(k*pmt+pv*j))/l); // same as forward
pv = -(k*pmt*r/j+fv)/s
pmt.= -(pv*s+fv)*j/k/r;
fv = -k*pmt*r/j+pv*s;

- Rob