PC-1211, PC-1250, etc. TVM

[robve]

(06-11-2024 10:32 PM)dm319 Wrote:  

(06-11-2024 07:13 PM)robve Wrote:  An improved version of the TVM program for SHARPs.

Great work!! Very impressive to fit all this in these lines of code. Also a very interesting performance profile. Eyeballing it, very close to the HP-37, sometimes identical. I will put it into the results table and I'm curious to run the dimension reduction on it again and see where it pops up. I'd be curious of a video of it in action solving tvm problems if you ever find the time.

Thanks. I still want to give the Albert's suggested modified Newton's method a try to compare efficiency of TVM rate calculations, but haven't had the time yet. Working in non-IEEE754 is a pain. Some things just don't work out numerically as we're stuck with 10 digits to store in variables (internal calculations are 12 digits). Results and behaviors of the C implementation don't necessarily translate, as can be expected, hence the challenge. The current method finds the root with a few NPMT evaluations. I was surprised that some roots are found in just 3 evaluations as you can see in the results.

Now, I'm curious how the SHARP PC-1421 fares on these TVM problems. I don't have access to one at the moment, but could try an emulator. It's a machine that supports TVM operations on the keypad but also directly within BASIC, i.e. setting TVM variable PV=x is just an assignment. Computing is done with e.g. COMP PMT that updates the PMT variable. It also has cash flow and IRR calculations. This should be very handy to write specialized financial applications (which is not my thing). It seems that the PC-1421 was not popular, perhaps because financial people don't write BASIC? Wrong target audience? The machine is a bit weird and that's why it intrigues me.

- Rob

[nickapos]

(06-09-2024 01:22 PM)rprosperi Wrote:  

(06-08-2024 09:21 PM)robve Wrote:  [snip]

If you want another bewildering TVM problem:
B=1 (begin mode)
N=40
PV=900
PMT=1000
FV=-1000
What's the interest rate?

In the real world this makes no sense, i.e. getting rich for nothing in return, but mathematically it is a valid TVM problem to solve. Try it on a calculator. Is the rate you get a whole number?

- Rob

Using an HP-27S, I get -50 exactly.

Prime v1 with latest firmware also returns -50

[dm319]

(06-12-2024 12:29 AM)robve Wrote:  Working in non-IEEE754 is a pain. Some things just don't work out numerically as we're stuck with 10 digits to store in variables (internal calculations are 12 digits).

I'm curious - can you elaborate on this? In C I am presuming you are using binary rather than decimal. When you mean IEEE754 do you mean floating point binary or decimal? I've always thought doing maths stuff in binary is a pain as we generally view the results in decimal, which makes it hard to determine accuracy in terms of the number of decimal points.

(06-09-2024 01:22 PM)rprosperi Wrote:  Using an HP-27S, I get -50 exactly.

In keeping with the emulated HP-30b I guess. I suspect the HP-27s is similar to the HP-17b and onwards, but I wouldn't mind more data...

(06-12-2024 05:20 AM)nickapos Wrote:  Prime v1 with latest firmware also returns -50

Exactly? That is a little disappointing.

[Albert Chan]

(07-11-2022 09:08 PM)Albert Chan Wrote:  I just realized there is no financial function called NPMT (we do have NPV, NFV)
We can define NPMT as a TVM function that work the same way, if time direction reversed.

NPMT = C*pv + C(n=-n)*fv + n*pmt

NPMT tends to be more "straight" than (NPV, NFV), with values inside (NPV, NFV)
Thus, f = NPMT/n = 0 is great for solving rate x, using Newton's method.

Just want to clear up meaning of NPMT. It is not "Net Payment" ... more like "N Payments".
Function used for rate search is f = NPMT/n (sometimes refer as npmt function)

Technically, time-symmetry is exactly what you expected when viewing film backwards:

{n,i,pv,pmt,fv} ↔ {-n,i,fv,-pmt,pv}

NPMT is time-symmetrical (same value with time-symmetry transform)
npmt is not. Because of n denominator, we had to negate it to preserve the same result.

{n,i,pv,pmt,fv} ↔ {-n,i,fv,-pmt,pv} --> negate --> {-n,i,-fv,pmt,-pv}

Warning: To count sign changes for possible solutions, we should use (n*pmt)
Time symmetry will give same sign changes either way, even if negated.

{pv,n*pmt,fv} ↔ {fv,(-n)*(-pmt),pv} = {fv,n*pmt,pv}

[robve]

(06-12-2024 12:45 PM)dm319 Wrote:  

(06-12-2024 12:29 AM)robve Wrote:  Working in non-IEEE754 is a pain. Some things just don't work out numerically as we're stuck with 10 digits to store in variables (internal calculations are 12 digits).

I'm curious - can you elaborate on this? In C I am presuming you are using binary rather than decimal. When you mean IEEE754 do you mean floating point binary or decimal? I've always thought doing maths stuff in binary is a pain as we generally view the results in decimal, which makes it hard to determine accuracy in terms of the number of decimal points.

It's about the internals. IEEE754 binary representations and operations must meet a set of requirements to ensure accuracy and to support +/-0 +/-inf and subnormals, guard and sticky bits/digits. That is not the case for 10 or 12 BCD vintage calculators (they may use one to three guard digits, but not consistently). Also, repeated internal rounding is not beneficial, actually it is "evil". Rounding problems lead to all sorts of numerical problems, hence IEEE754 offers a set of rounding flags to control this. IEEE754 also defines decimal representations that must meet this set of requirements. IEEE754 operations in an implementation must be correctly rounded. When considering implementations, implementing IEEE754 surely takes more effort. I did this in Z80 for Forth850 to support IEEE754 single precision with all requirements including rounding modes (but statically) except for subnormals (maybe sometime...).

To determine accuracy, I don't tend to think in the number of decimal places, but epsilon relative or absolute error. Decimal places do offer a quick insight, but it is useless otherwise. A silly example is 1.499999999 and 1.500000000 which only differ in 1E-9 but share no common digits except the first 1. Now 1.499999999 and 1.49999999 are further apart (9E-9).

One benefit of BCD is for accurate financial calculations perhaps, because amounts such as $0.01 are represented exactly whereas it is not in binary float (it's not 0.01 exactly but 0.010000000000000000208 in double precision). This matters for long running sums and in accumulation loops and when comparing totals. On the other hand, why people claim this but don't scale amounts by a factor 100 to avoid this is beyond me.

- Rob

[Albert Chan]

(06-12-2024 12:45 PM)dm319 Wrote:  I've always thought doing maths stuff in binary is a pain as we generally view the results in decimal,
which makes it hard to determine accuracy in terms of the number of decimal points.

Actually, error analysis is simpler if algorithm use binary math.
Example, decimal with 3 digits precision:

0.999 --> 1 ULP = relative error 0.001 / 0.999 ≈ 0.001
0.100 --> 1 ULP = relative error 0.001 / 0.100 ≈ 0.010

With possible 10 fold difference in relative error, some algorithm may not work at all!
For binary, factor is at most 2 (only in extreme case), and we may ignore this.

Regarding bin ↔ dec ↔ bin conversion errors, if you stick with integers, there is no difference.
For money, instead of using dollars as unit, use cents!

https://www.lua.org/pil/2.3.html Wrote:There is a widespread misconception about floating-point arithmetic errors and some people fear that even a simple increment can go weird with floating-point numbers. The fact is that, when you use a double to represent an integer, there is no rounding error at all (unless the number is greater than 100,000,000,000,000). Specifically, a Lua number can represent any long integer without rounding problems. Moreover, most modern CPUs do floating-point arithmetic as fast as (or even faster than) integer arithmetic.

And, for equivalent precision, binary setup tends to run much, much faster.
This is the reason mpmath choose Python integer type, over Decimal type

https://fredrikj.net/blog/2017/10/mpmath...ospective/ Wrote:When I started sympy.numerics, Python already had the arbitrary-precision Decimal type in the standard library, but I found that it was a lot more efficient to emulate binary floating-point arithmetic using Python's arbitrary-precision long integers. It was especially fast to do fixed-point arithmetic this way: an addition is just an addition, and a multiplication (a * b) >> shift requires just one extra operation, whereas a floating-point operation with rounding requires something like 50 Python bytecode instructions. Implementing core operations (like elementary functions) mostly using fixed-point arithmetic internally made it possible to achieve reasonable performance ...

BTW, Decimal type may "cheat" using binary math!
Example, IntelDecimal128, used in Free42/Plus42

https://www.cl.cam.ac.uk/~jrh13/slides/a...slides.pdf

[dm319]

Very intereting Rob and Albert.

That's a good point about using cents. I only realised yesterday that Excel and LO Calc use double binary rather than decimal, which surprised me because I know Excel is established in finance. It looks like the display value is rounded, so I'm sure most accountants wouldn't notice this.

Also interesting re: the rounding. Is the rounding done beyond the last digit? I.e. if you're rounding to the precision available, that implies you've calculated to beyond that?

I'm still slightly unsure I get the concept of ULP 'accuracy'. Whereas, as you mention Rob, absolute and relative error seems a sensible approach.

[Albert Chan]

(06-22-2024 12:12 PM)dm319 Wrote:  I only realised yesterday that Excel and LO Calc use double binary rather than decimal ...

Do you really know how your spreadsheet works?

Oops, XL did it again

Quote:Also interesting re: the rounding. Is the rounding done beyond the last digit? I.e. if you're rounding to the precision available, that implies you've calculated to beyond that?

Not sure about decimal, but for binary, all you need is 3 extra bits: guard, round, and sticky

[robve]

Updated Newton and Hybrid implementations with SHARP PC-1403H TVM test results.

The improved and combined SHARP BASIC TVM program listing includes both parts for the Hybrid and Newton rate solver methods:

Code:

' 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

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 25
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
25 G=1,C=0 : IF P+F=-N*M LET I=0 : RETURN
' solve Y=NPMT=0 for rate I% pick initial nonzero best guess and perform newton step
26 V=P+B*M,W=F-B*M
27 Y=(V+W)/N,W=W-V-Y,V=(M+Y)/W,W=(N*N-1)*Y*V/W,I=100*V*(W-3)/(W-1.5)
28 IF ABS I<1E-9 LET I=1E-5*(SGN I+(I=0))
29 IF ABS I>99 LET I=99*SGN I
30 GOSUB 38 : IF Y=0 RETURN

' USE THIS PART FOR HYBRID:
31 V=I,I=I-(K*M-(P+F*S)*J/R)*100/(B*M-(P+F)*(1+N*S/(R+R/J))/R-F)
32 GOSUB 38 : IF Y=0 OR Y=W RETURN
' secant loop
33 T=I,I=I-(I-V)*Y/(Y-W),V=T
34 GOSUB 38 : IF Y<>0 IF ABS Y<ABS W GOTO 33
35 IF Y<>0 IF Y<>W IF G LET G=G-1 : GOTO 33

' USE THIS PART FOR NEWTON:
31 V=I,I=I-(K*M-(P+F*S)*J/R)*100/(B*M-(P+F)*(1+N*S/(R+R/J))/R-F)
34 GOSUB 38 : IF Y<>0 IF ABS Y<ABS W GOTO 31
35 IF Y<>0 IF Y<>W IF G LET G=G-1 : GOTO 31

' pick previous I% after overstepping
36 IF Y<>0 LET I=V
' round very small I% to zero
37 I=(ABS I>1E-9)*I : RETURN
' compute Y=NPMT and save previous Y in W
38 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

The Hybrid version executes faster on this machine and its overall accuracy is at least as good as Newton or better (for the set of test cases and other tests that were not included in the results because of similarity).

Add these four lines to compute amortization with <period> DEF-A:

Code:

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

- Rob