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