Threaded Mode | Linear Mode

Dave Britten · 03-30-2021, 04:58 PM

In light of all the recent Sharp enthusiasm, here's my attempt at a TVM solver for the PC-1211, and compatible models like the PC-1250, and no doubt many others (and their corresponding Tandy/Radio Shack versions).

The math is based largely on the 17BII manual's formulas, this TVM program for serveral different HPs, and some of the mathematical notes found here.

The PC-1211 and PC-1250 don't have hyperbolic functions for easy, accurate calculation of e^x-1 for |x| close to zero, but using the infinite series for ln(1+x) does offer some accuracy improvement.

Usage is fairly straight-forward. On the PC-1211, you'll need to switch to DEF mode. On other models like the PC-1250, there's a dedicated DEF key that you press instead of Shift. Set B=1 for begin mode, otherwise set B=0. These are the key assignments - for storing variables, first enter the value into the display, then press the appropriate key:

Shift Z: Store n
Shift X: Store i%
Shift C: Store PV
Shift V: Store PMT
Shift B: Store FV

Shift A: Compute n
Shift S: Compute i%
Shift D: Compute PV
Shift F: Compute PMT
Shift G: Compute FV

Code:

1 "Z":AREAD N

2 PRINT "N=";N:END

3 "X":AREAD I

4 PRINT "I%=";I:END

5 "C":AREAD P

6 PRINT "PV=";P:END

7 "V":AREAD M

8 PRINT "PMT=";M:END

9 "B":AREAD F

10 PRINT "FV=";F:END

11 "A":GOSUB 100:N=LN ((MK/J-F)/(MK/J+P))/LN (1+J):GOTO2

21 "S":IF I=0LET I=.01

22 GOSUB 28:IF Y=0GOTO 4

23 V=I:I=I+.01:GOSUB 28

24 T=I:I=I-(I-V)Y/(Y-W):V=T:GOSUB 28

25 IF (Y<>0)*(Y<>W)GOTO 24

26 GOTO 4

28 W=Y:GOSUB 100:Y=P+KMU+FS:RETURN

31 "D":GOSUB 100:P=-KMU-FS:GOTO 6

41 "F":GOSUB 100:M=(-FS-P)/K/U:GOTO 8

51 "G":GOSUB 100:F=-(MKU+P)/S:GOTO 10

100 J=.01I:S=(1+J)^-N:U=(1-S)/J:K=1+JB:RETURN

Alternative lines for greater accuracy at the cost of speed:

100 J=.01I:IF ABS J>=1 LET L=LN (1+J):GOTO 105

101 S=J:C=2:X=-J:Z=X

102 Z=Z*X:L=S-Z/C:IF L=S GOTO 105

103 C=C+1:S=L:GOTO 102

105 S=1/EXP (NL):U=(1-S)/J:K=1+JB:RETURN

Results are generally pretty good, and don't seem to differ by more than a few thousandths of a penny in the few sample problems that I tried. I'm not sure if the Sharps use binary rather than decimal, but I did notice that solving for PMT in a problem where FV=0, then immediately solving for FV could yield small non-zero values on the order of 1E-5.

robve · 03-30-2021, 05:07 PM

Thanks for sharing!

(03-30-2021 04:58 PM)Dave Britten Wrote: I'm not sure if the Sharps use binary rather than decimal, but I did notice that solving for PMT in a problem where FV=0, then immediately solving for FV could yield small non-zero values on the order of 1E-5.

FYI. SHARP pocket computers use 10 digits BCD to display values and store values in variables and 12 digits internally to evaluate expressions.

For example (SQR is sqrt):

SQR 2 * SQR 2
2

X=SQR 2
X*X
1.999999999

There is a small penalty to use variables to store (temporary) results. It is sometimes beneficial to "inline" expressions to preserve or increase accuracy.

- Rob

Dave Britten · 04-01-2021, 04:44 PM

(03-30-2021 05:07 PM)robve Wrote: Thanks for sharing!

(03-30-2021 04:58 PM)Dave Britten Wrote: I'm not sure if the Sharps use binary rather than decimal, but I did notice that solving for PMT in a problem where FV=0, then immediately solving for FV could yield small non-zero values on the order of 1E-5.

FYI. SHARP pocket computers use 10 digits BCD to display values and store values in variables and 12 digits internally to evaluate expressions.

Good to know!

Interestingly, I tested the same TVM program on the EL-5500III that I received in the mail today (this model supports two-character variable names, so explicit * had to be inserted in most of the implied multiplications), and it seems to have better accuracy. Example:

N=360
I=4.75/12 (3.958333333E-01)
PV=141000
FV=0

Solve for PMT, on both machines:
PMT=-735.5227446

Then immediately solve for FV without entering anything else.

PC-1250A (Tandy PC-3):
FV=3.731560805E-05

EL-5500III:
FV=0.

robve · (This post was last modified: 05-29-2024 01:55 AM by robve.)

(04-01-2021 04:44 PM)Dave Britten Wrote: Example:

N=360
I=4.75/12 (3.958333333E-01)
PV=141000
FV=0

Solve for PMT, on both machines:
PMT=-735.5227446

Then immediately solve for FV without entering anything else.

PC-1250A (Tandy PC-3):
FV=3.731560805E-05

EL-5500III:
FV=0.

I also get FV=0 with a PC-1350 with PMT=-735.5227446 when B=0 (pay at end) and with B=1 (pay at begin) I get PMT=-732.6227796 reported for DEF-V. This value is also computed with the alternative 100-105 lines.

The S-BASIC interpreter and the CPU on the EL-5500III and PC-1350 are (almost) the same, the forensics: http://www.rskey.org/~mwsebastian/miscprj/results.htm see also http://basic.hopto.org/basic/manual/basic-compare.pdf

The BASIC of the PC-1211 is an earlier S'-BASIC. The PC-1250 is derived from that early version, see: https://sharppocketcomputers.com/#basic

- Rob

robve · 05-29-2024, 09:17 PM

(03-30-2021 04:58 PM)Dave Britten Wrote: In light of all the recent Sharp enthusiasm, here's my attempt at a TVM solver for the PC-1211, and compatible models like the PC-1250, and no doubt many others (and their corresponding Tandy/Radio Shack versions).

In light of all the recent exciting discussions about TVM contributions, I looked back at your excellent post of the TVM program for SHARPs as I wanted to play around with it a bit more and compare to HPs etc.

I found an issue with the program where the rate calculation would not converge, which might also affect similar programs on other machines with 10 digit precision. The problem is roundoff. The Newton iteration keeps updating the estimated rate I but Y (the root) never reaches zero for this case:

B=1 (begin mode)
N=40
PMT=-40
PV=900
FV=-1000

The rate should be 4.753367189% but it won't get there because Y (the root) hovers around 0.000000171.

Coming up with an effective and safe convergence criterium can be tricky. A possible workaround is to compare the estimated rate to the previously estimated rate:

25 IF Y<>W AND 100-V/I<>99 GOTO 24

where 100-V/I equals 99 if V is close to I up to 9 digits. It's a bit of a risky move since we could hit a rate plateau that is not close to the root (so perhaps test if we are at least somewhat close to the root also?). Anecdotally, this change appears to work well, though I haven't tried a lot of cases and the accuracy of the estimated rate may be slightly lower, because I get 4.753367191% instead of 4.753367189%.

Perhaps the following is a worthy improvement:

25 IF Y<>W AND 10-V/I<>9 GOTO 24

This check works if V is close to I when V<I but not when V>I. But this might in fact still be alright when the Newton iteration hops back and forth around the root as we would actually expect to observe with Newton.

This version gives the exact rate 4.753367189%. Changing 10-V/I<>9 to 10-I/V<>9 (a forced "hop around the root") converges with 4.753367188%.

I've also slightly improved and compacted your program by computing with DEF-C followed by DEF-N or DEF-J or DEF-V or DEF-M or DEF-F to compute N, I%, PV, PMT or FV, respectively:

Code:

' TVM

' To switch begin mode on/off: DEF-B

' To enter values: DEF-N or DEF-J or DEF-V or DEF-M or DEF-F

' To calculate: press DEF-C (beep) then DEF-N or DEF-J or DEF-V or DEF-M or DEF-F

10 "B" B=B=0 : PRINT "BGN=";MID$("NY",B+1,1) : END

11 "C" C=C=0 : BEEP C : END

12 "N" AREAD N : IF C GOSUB 39 : N=LN((M*K/J-F)/(M*K/J+P))/LN(1+J)

13 PRINT "N=";N : END

14 "J" AREAD I : IF C GOSUB 30

15 PRINT "I%=";I : END

16 "V" AREAD P : IF C GOSUB 39 : P=-K*M*U-F*S

17 PRINT "PV=";P : END

18 "M" AREAD M : IF C GOSUB 39 : M=(-F*S-P)/K/U

19 PRINT "PMT=";M : END

20 "F" AREAD F : IF C GOSUB 39 : F=-(M*K*U+P)/S

21 PRINT "FV=";F : END

30 IF I=0 LET I=.01

31 GOSUB 36 : IF Y=0 RETURN

32 V=I,I=I+.01 : GOSUB 36

33 T=I,I=I-(I-V)*Y/(Y-W),V=T : GOSUB 36

34 IF Y<>W AND 10-V/I<>9 GOTO 33

35 RETURN

36 W=Y : GOSUB 39 : Y=P+K*M*U+F*S : RETURN

39 J=.01*I,S=(1+J)^-N,U=(1-S)/J,K=1+J*B,C=0 : RETURN

I left out the alternative to compute S=(1+J)^-N with the series for LN1(J)=LN(1+J) for small J.

- Rob

dm319 · 05-31-2024, 01:11 PM

(05-29-2024 09:17 PM)robve Wrote: In light of all the recent exciting discussions about TVM contributions

If this relates to any of my posts, then I am quoting this to my wife.

(05-29-2024 09:17 PM)robve Wrote: B=1 (begin mode)
N=40
PMT=-40
PV=900
FV=-1000

Interesting, I wonder if this will throw a problem to any of the solvers. Seems to get there ok on my 12cp, but it does seem to take a while...

Albert Chan · (This post was last modified: 06-01-2024 05:07 PM by Albert Chan.)

(05-29-2024 09:17 PM)robve Wrote: 36 W=Y : GOSUB 39 : Y=P+K*M*U+F*S : RETURN

Instead of NPV = 0, solve NPMT = 0 for rate may be better.

10 B=1 @ N=40 @ P=900 @ M=-40 @ F=-1000
20 DEF FNF(I) @ S=(1+I)^(-N) @ U=(1-S)/I @ K=1+B*I
30 FNF=P+K*M*U+F*S @ END DEF ! NPV
40 V=.01 @ W=FNF(V) @ I=.02 @ Y=FNF(I)
50 H=I-V @ V=I @ I=I-H*Y/(Y-W) @ W=Y @ Y=FNF(I)
60 DISP 100*I,Y @ IF V/I+I/V<>2 THEN 50

>run
3.55877069046      -223.49575008
4.3407684864        -68.532154113
4.68660450765      -10.546204073
4.74950349013      -.604877167
4.75333056057      -.005731135
4.75336716844      -.000003149
4.75336718857      -.000000001

>30 FNF=(P+F*S)/U+K*M @ END DEF ! NPMT
>run
4.7203131853        -.2913219635
4.75310902569      -.0022751338
4.75336716699      -.0000001902
4.75336718857      -.0000000001

For rate convergence, I use ratio, r + 1/r = 2
If r = 1-ε, r+1/r = 2 + ε^2 + ε^3 + ε^4 + ...

robve · 06-01-2024, 01:07 PM

(05-31-2024 10:39 PM)Albert Chan Wrote: Instead of NPV = 0, solve NPMT = 0 for rate may be better.

Nice. This seems to work better. Do you have a reference to an article or implementation?

(05-31-2024 10:39 PM)Albert Chan Wrote: 60 DISP 100*I,Y @ IF V/I+I/V<>2 THEN 50

The termination test V/I+I/V<>2 may not be accurate when solving for NPV=0 or for NPMT=0. For example, solving NPV=0 in 10 digits BCD precision for the TVM example gives I=4.753367173 and previous guess V=4.753337033 because V/I=.9999936592 and I/V=1.000006341 sum to 2.000000000 with 10 digits. Solving for NPMT=0 gives I=4.753367191 and V=4.753366979, with I being very close though.

But when I try 1+V/I<>2 I get the exact rate I% when solving for NPV=0 or for NPMT=0. Also with 10-V/I<>9 or the similar 9+V/I<>10 to shift out irrelevant last digit to mitigate rounding errors. Or ABS((V-I)/I)<1E-9 which checks for a specified relative error 1E-9. Or check if |Y|<E for Y with error E but this won't terminate when we desire a reasonable E=1E-9 or even E=1E-8 to solve for Y=NPV=0 with 10 digits BCD evaluation (i.e. on a decent calc). Now, the test set is too small to make sweeping conclusions. It just illustrates that the termination criterium is important to obtain some level of accuracy while avoiding non-termination. Bisection terminates, which might be a better but slower alternative.

- Rob

Albert Chan · (This post was last modified: 06-01-2024 05:17 PM by Albert Chan.)

Hi, robve

Perhaps make rate search true Newton's method?

10 B=1 @ N=40 @ P=900 @ M=-40 @ F=-1000
20 DEF FNF(I) @ S=(1+I)^(-N) @ U=(1-S)/I @ K=1+B*I
25 D=(P+F)/(1-S)*(1-N*S/(1-S+U))+(B*M-F)
30 FNF=(P+F)/U-F*I+K*M @ END DEF ! FNF=NPMT
40 I=.01 @ Y=FNF(I)
50 H=-Y/D @ I=I+H @ Y=FNF(I) ! Newton's method
60 DISP 100*I,Y @ IF I+H*H<>I THEN 50

>run
4.707227645        -.4066650546
4.75336179241     -.0000475552
4.75336718857       0

And time-reversed setup. (barring rounding errors, same FNF(I))

>15 VARSWAP P,F @ P=-P @ F=-F @ N=-N
>run
4.70722764501     -.4066650545
4.75336179241     -.0000475552
4.75336718857       0

Albert Chan · 06-01-2024, 09:47 PM

(06-01-2024 01:07 PM)robve Wrote: The termination test V/I+I/V<>2 may not be accurate when solving for NPV=0 or for NPMT=0 ...

I rather early termination than never terminate.

Perhaps perform an extra Newton's correction after V/I+I/V==2 ?
Code already has points (V,W), (I,Y). Why waste them?

Albert Chan · 06-01-2024, 10:13 PM

Plus42 uses NPMT formula and Newton's method for rate too.

Guess is picked so that each iteration, f always get better (closer to zero)
Termination criteria is that if f=0 or it change sign, or f not getting better.

Since ULP's are finite, even if f is at plateau, it always terminate.

https://raw.githubusercontent.com/thomas...mmands9.cc

Code:

if (p_isnan(f))

    return ERR_NO_SOLUTION_FOUND;

if (f == 0 || (f > 0) != (f0 > 0))

    break;

if (fabs(f) >= fabs(f0)) {

    if (i == i - eps / 1000) {

        i -= eps / 2;

        break;

    } else

        return ERR_NO_SOLUTION_FOUND;

}

Albert Chan · 06-02-2024, 06:37 PM

(05-29-2024 09:17 PM)robve Wrote: B=1 (begin mode)
N=40
PMT=-40
PV=900
FV=-1000

We can make this to END mode, by shift up payment: {PV,FV} ← {PV+PMT, FV-PMT}

N=40, PV=860, PMT=-40, FV=-960

We can split up into 2 loans (combined, we get back original loan)

#1: N=40, PV=860-x, PMT=-40, FV=-(860-x)
#2: N=40, PV=x, PMT=0, FV=(860-x)+(-960) = -(x+100)

Unlike original loan, splitted loans have have direct formula for rate.

From 1st loan, rate = pmt/(x-pv) = 40/(860-x)
From 2nd loan, x*(1+rate)^N = x - (PV+FV)

--> x*((1+rate)^N-1) + (PV+FV) = 0
--> f = x + (PV+FV) / ((1+rate)^N-1) = 0      // flatten curve to make Newton's method efficient

If x=0, 1st loan rate = 40/860 = 4.65%, 2nd loan rate = ∞ --> (x>0, rate > 4.65%)

f(x=10) = 10 + (-100) / ((1+40/850)^40-1) = -8.89363
f(x=20) = 20 + (-100) / ((1+40/840)^40-1) = 1.58000

Interpolate for f=0, we get x = 18.4914 --> rate = 4.75337% (all digits correct)

Below code assumed pmt ≠ 0 (we have direct rate formula for pmt=0)
We use smaller size edge rates [pmt/-pv, pmt/fv] as guess, same as Plus42.
For this example, rate guess = -40 / (-1000+40) ≈ 4.16667%

10 B=1 @ N=40 @ P=900 @ M=-40 @ F=-1000
15 IF B THEN P=P+M @ F=F-M ! Now, B=0
20 DEF FNF(X) @ D=X-P @ I=M/D @ S=(1+I)^N
25 D=1 - (P+F)/(S-1)^2 * N*S/(1+I) * -I/D
30 FNF=X+(P+F)/(S-1) @ END DEF
40 X=(ABS(P)<ABS(F))*(P+F) @ Y=FNF(X)
50 H=-Y/D @ X=X+H @ Y=FNF(X)
60 DISP 100*I,X,Y @ IF X+H*H<>X THEN 50

>run
4.75245128279      18.32908914         -.1698427457
4.75336718645      18.4912666946      -.000000398
4.75336718859      18.4912670746      .0000000071

Interestingly, converge rate seems to match NPMT formula with Newton's method.

robve · (This post was last modified: 06-04-2024 02:11 AM by robve.)

(06-01-2024 10:13 PM)Albert Chan Wrote: Plus42 uses NPMT formula and Newton's method for rate too.

Guess is picked so that each iteration, f always get better (closer to zero)
Termination criteria is that if f=0 or it change sign, or f not getting better.

Since ULP's are finite, even if f is at plateau, it always terminate.

Indeed, that is an interesting approach. We could improve this a bit by iterating one more time after detecting a sign change or when f does not get better. We could also compare f to f0 to use the best rate corresponding to f or f0 whichever is closer to zero.

In the Sharp code this would be:

Code:

10 "B" B=B=0 : PRINT "BGN=";MID$("NY",B+1,1) : END

11 "C" C=C=0 : BEEP C : END

12 "N" AREAD N : IF C GOSUB 39 : N=LN((M*K/J-F)/(M*K/J+P))/LN(1+J)

13 PRINT "N=";N : END

14 "J" AREAD I : IF C GOSUB 30

15 PRINT "I%=";I : END

16 "V" AREAD P : IF C GOSUB 39 : P=K*M*(S-1)/J-F*S

17 PRINT "PV=";P : END

18 "M" AREAD M : IF C GOSUB 39 : M=(F*S+P)*J/K/(S-1)

19 PRINT "PMT=";M : END

20 "F" AREAD F : IF C GOSUB 39 : F=(K*M*(S-1)/J-P)/S

21 PRINT "FV=";F : END

' solve Y=NPMT=0 for rate I

30 IF I=0 LET I=.01

31 GOSUB 38 : IF Y=0 RETURN

32 G=1,V=I,I=I+.01 : GOSUB 38

33 T=I,I=I-(I-V)*Y/(Y-W),V=T : GOSUB 38

34 IF SGN Y=SGN W IF ABS Y<ABS W GOTO 33

35 IF Y<>0 IF W<>0 IF G LET G=0 GOTO 33

36 IF ABS Y>ABS W LET I=V

37 RETURN

38 W=Y : GOSUB 39 : Y=(P+F*S)*J/(1-S)+K*M : RETURN

' compute

39 J=.01*I,S=(1+J)^-N,K=1+J*B,C=0 : RETURN

where Y and W are f and f0, respectively.

With this, a more accurate rate is obtain such as for PMT=-40 for the example. For PMT=-400 the rate is 80% and we get very close with 80.00000002%.

I've also changed the code to substitute U=(1-S)/J inline, because when entering the rate I%=80 to compute FV we get a very bad FV=0 instead of FV=-1000 because of catastrophic cancellation in U=(1-S)/J. With the change we get FV=-893.8463964 instead of -1000. As we all know, this catastrophic cancellation can be avoided if we use LN1(x) and EXP1(x), which aren't built-in on SHARPs so we should roll out our own series computations, i.e. the alternative part that Dave has in his code computes LN1 but EXP1 is also required. This is a nice example why this matters (a lot).

- Rob

EDIT: correct typo

Albert Chan · 06-02-2024, 11:18 PM

(06-02-2024 06:37 PM)Albert Chan Wrote: rate = pmt/(x-pv)
f = x + (PV+FV) / ((1+rate)^N-1) = 0 // flatten curve to make Newton's method efficient
...
Interestingly, converge rate seems to match NPMT formula with Newton's method.

I just checked. Barring rounding errors, both ways are equivalent!

I reuse variables from HP71B code. Let I = initial rate, I2 the next iteration.
H is newton correction for X (i.e. X --> X+H)

I = M/(X-P)
I2 = M/((X-P)+H) = I/(1+H/(X-P)) = 1/(1/I+H/M)

From I, we have X = M/I+P

Y = FNF(X)
Y*I = X*I + (P+F)/(S-1)*I = (P+F) * I/(S-1) + P*I + M

This is just NPMT formula, time-reversed!
Let Z = Y*I = NPMT formula

(dZ/dI)/Z
= d(ln(Y*I))/dI
= 1/I + (dY/dI)/Y                      -- H = -Y/(dY/dX)
= 1/I + (-1/H) * (dX/dI)            -- (dY/dI) = (dY/dX) * (dX/dI)
= 1/I + (-1/H) * (-M/I^2)          -- (dX/dI) = d(M/I+P)/dI
= (I+M/H)/I^2

I - Z/(dZ/dI) = 1/(1/I+H/M)       -- = I2. Doing dY/dX vs dZ/dI are equivalent!

Albert Chan · 06-03-2024, 11:38 AM

(06-02-2024 09:46 PM)robve Wrote: We could improve this a bit by iterating one more time after detecting a sign change or when f does not get better. We could also compare f to f0 to use the best rate corresponding to f or f0 whichever is closer to zero.

Not necessarily.

NPMT formula is strictly concave or convex function.
see TVM solve for interest rate, revisited, post 7, 24, 25.

If we start rate from the edge, NPMT Newton convergence is always 1 sided.
Each iteration is an improvement over previously, with quadratic convergence.

Both situations you mentioned technically should not happen.
The fact that they do suggested f and f0 numbers are suspect.

If we overshoot with f having opposite sign, this is as good as we can get.

Regarding f or f0, both same sign, but rate seems to converge, we don't know which is better.
That's why Plus42 pick the center of 2 "best" rate, to keep error radius small.

(05-31-2024 10:39 PM)Albert Chan Wrote: 10 B=1 @ N=40 @ P=900 @ M=-40 @ F=-1000
20 DEF FNF(I) @ S=(1+I)^(-N) @ U=(1-S)/I @ K=1+B*I
30 FNF=P+K*M*U+F*S @ END DEF ! NPV
40 V=.01 @ W=FNF(V) @ I=.02 @ Y=FNF(I)
50 H=I-V @ V=I @ I=I-H*Y/(Y-W) @ W=Y @ Y=FNF(I)
60 DISP 100*I,Y @ IF V/I+I/V<>2 THEN 50

>run
3.55877069046      -223.49575008
4.3407684864        -68.532154113
4.68660450765      -10.546204073
4.74950349013      -.604877167
4.75333056057      -.005731135
4.75336716844      -.000003149
4.75336718857      -.000000001

Above example, NPV=0 for rate, we have 1 sided convergence.
But we may lose this nice feature, with time-reversed setup.

>15 VARSWAP P,F @ P=-P @ F=-F @ N=-N
>run
11.3570289363      -36923.539379
2.35993812666      1387.63023153
2.68581267787      1287.38125489
6.87064131886      -3585.77746844
3.79135227821      769.04731268
4.33514429757      379.62865086
4.86526534953      -115.24756927
4.7418099072        11.55575603
4.75306056619      .30740914
4.75336803852      -.00085219
4.75336718851      .00000007

NPV=0 for rate does not handle well with negative N
For the same reason, NFV=0 for rate have problem with positive N.

robve · (This post was last modified: 06-04-2024 12:15 AM by robve.)

(06-03-2024 11:38 AM)Albert Chan Wrote:
(06-02-2024 09:46 PM)robve Wrote: We could improve this a bit by iterating one more time after detecting a sign change or when f does not get better. We could also compare f to f0 to use the best rate corresponding to f or f0 whichever is closer to zero.

Not necessarily.

First, I neglected to add zero tests in the updated version for f and f0 that get us to the root immediately. Second, finite precision especially BCD isn't behaving as smoothly. If we're not exactly at the root yet, it may not hurt to do a step when we keep the previous best guess to compare to later. Something like this:
34 IF SGN Y=SGN W IF ABS Y<ABS W GOTO 33 ' convergence?
35 IF Y<>0 IF W<>0 IF G LET G=0 GOTO 33 ' try one more time
36 IF ABS Y>ABS W LET I=V ' pick better one
37 RETURN

This is one way I get better rates up to the last digit (10 digit BCD) that is off less than without this approach. But it may not always work. It appears to produce results that are often "retractable", meaning computing PMT with the computed rate gives the PMT value that we started with up to one ulp.

- Rob

robve · (This post was last modified: 06-04-2024 12:51 AM by robve.)

I flipped the test lines 35 and 36 to avoid an issue when retrying. The point is really about picking a better rate when reasonably possible:

34 IF SGN Y=SGN W IF ABS Y<ABS W GOTO 33
35 IF Y<>0 IF W<>0 IF G LET G=0 : GOTO 33
36 IF ABS Y>ABS W LET I=V
37 RETURN

I actually also played a bit with the idea of a weighted average, something like:

35 IF Y<>0 IF W<>0 IF G LET G=0 GOTO 33
36 IF Y<>0 OR W<>0 IF SGN Y<>SGN W LET I=(I*ABS Y+V*ABS W)/(ABS Y+ABS W)
37 RETURN

But it's not working as well as I had expected.

Note that SHARP IF ... IF ... LET is more efficient and faster than IF ... AND ... LET.

EDIT: OK, just after writing this I see right away why it's not working well: if Y=0 then I=V when I should be picked. Duh... So it should be:
36 IF SGN Y<>SGN W LET I=(I*ABS W+V*ABS Y)/(ABS Y+ABS W) and maybe get rid of the extra iteration, although I get better results for some cases with it, such as I%=80 for PMT=-400 which is spot on.

- Rob

Dave Britten · 06-04-2024, 07:33 PM

Interesting stuff, Rob! I will confess that when I wrote this, I mostly tested it with very banal TVM problems, and I think maybe the standard "Susan problem", but didn't go so far as to find exotic edge cases that could break it. So I'm not too surprised someone has been able to find a few. Smile

robve · 06-04-2024, 11:35 PM

Albert, we're mixing methods and various different code snippets which becomes a bit of a confusing mess. So I was mistaken to believe that a sign change suffices even with one more iteration I had suggested (but "not necessarily", do you mean more than one or none?). This is not right, to just stop after a sign change, at least not for the method based on the earlier post but modified to solve NPMT=0 rather than NPV=0. Counter examples are not hard to find.

Take for example:
B=1
PMT=1,000
PV=900
FV=-1,000
N=40
which requires two more iterations to get I%=-50% correct on the SHARP (10 digits BCD). Milage may vary depending on precision and float representations.

Another example "how to be a millionaire":
B=0
PMT=-1,000
PV=0
FV=1,000,000
N=360
which gives I%/period=0.498036584 but that takes nine more iterations after the sign change to get there.

Granted, this is with solving NPMT=0 combined with the original approach to iterate until sign change f and f0 and as long as the rate improves i<i0 which is:

32 G=9,V=I,I=I+.01 : GOSUB 38
33 T=I,I=I-(I-V)*Y/(Y-W),V=T : GOSUB 38
34 IF SGN Y=SGN W IF ABS Y<ABS W GOTO 33

Setting G=9 as above allows up to 9 more iterations (or more if i<i0 again) to always terminate. We also want to avoid Y=W that errors out (the slope hits zero). And when the final iteration produces opposite f and f0 signs we can take a weighted average as a final interpolation of the root:

35 IF Y<>0 IF W<>0 IF Y<>W IF G LET G=G-1 : GOTO 33
36 IF SGN Y<>SGN W LET I=(I*ABS W+V*ABS Y)/(ABS Y+ABS W)
37 RETURN

This runs fairly efficiently, as far as vintage 80s CMOS 8 bit calculators are concerned i.e. a couple of seconds at most.

I know, there are many other ways to implement a solver with Newton, Halley, Ostrowski or others. On the other hand, it is nice to have a small and simple program to get accurate answers while guaranteeing termination.

- Rob

Albert Chan · 06-05-2024, 01:05 AM

Hi, robve

Conditions required to make sign change suffices, without more iterations.

1. use NPMT formula (equivalent splitted loan method also OK)
2. use true derivative, not secant line slope.
3. use edge rate, in END mode, either [M/-P, M/F]

I modified my 2nd post (#9), with Plus42 "smaller" edge as guess.
I also make calculations more accurate, just to show it is hard to have f overshoot.

10 INPUT "B,N,P,M,F? ";B,N,P,M,F @ B=B*M
20 DEF FNF(I) @ S=-EXPM1(LOGP1(I)*-N)
25 D=(P+F)/S*(1-N*(1-S)/(S+S/I))+(B-F)
30 FNF=(P+F)*I/S+(B-F)*I+M @ END DEF ! FNF=NPMT
40 X=(ABS(P)<ABS(F))*(P+F)-B @ I=M/(X-P) @ Y=FNF(I)
50 H=-Y/D @ I=I+H @ Y=FNF(I) ! Newton's method
60 DISP 100*I,Y @ IF I+H*H<>I THEN 50

Code:

>run

B,N,P,M,F? 1,40,900,-40,-1000

 4.7524512828        -.0080716937

 4.75336718645       -.0000000187

 4.75336718857        0

>run

B,N,P,M,F? 1,40,900,-400,-1000

 79.9999999644       -.000000183

 80.000000001         0  

>run

B,N,P,M,F? 1,40,900,1000,-1000

-50                   0

>run

B,N,P,M,F? 0,360,0,-1000,1e6

 .313028391458        504.47068745

 .47127071031         62.96905733

 .497374290777        1.51994785

 .498036166886        .00095669

 .498036584011        0