An improved version of the TVM program for SHARPs.
I kept it compact (22 lines of BASIC) and made it faster than the previous version. It takes only a few seconds to compute interest rates (faster than the SHARP EL-735). This version is also fully tested (on a SHARP PC-1403H).
I wanted to keep it small while returning accurate results, at least as accurate as possible with a 10 digit BASIC calculator that has limitations. This update uses LNP1 and EXPM1 formulas by Dieter and the initial interest guess by Albert. Also I%=0 is now permitted as input and is returned when the computed interest is zero or very small. This required a couple of tricks to avoid computation errors. There is no error trapping. When a computation fails then the result is ERR. Some explanations are included in the program.
How to use:
- Switch begin mode on/off: DEF-B
- Enter values: DEF-N or DEF-J (I%) or DEF-V (PV) or DEF-M (PMT) or DEF-F (FV)
- Calculate: press DEF-C (beep) then DEF-N or DEF-J or DEF-V or DEF-M or DEF-F
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 : IF N=0 LET N=LN((K*M-F*J)/(K*M+P*J))/L
13 PRINT "N=";N : END
14 "J" AREAD I : IF C GOSUB 29
15 PRINT "I%=";I : END
16 "V" AREAD P : IF C GOSUB 39 : P=K*M*R/J-F*S
17 PRINT "PV=";P : END
18 "M" AREAD M : IF C GOSUB 39 : M=(P+F*S)*J/K/R
19 PRINT "PMT=";M : END
20 "F" AREAD F : IF C GOSUB 39 : F=(K*M*R/J-P)/S
21 PRINT "FV=";F : END
' 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 secant method, pick initial nonzero guess
30 G=3,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 I=I+(I=0)*1E-5 : GOSUB 37 : V=I,I=I/2 : GOSUB 37
' secant loop with count down to ensure termination
33 T=I,I=I-(I-V)*Y/(Y-W),V=T : GOSUB 37
34 IF SGN Y=SGN W IF ABS Y<ABS W GOTO 33
35 IF Y<>0 IF Y<>W IF G LET G=G-1 : GOTO 33
' round very small I% to zero
36 I=(ABS I>1E-9)*I : RETURN
' compute Y=NPMT and save previous Y in W
37 W=Y : GOSUB 40 : Y=K*M-(P+F*S)*J/R : 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
Notes:
- line 30: a retry counter G is still necessary on this calculator, but G=3 suffices instead of G=9 with the improved initial guess
- line 31: we must avoid initializing I%=0 even when pv+fv+n*pmt is nonzero, so we assign I%=1E-5
- line 36: small |I%|<1E-9 are zero, based on experiments, this threshold may vary on other machines' MachEps
- line 39: if I%=0 then assign variables so that P, M, F and N are computed correctly with J/R=-1/N and K=R=S=1 and direct formula for N while avoiding div by zero when M=0
Add these four lines for amortization with
period DEF-A after entering a TVM problem:
Code:
' amortization for given period and next, accumulated int and prn are T and V, remaining balance R
50 "A" AREAD A : J=.01*I,R=P-F,T=0,V=0
51 FOR K=1 TO N : S=R*J,Q=-M-S,R=R-Q
52 IF K>=A LET T=T+S,V=V+Q : PRINT STR$ K;" ";INT(100*S+.5)/100;"int ";INT(100*Q+.5)/100;"prn"
53 NEXT K : END
The results for the 12 TVM test examples look alright to me:
Code:
| # | Ref | N | I%YR | PV | PMT | FV | P/YR | End |
|----|------------|--------------|-----------|----------|--------------|-----------|------|-----|
| 1 | DM | 38 x 12 | 5.25% | 270'000 | ? | 0 | 12 | yes |
PMT=-1368.148355
| 1b | DM | 38 x 12 | ? | 270'000 | -14'584/12 | 0 | 12 | yes |
I%YR=4.373218369
| 2 | SlideRule | 360 | 15% → 12% | 100'000 | ?-? | 0 | 12 | yes |
-1264.444022 - -1028.612597 = -235.831425 => PV=6803.092173
| 3 | Kahan 1983 | 60x60x24x365 | 10% | 0 | -0.01 | ? | =N | yes |
FV=331667.0067
| 4 | DM | 480 | 0 → ? | 100'000 | ?→ PMT | 0 | 12 | yes |
I%=0
| 5 | Dieter | 10 | ? | 50 | -30 | 400 | 1 | yes |
I%=14.43587133
| 6 | Dieter | 10 | ? | 50 | -30 | 80 | 1 | yes |
I%=ERR
| 7 | A Chan | 10 | ? | -100 | 10 | 1e-10 | 12 | yes |
I%=0
| 8 | Miguel | 32 | ? | -999'999 | 0 | 1e6 | 1 | yes |
I%=3.12500161E-06
| 9 | DM | ? | 25 | 100000 | -2083.333334 | 0 | 12 | yes |
N=1040.639029 this one is not that good due to limited precision of this machine
| 10 | DM | ? | 25 | 100000 | -2040.816327 | 0 | 12 | no |
N=1139.586234 also this is is not that good due to limited precision of this machine
| 11 | robve | 60x24x365 | 1/6% → ? | 0 | -0.01 | ?→ FV | =N | yes |
FV=5260.382453 => I%N=1.666673857E-01
| 12 | robve | 40 | ? → I%YR | 900 | -400 | -1000 → ? | 1 | no |
I%=80 => FV=-900
Results for the C implementation of the same program with IEEE754 double precision, where I% is returned (to get I%YR multiply I% by 12 when P/YR=12 or multiply by N for tests 3 and 11). To compute the rate, the number of npmt function evaluations for root finding and retry counter G initialized to 2 are shown as evals/G. The retry counter G matters for TVM problems that dance around the root, such as n=40 pv=900 pmt=1000 fv=-1000 in BGN mode and n=360 pv=0 pmt=-1000 fv=1e+06 (not shown):
Code:
======== 1
TVM n=456 i%=0.4375% pv=270000 pmt=0 fv=0
pmt=-1368.1483553505775
======== 1b
TVM n=456 i%=0% pv=270000 pmt=-1215.33 fv=0
i%=0.182956 y=-341.776
3/2 i%= 0.364314888635 y= -0.244321551002
4/2 i%= 0.36444462776 y= 0.0198834207638
5/1 i%= 0.36443486391 y=-9.14124257179e-07
6/0 i%= 0.364434864359 y=-3.41060513165e-12
7/0 i%= 0.364434864359 y= 0
i%=0.364434864359157562%
======== 2
TVM n=360 i%=1.25% pv=100000 pmt=0 fv=0
pmt=-1264.44402156504339
TVM n=360 i%=1% pv=100000 pmt=0 fv=0
pmt=-1028.61259692550448
TVM n=36 i%=1.25% pv=0 pmt=-235.831 fv=0
pv=6803.09216198562535
======== 3
TVM n=3.1536e+07 i%=3.17098e-07% pv=0 pmt=-0.01 fv=0
fv=331667.006690777082
======== 4
TVM n=480 i%=0% pv=100000 pmt=0 fv=0
pmt=-208.333333333333343
TVM n=480 i%=0% pv=100000 pmt=-208.333 fv=0
i%=0%
======== 5
TVM n=10 i%=0% pv=50 pmt=-30 fv=400
i%=7.17794 y=5.88612
3/2 i%= 14.4161676678 y= 0.0120456394605
4/2 i%= 14.4310106992 y= 0.00296966661476
5/2 i%= 14.4358673527 y= 2.42831893971e-06
6/2 i%= 14.4358713273 y= 4.90221196969e-10
7/2 i%= 14.4358713281 y= 0
i%=14.4358713280799602%
======== 6
TVM n=10 i%=0% pv=50 pmt=-30 fv=80
i%=-50.9714 y=10.8304
3/2 i%= nan y= nan
i%=nan%
======== 7
TVM n=10 i%=0% pv=-100 pmt=10 fv=1e-10
i%=9.09172e-12 y=5.00044e-12
3/2 i%= 1.81802059632e-11 y= 0
i%=1.81802059632393454e-11%
======== 8
TVM n=32 i%=0% pv=-999999 pmt=0 fv=1e+06
i%=1.5625e-06 y=0.015625
3/2 i%= 3.12500161133e-06 y= 0
i%=3.12500161132921604e-06%
======== 9
TVM n=0 i%=2.08333% pv=100000 pmt=-2083.33 fv=0
n=1060.30340666380812
======== 10
TVM n=0 i%=2.08333% pv=100000 pmt=-2040.82 fv=0 BGN
n=1076.31957397602309
======== 11
TVM n=525600 i%=3.17098e-07% pv=0 pmt=-0.01 fv=0
fv=5260.38242600032663
TVM n=525600 i%=0% pv=0 pmt=-0.01 fv=5260.38
i%=1.58549e-07 y=4.1684e-06
3/2 i%= 3.17097919838e-07 y= 1.73472347598e-18
4/2 i%= 3.17097919838e-07 y= 0
i%=3.17097919837667697e-07%
======== 12
TVM n=40 i%=0% pv=900 pmt=-400 fv=-1000 BGN
i%=25.7127 y=-271.439
3/2 i%= 79.9993961671 y= -0.00301916932335
4/2 i%= 79.9999999949 y=-3.02858325085e-08
5/2 i%= 80.000000001 y=-1.13686837722e-13
6/2 i%= 80.000000001 y= 0
i%=80.0000000009845138%
TVM n=40 i%=80% pv=900 pmt=-400 fv=0 BGN
fv=-999.999523662601518
The C program and a script to run the tests:
Code:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
double fsign(double x) { return copysign(x != 0, x); }
double n, pv, pmt, fv, i, b = 0; // TVM parameters
double j, k, l, r, s; // dependent variables
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;
}
void comp()
{
if (i == 0.0) {
j = k = s = 1.0;
r = -n;
if (n == 0.0)
n = -(pv + fv)/pmt;
return;
}
npmt();
}
double iguess()
{
// https://www.hpmuseum.org/forum/thread-12941-post-160018.html#pid139867
double a = (fv + pv)/n;
double b = (fv - pv) - a;
double c = (pmt + a)/b;
b = (n*n-1)*a*c/b;
return 100*c*(b - 3)/(b - 1.5);
}
void solve()
{
int m = 1; // count number of evals
int g = 2; // retries after overshooting
double v, w, y;
if (pv + fv + n*pmt == 0) {
i = 0;
return;
}
i = iguess();
if (i == 0)
i = 1e-5;
y = npmt();
v = i;
i /= 2.0;
w = y;
y = npmt();
++m;
printf("i%%=%g\ty=%g\n", i, y);
do {
do {
double t = i;
i -= (i - v)*y/(y - w);
v = t;
if (i != 0) { // to avoid I%=0 (not included in the BASIC version)
w = y;
y = npmt();
}
printf("%2d/%d i%%=%18.12g y=%18.12g\n", ++m, g, i, y);
if (isnan(i))
return;
} while (fsign(y) == fsign(w) && fabs(y) < fabs(w));
} while (y != 0.0 && y != w && g-- > 0);
if (fabs(i) < 1e-16)
i = 0;
}
int main(int argc, char **argv)
{
n = strtod(argv[1], NULL);
i = strtod(argv[2], NULL);
pv = strtod(argv[3], NULL);
pmt = strtod(argv[4], NULL);
fv = strtod(argv[5], NULL);
if (argc > 6)
b = strtod(argv[6], NULL) != 0.0;
printf("\nTVM n=%g i%%=%g%% pv=%g pmt=%g fv=%g", n, i, pv, pmt, fv);
if (b)
printf(" BGN");
printf("\n");
if (*argv[2] == '?') {
solve();
printf("i%%=%.18g%%\n", i);
} else {
comp();
if (*argv[1] == '?')
printf("n=%.18g\n", (n != 0.0 ? n : log((k*pmt-fv*j)/(k*pmt+pv*j))/l));
else if (*argv[3] == '?')
printf("pv=%.18g\n", k*pmt*r/j-fv*s);
else if (*argv[4] == '?')
printf("pmt=%.18g\n", (fv*s+pv)*j/k/r);
else if (*argv[5] == '?')
printf("fv=%.18g\n", (k*pmt*r/j-pv)/s);
}
}
Code:
./tvm 40 ? 900 -40 -1000 1
./tvm 40 ? 900 1000 -1000 1
./tvm 40 ? 900 -400 -1000 1
./tvm 31536000 ? 0 -0.01 331667.0067
./tvm 525600 ? 0 -0.01 5282.36776893622391
./tvm 525600 ? 0 -0.01 5260.38242600032663
./tvm 32 ? -999999 0 1000000
./tvm 360 ? 0 -1000 1000000
echo
echo ======== 1
./tvm 456 0.4375 270000 ? 0
echo
echo ======== 1b
./tvm 456 ? 270000 -1215.3333333333 0
echo
echo ======== 2
./tvm 360 1.25 100000 ? 0
./tvm 360 1 100000 ? 0
./tvm 36 1.25 ? -235.83142463954 0
echo
echo ======== 3
./tvm 31536000 3.1709791983765e-07 0 -0.01 ?
echo
echo ======== 4
./tvm 480 0 100000 ? 0
./tvm 480 ? 100000 -208.333333333333343 0
echo
echo ======== 5
./tvm 10 ? 50 -30 400
echo
echo ======== 6
./tvm 10 ? 50 -30 80
echo
echo ======== 7
./tvm 10 ? -100 10 1e-10
echo
echo ======== 8
./tvm 32 ? -999999 0 1e6
echo
echo ======== 9
./tvm ? 2.08333333333333333 100000 -2083.333334 0
echo
echo ======== 10
./tvm ? 2.08333333333333333 100000 -2040.816327 0 1
echo
echo ======== 11
./tvm 525600 3.1709791983765e-07 0 -0.01 ?
./tvm 525600 ? 0 -0.01 5260.38242600032663
echo
echo ======== 12
./tvm 40 ? 900 -400 -1000 1
./tvm 40 80.0000000009844996 900 -400 ? 1
- Rob