Threaded Mode | Linear Mode

SlideRule · 04-04-2022, 01:54 PM

Find attached a two-page extract from Fundamentals of Infrastructure Engineering, Appendix A (pgs. 445 & 447) with respect to Computing the Discount Rate utilizing a generic calculator adaptation of the Newton-Raphson algorithm. The article suggests this listing needs only minor alterations for RPN deployment. We all await the response from the readership to this challenge.

[attachment=10502] [attachment=10503]

BEST!
SlideRule

Thomas Klemm · (This post was last modified: 04-09-2022 08:00 AM by Thomas Klemm.)

There is a bug in the listing.
When calculating

$
- n (i + 1)^{n - 1}
$

the multiplication by $ n $ is missing.

In the third column we find:

Code:

rcl 9           2

aˣ rcl 4 =      3

sto - chs       1

Instead it should be:

Code:

rcl 9           2

aˣ rcl 4 =      3

rcl × 7         2

sto - chs       1

Quote:The machine should have a capacity of 98 keystrokes to perform this routine.

With this fix we are at 100 keystrokes.

Instead of using this suggestion I would rather follow how this is done in the HP-80.

Initial Value

\[i_0 = \frac{2(n-P)}{n(n+1)}\]

Newton Iteration

\[
\begin{align}
Z &= (1 + i)^{-n} \\
\\
i_\text{next} &= i \left[ 1 - \frac{iP - (1 - Z)}{1 - Z - \frac{niZ}{1 + i}} \right]
\end{align}
\]

Register

R0: n
R1: P = PV / PMT
R2: i
R3: u = 1 + i
R4: Z = u^-n
R5: v = Z - 1

Program

This program is for the HP-25, but can easily be adapted for other models:

Code:

 24 00    : RCL 0            ; n

 24 01    : RCL 1            ; P

 41       : -                ;

 02       : 2                ;

 61       : *                ;

 24 00    : RCL 0            ;

 31       : ENTER            ;

 15 02    : g x^2            ;

 51       : +                ;

 71       : /                ; 2 * (n - P) / (n^2 + n)

 23 02    : STO 2            ; i

 01       : 1                ;

 51       : +                ;

 23 03    : STO 3            ; u = 1 + i

 24 00    : RCL 0            ; n

 32       : CHS              ; -n

 14 03    : f y^x            ;

 23 04    : STO 4            ; Z = u^(-n)

 01       : 1                ;

 41       : -                ;

 23 05    : STO 5            ; v = Z - 1

 24 02    : RCL 2            ; i

 24 01    : RCL 1            ; P

 61       : *                ;

 51       : +                ; i * P + v

 24 00    : RCL 0            ; n

 24 02    : RCL 2            ; i

 61       : *                ; n * i

 24 04    : RCL 4            ; Z

 61       : *                ; n * i * Z

 24 03    : RCL 3            ; u

 71       : /                ; n * i * Z / u

 24 05    : RCL 5            ; v

 51       : +                ; v + n * i * Z / u

 71       : /                ; (i * P + v)/(v + n * i * Z / u)

 01       : 1                ;                                     

 51       : +                ; 1 + (i * P + v)/(v + n * i * Z / u)

 24 02    : RCL 2            ; i

 61       : *                ;

 74       : R/S              ; i_next

 13 11    : GTO 11           ;

Quote:The procedure is terminated when two successive values differ by less than a pre-set amount. This test can easily be built-in to the program by a conditional if-then-else command sequence if desired and if sufficient capacity is available. A round-off command should be programmed just prior to the test, of course.

Example: True Interest Rate of a Loan

What is the true interest rate on a 36 month loan for $3000 having monthly payments of $100?

FIX 6

36
STO 0

3000
ENTER
100
/
STO 1

CLEAR PRGM
R/S
0.010190

R/S
0.010207

R/S
0.010207

Thus, after only three iterations, we get the result.
We can now calculate the annual percentage rate:

FIX 2

1200
×
12.25

References

Albert Chan · 04-08-2022, 10:21 PM

Although time flow only from the present to the future, formula does not care.
We can use the *same* code, getting interest rate from PMT, FV (now, with PV=0)
Just have the code run with negative input, traveling back in time Big Grin

Example, lets try n=36, PMT=100, FV=5000, find i

lua> n = -36 -- note the negative sign, for both n and p
lua> p = -5000/100
lua> i = 2*(n-p)/(n*(n+1))
lua> i
0.022222222222222223

Guess i formula is not designed for FV, but lets use it anyway.
With better guess_i(n,pv,pmt,fv), we get guess i = 0.0179687

lua> z = (1+i)^-n
lua> i = i - i*(i*p-(1-z))/(1-z-n*i*z/(1+i))
lua> i
0.018164680329101488
lua> z = (1+i)^-n
lua> i = i - i*(i*p-(1-z))/(1-z-n*i*z/(1+i))
lua> i
0.017958086653757012
lua> z = (1+i)^-n
lua> i = i - i*(i*p-(1-z))/(1-z-n*i*z/(1+i))
lua> i
0.017957584809347397

Albert Chan · (This post was last modified: 04-09-2022 01:24 PM by Albert Chan.)

(10-16-2020 04:02 PM)Albert Chan Wrote: XCas> C := I*N / (1 - (1+I)^-N) // C = |N*PMT/PV|, "compounding factor"
XCas> series(C,I,polynom)

$1
+\frac{I(N+1)}{2}
+\frac{I^2 (N^2-1)}{12}
+\frac{I^3 (-N^2+1)}{24}
+\frac{I^4 (-N^4+20N^2-19)}{720}
+\frac{I^5 (N^4-10N^2+9)}{480}$

This may be a better rate estimate, by dropping compounding factor O(I^2)
With previous defined P, solve for I, we have:

$\displaystyle I ≈ \frac{2\;(N-P)}{P\;(N+1)}$

For previous 2 examples, formula give slightly better guess rate.
Newton steps (denoted by →) converge faster.

I(N= 36, P= 3000/100) ≈ 0.010811 → 0.010203 → 0.010207
I(N=-36, P=-5000/100) ≈ 0.016000 → 0.018003 → 0.017958

Albert Chan · 04-09-2022, 05:47 PM

(02-13-2020 10:42 AM)Gamo Wrote: To find the guess the HP-55 programs book used these formulas.

1. Direct Reduction Loan Interest Rate

[1 / (PV / PMT)] - [(1 / n^2) x (PV / PMT)]

Using our symbols here, we have:

$\displaystyle I ≈ \frac{1}{P} - \frac{P}{N^2}$

This is even better, and work well with big N.

With above guess rate formula, we have:

I(N= 36, P= 3000/100) ≈ 0.010185 → 0.010207
I(N=-36, P=-5000/100) ≈ 0.018580 → 0.017962 → 0.017958

Does anyone know how this is derived ?

This is as good as first pade approximation of reverted C, but less complicated.

$\displaystyle I ≈ \frac{6\;(N-P)}
{N\;(N-1) + 2P\;(N+2)}$

With above guess rate formula, we have:

I(N= 36, P= 3000/100) ≈ 0.010169 → 0.010207
I(N=-36, P=-5000/100) ≈ 0.017751 → 0.017958

Albert Chan · 04-10-2022, 04:03 AM

(04-09-2022 05:47 PM)Albert Chan Wrote: $\displaystyle I ≈ \frac{1}{P} - \frac{P}{N^2}$

This is even better, and work well with big N.
...
Does anyone know how this is derived ?

(10-16-2020 04:02 PM)Albert Chan Wrote: XCas> C := I*N / (1 - (1+I)^-N) // C = |N*PMT/PV|, "compounding factor"

This is my guess on how formula is derived, using continuous compounding.
With small I, (1+I)^-N = exp(ln(1+I)*-N) ≈ exp(-I*N)

Let x = I*N, and solve for x

C = N/P ≈ x / (1 - e^-x)
x = I*N = C + W(-C*exp(-C))            // W has 2 branches, but it does not matter here

Let y = -C*exp(-C), solve for I

I ≈ 1/P + W(y)/N

If we can show -W(y) ≈ 1/C, we have I ≈ 1/P - 1/(CN) = 1/P - P/N^2, and prove is done.

There are 2 LambertW branches, but we only want estimate when C close to 1
Both branches behave about the same, because y is very close to -1/e

C = 1 + ε      → y = 1/e * (-1 + ε^2/2 - ε^3/3 + ε^4/8 - ...)

I borrowed lyuka e^W estimation formula, replace 0.3 by 1/3, for Puiseux series.

e^W(y) ≈ 1/e + sqrt ((2/e)*(y+1/e)) + (y+1/e)/3

-W(y)
= 1 - ln(e * e^W(y))
= 1 - ln(1 + sqrt(2*(e*y+1)) + (e*y+1)/3)
= 1 - ln(1 + ε - 1/6*ε^2 + ...)
= 1 - 2*atanh(ε/2 - 1/3*ε^2 + ...)
≈ 1 - ε + 2/3*ε^2
≈ 1/(1+ε)
= 1/C            QED

To confirm math is correct, at least numerically.

lua> z = 1e-4
lua> c = 1 + z
lua> y = -c*exp(-c)
lua> -W(y), 1 - z + 2/3*z^2, 1/c
0.9999000066665401      0.9999000066666667      0.9999000099990001

Thomas Klemm · (This post was last modified: 04-12-2022 11:16 PM by Thomas Klemm.)

After a closer inspection of the assembler code I noticed that in fact the following formula is used:

Newton Iteration

\[
\begin{align}
Z &= (1 + i)^{-n} \\
\\
i_\text{next} &= i \left[ 1 + \frac{P + \frac{Z - 1}{i}}{\frac{Z - 1}{i} + \frac{nZ}{1 + i}} \right]
\end{align}
\]

This saves two multiplications by $ i $ at the cost of one division.

Register

R0: n
R1: P = PV / PMT
R2: i
R3: u = 1 + i
R4: Z = u^-n
R5: v = (Z - 1) / i

Program

It allows to reduce the program by 2 steps:

Code:

 24 00    : RCL 0

 24 01    : RCL 1

 41       : -

 02       : 2

 61       : *

 24 00    : RCL 0

 31       : ENTER

 15 02    : g x^2

 51       : +

 71       : /

 23 02    : STO 2

 01       : 1

 51       : +

 23 03    : STO 3

 24 00    : RCL 0

 32       : CHS

 14 03    : f y^x

 23 04    : STO 4

 01       : 1

 41       : -

 24 02    : RCL 2

 71       : /

 23 05    : STO 5

 24 01    : RCL 1

 51       : +

 24 00    : RCL 0

 24 04    : RCL 4

 61       : *

 24 03    : RCL 3

 71       : /

 24 05    : RCL 5

 51       : +

 71       : /

 01       : 1

 51       : +

 24 02    : RCL 2

 61       : *

 74       : R/S

 13 11    : GTO 11

Eddie W. Shore · (This post was last modified: 04-10-2022 11:29 PM by Eddie W. Shore.)

HP 32S/HP 32SII code:

N = n (R0)
P = PV/PMT (R1)
I = i (R2)
U = 1+i (R3)
Z = u^(-n) (R4)
V = Z - 1 (R5)

1. Press XEQ I, enter N, enter PV
N = n; P = PV/PMT (or)
N = -n; P = -FV/PMT
2. Press XEQ J
3. Keep pressing R/S until desired accuracy is reached

Code:

I01 LBL I

I02 INPUT N

I03 INPUT P

I04 - 

I05 2

I06 ×

I07 RCL N

I08 ENTER

I09 x^2

I10 +

I11 ÷

I12 STOP

J01 LBL J

J02 STO I

J03 1

J04 +

J05 STO U

J06 RCL N

J07 +/-

J08 y^x

J09 STO Z

J10 1

J11 -

J12 STO V

J13 RCL I

J14 RCL× P

J15 +

J16 RCL N

J17 RCL× I

J18 RCL× Z

J19 RCL÷ U

J20 RCL+V

J21 ÷

J22 

J23 +

J24 RCL× I

J25 STOP

J26 GTO J

Size:
I: 18.0 bytes, Checksum 32SII: 82BD
J: 39.0 bytes, Checksum 32SII: 3EBB
57 bytes total

Albert Chan · 04-11-2022, 01:18 AM

(04-10-2022 12:19 PM)Thomas Klemm Wrote: Newton Iteration
\[
\begin{align}
Z &= (1 + i)^{-n} \\
\\
i_\text{next} &= i \left[ 1 + \frac{P + \frac{Z - 1}{i}}{\frac{Z - 1}{i} + \frac{nZ}{1 + i}} \right]
\end{align}
\]

When interest rate converging to true rate, numerator goes 0, or (Z-1)/i = -P

I was reading https://brownmath.com/bsci/loan.htm#Newton.
It's rate newton formula actually replaced denominator (Z-1)/i by -P

The formula quoted above is good when guess rate is under-estimated.
This explained why it has i0 = 2*(n-p)/(n*(n+1)) instead of 2*(n-p)/(p*(n+1))

Too high a guess can make Newton's step diverge.
Example, with N=36, P=30, guess rate of 0.08, it barely able to converge.

0.08 → -0.071887 → -0.041499 → -0.014129 → 0.003911 → 0.009742 → 0.010205 → 0.010207

"-P" version had the opposite problem, and preferred over-estimated guess.
Again, with N=36, P=30, guess rate of 0.004, it totally diverged

0.004 → -0.009120 → -0.003414 → -0.000780 → -0.000059 ...

A better option is to solve NPMT=0 instead (now used for rate search in Plus42, see tvm_eq())
From NPMT=0 based loan_rate(n, pv, pmt, fv, i), simplified with FV=0, we get this:

$ \displaystyle \large
i_\text{next} = i \left[
1 -
\frac{(\frac{Z - 1}{i})/P + 1}
{(\frac{nZ}{1 + i})/(\frac{Z - 1}{i}) +1} \right] $

This version, convergence is slightly faster, and *very* stable for guess i
Again, with N=36, P=30, but ridiculous rate guess:

10⁺¹⁰ → 0.033333 → 0.012232 → 0.010228 → 0.010207
10⁻¹⁰ → 0.000231 → 0.010782 → 0.010209 → 0.010207

Thomas Klemm · 04-11-2022, 02:37 AM

(04-10-2022 11:29 PM)Eddie W. Shore Wrote:

Code:

J20 RCL+ V J21 ÷ J22 J23 + J24 RCL× I

There's a 1 missing in line J22.

Using my 2nd program could possibly save a byte or two.
You could also consider using the solver.
This made me wonder if there's a programmable HP calculator with recall arithmetic but without a solver (generic or specific for the TVM problem).

Albert Chan · (This post was last modified: 04-12-2022 03:53 PM by Albert Chan.)

(04-11-2022 01:18 AM)Albert Chan Wrote: The formula [OP] is good when guess rate is under-estimated.
This explained why it has i0 = 2*(n-p)/(n*(n+1)) instead of 2*(n-p)/(p*(n+1))

To be more precise, OP formula preferred rate magnitude under-estimated.

Plotting errors of both rate estimate, error(n version) ≈ -2 * error(p version)
In other words, this is much better rate estimate.

$\displaystyle I_0 = \frac{2\;(N-P)}{N+1}
\left( \frac{{1 \over N} + {2 \over P}}{3} \right)
$

Examples:

I₀(N= 36, P= 30) = .010210 // true rate = .010207
I₀(N=-36, P=-50) = .018074 // true rate = .017958

Nice. But, for OP formula, we wanted guess rate under-estimated.
Perhaps give up some accuracy for safety in convergence, and make ratio 1:1 ?
While we are at it, why not assume N+1 ≈ N, to simplify further ?

$\displaystyle I_0 = \frac{2\;(N-P)}{N}
\left( \frac{{1 \over N} + {1 \over P}}{2} \right)
= \frac{1}{P} - \frac{P}{N^2}
$

This may be how I₀ = 1/P - P/N² comes from Smile

Update: Perhaps formula designed also to match asymptote, when C is big.
When compounding factor C is big, literally all payments goes to paying interest.

C = N/P = I*N/(1-(1+I)^-N) ≈ N*I → I = 1/P

Rate formula matched this behavior: I₀ = 1/P - P/N^2 = (1 - 1/C^2)/P ≈ 1/P

For better estimate, we can use continuous compouding for I, derived earlier.

(04-10-2022 04:03 AM)Albert Chan Wrote: Let y = -C*exp(-C), solve for I

I ≈ 1/P + W(y)/N

When C = N/P is big, y is small, we have W(y) ≈ y

I ≈ C/N + (-C*exp(-C))/N = (1-exp(-N/P)) / P

Albert Chan · 04-11-2022, 03:57 PM

(04-11-2022 03:11 PM)Albert Chan Wrote: Plotting errors of both rate estimate, error(n version) ≈ -2 * error(p version)

Instead of looking at plots, lets look at reverted series of I, from C = N/P
More specifically, from u = (C-1)/(N+1)

From A Series of Interest, by David W. Cantrell

I = (2*u)*(1 - (N-1)/3*u * (1 - (2*N+1)/3*u * (1 - (11*N+7)/15*u * (1 - ...

Assume 2 terms of I is good enough for estimating true I
Assume N big enough, so that N+1 ≈ N-1 ≈ N
Let C = 1 + ε

I/(2u) ≈ 1 − (N/3) * (ε/N) = 1 − 1/3*ε
I/(2u/C) ≈ (1+ε) * (1-ε/3) ≈ 1 + 2/3*ε

Rate guess, p version is (2u), n version is (2u/C), rel error ratio ≈ -2 QED

Thomas Klemm · 04-11-2022, 08:48 PM

Comparing the code in HP Journal to the US patent, I noticed something slightly interesting:
In the US patent DOWN ROTATE is used twice in several places, while in HP Journal a JSB DR2 is used.

Here are the corresponding code snippets for comparison:

US Patent

Code:

L3222:   .....1...1      → L3004               JSB XTY

L3223:   11..1.1...                            DOWN ROTATE

L3224:   .1..1.1...                            C → STACK

L3225:   111..1.1.1      → L3345               JSB MPY

L3226:   ..1...111.                            B → C[W]

L3227:   11..1.1...                            DOWN ROTATE

L3230:   11..1.1...                            DOWN ROTATE

L3231:   111...11.1      → L3343               JSB DIV

L3232:   11..1.1...                            DOWN ROTATE

L3233:   11..1.1...                            DOWN ROTATE

L3234:   .1...11..1      → L3106               JSB ONE

L3235:   .11.1.1..1      → L3152               JSB SUB

L3236:   1.1.1.1...                            M → C[W]

L3237:   111...11.1      → L3343               JSB DIV

L3240:   11..1.1...                            DOWN ROTATE

L3241:   11..1.1...                            DOWN ROTATE

L3242:   .11.1.1..1      → L3152               JSB SUB

HP Journal

Code:

L3224:   .....1...1      → L3004               JSB XTY

L3225:   11..1.1...                            DOWN ROTATE

L3226:   .1..1.1...                            C → STACK

L3227:   111..1...1      → L3344               JSB MPY

L3230:   ..1...111.                            B → C[W]

L3231:   ...1.1.1.1      → L3025               JSB DR2

L3232:   111...11.1      → L3343               JSB DIV

L3233:   ...1.1.1.1      → L3025               JSB DR2

L3234:   .1...11..1      → L3106               JSB ONE

L3235:   .11.1.1..1      → L3152               JSB SUB

L3236:   1.1.1.1...                            M → C[W]

L3237:   111...11.1      → L3343               JSB DIV

L3240:   ...1.1.1.1      → L3025               JSB DR2

L3241:   .11.1.1..1      → L3152               JSB SUB

While I couldn't find the DR2 label in ROM 3 of the US patent, I was able to find similar code in ROM 5:

Code:

L5071:   11..1.1...                    DOWN3 : DOWN ROTATE

L5072:   11..1.1...                    DOWN2 : DOWN ROTATE

L5073:   11..1.1...                            DOWN ROTATE

L5074:   ....11....                            RETURN

Both DOWN2 and DOWN3 are used in this ROM.

I wasn't aware of several ROM versions of the HP-80.

Based on Bernhard's statement I assume that the code from the US patent is used in the HP-80:

Quote:I compared my HP-80 ROM code and found, that is was correct.

While some of the ROMs can be downloaded, the HP-80 unfortunately is missing.
Therefore I can not verify that.

If anyone can shed some light on the subject I would be very grateful.

Thomas Klemm · 04-12-2022, 11:11 PM

(04-10-2022 12:19 PM)Thomas Klemm Wrote: After a closer inspection of the assembler code I noticed that in fact the following formula is used:

For this I wrote a simulator for the CPU, similar to the Code Analyzer for HP Calculators:

Code:

@trace

def ADD(): # 03153

    global A, B, C, D, E, F, M

    B = A

    A = C = A + C

@trace

def SUB(): # 03152

    global A, B, C, D, E, F, M

    B = A

    A = C = A - C

@trace

def MPY(): # 03345

    global A, B, C, D, E, F, M

    B = A

    A = C = A * C

@trace

def DIV(): # 03343

    global A, B, C, D, E, F, M

    B = C

    A = C = A / C

@trace

def XTY(): # 03004

    global A, B, C, D, E, F, M

    B = A

    D = C

    A = C = A**C

@trace

def ONE(): # 03106

    global A, B, C, D, E, F, M

    A, C = C, 1

@trace

def CSN(): # 03321

    global A, B, C, D, E, F, M

    E = - E

@trace

def DOWN_ROTATE():

    global A, B, C, D, E, F, M

    C, D, E, F = D, E, F, C

@trace

def STACK_A():

    global A, B, C, D, E, F, M

    A, D, E = D, E, F

@trace

def C_STACK():

    global A, B, C, D, E, F, M

    D, E, F = C, D, E

@trace

def A_EXCHANGE_B():

    global A, B, C, D, E, F, M

    A, B = B, A

@trace

def B_EXCHANGE_C():

    global A, B, C, D, E, F, M

    B, C = C, B

@trace

def C_EXCHANGE_M():

    global A, B, C, D, E, F, M

    C, M = M, C

@trace

def B_C():

    global A, B, C, D, E, F, M

    C = B

@trace

def C_A():

    global A, B, C, D, E, F, M

    A = C

@trace

def M_C():

    global A, B, C, D, E, F, M

    C = M

@trace

def TEN6(): # 03354

    global A, B, C, D, E, F, M

    C = M

    C *= 100

    A *= 10000

def trace(operation):

    def fmt(s):

        return str(s)

    def wrapper():

        operation()

        if SYMBOLIC:

            print((

              f"A: {fmt(A):<45s}    "

              f"B: {fmt(B):<45s}    "

              f"C: {fmt(C):<45s}    "

              f"D: {fmt(D):<45s}    "

              f"E: {fmt(E):<45s}    "

              f"F: {fmt(F):<45s}    "

              f"M: {fmt(M):<45s}    "

              f": {operation.__name__}  "

            ))

        else:

            print((

              f"A: {A:< 20.12e}    "

              f"B: {B:< 20.12e}    "

              f"C: {C:< 20.12e}    "

              f"D: {D:< 20.12e}    "

              f"E: {E:< 20.12e}    "

              f"F: {F:< 20.12e}    "

              f"M: {M:< 20.12e}    "

              f": {operation.__name__}  "

            ))

    return wrapper

First we create some symbols:

Code:

from sympy import symbols, simplify, expand, latex

n, PMT, PV, P, i, Z, Δi = symbols("n PMT PV P i Z Δi")

Then we reset the registers:

Code:

A, B, C, D, E, F, M = [0] * 7

SYMBOLIC = True

And now we just follow the assembler code.
From time to time, we reassign variables to keep the expressions simple.

PS: I had to split this post since it was too big.

Thomas Klemm · (This post was last modified: 04-13-2022 01:26 AM by Thomas Klemm.)

Code:

C, D, E = PV, PMT, n

C_STACK()

STACK_A()

DOWN_ROTATE()

DIV()

A: 0 B: 0 C: PV D: PV E: PMT F: n M: 0 : C_STACK
A: PV B: 0 C: PV D: PMT E: n F: n M: 0 : STACK_A
A: PV B: 0 C: PMT D: n E: n F: PV M: 0 : DOWN_ROTATE
A: PV/PMT B: PMT C: PV/PMT D: n E: n F: PV M: 0 : DIV

Code:

print(latex(C))

\frac{PV}{PMT}

\[
P = \frac{PV}{PMT}
\]

Code:

A = C = P

C_EXCHANGE_M()

DOWN_ROTATE()

C_STACK()

ONE()

ADD()

A_EXCHANGE_B()

MPY()

C_EXCHANGE_M()

STACK_A()

C_STACK()

CSN()

DOWN_ROTATE()

SUB()

ADD()

C_EXCHANGE_M()

DIV()

A: P B: PMT C: 0 D: n E: n F: PV M: P : C_EXCHANGE_M
A: P B: PMT C: n D: n E: PV F: 0 M: P : DOWN_ROTATE
A: P B: PMT C: n D: n E: n F: PV M: P : C_STACK
A: n B: PMT C: 1 D: n E: n F: PV M: P : ONE
A: n + 1 B: n C: n + 1 D: n E: n F: PV M: P : ADD
A: n B: n + 1 C: n + 1 D: n E: n F: PV M: P : A_EXCHANGE_B
A: n*(n + 1) B: n C: n*(n + 1) D: n E: n F: PV M: P : MPY
A: n*(n + 1) B: n C: P D: n E: n F: PV M: n*(n + 1) : C_EXCHANGE_M
A: n B: n C: P D: n E: PV F: PV M: n*(n + 1) : STACK_A
A: n B: n C: P D: P E: n F: PV M: n*(n + 1) : C_STACK
A: n B: n C: P D: P E: -n F: PV M: n*(n + 1) : CSN
A: n B: n C: P D: -n E: PV F: P M: n*(n + 1) : DOWN_ROTATE
A: -P + n B: n C: -P + n D: -n E: PV F: P M: n*(n + 1) : SUB
A: -2*P + 2*n B: -P + n C: -2*P + 2*n D: -n E: PV F: P M: n*(n + 1) : ADD
A: -2*P + 2*n B: -P + n C: n*(n + 1) D: -n E: PV F: P M: -2*P + 2*n : C_EXCHANGE_M
A: (-2*P + 2*n)/(n*(n + 1)) B: n*(n + 1) C: (-2*P + 2*n)/(n*(n + 1)) D: -n E: PV F: P M: -2*P + 2*n : DIV

Code:

print(latex(simplify(C)))

\frac{2 \left(- P + n\right)}{n \left(n + 1\right)}

\[
i = \frac{2 \left(- P + n\right)}{n \left(n + 1\right)}
\]

Code:

A = C = i

C_EXCHANGE_M()

DOWN_ROTATE()

C_EXCHANGE_M()

STACK_A()

ONE()

ADD()

C_STACK()

C_STACK()

B_C()

C_EXCHANGE_M()

XTY()

A: i B: n*(n + 1) C: -2*P + 2*n D: -n E: PV F: P M: i : C_EXCHANGE_M
A: i B: n*(n + 1) C: -n D: PV E: P F: -2*P + 2*n M: i : DOWN_ROTATE
A: i B: n*(n + 1) C: i D: PV E: P F: -2*P + 2*n M: -n : C_EXCHANGE_M
A: PV B: n*(n + 1) C: i D: P E: -2*P + 2*n F: -2*P + 2*n M: -n : STACK_A
A: i B: n*(n + 1) C: 1 D: P E: -2*P + 2*n F: -2*P + 2*n M: -n : ONE
A: i + 1 B: i C: i + 1 D: P E: -2*P + 2*n F: -2*P + 2*n M: -n : ADD
A: i + 1 B: i C: i + 1 D: i + 1 E: P F: -2*P + 2*n M: -n : C_STACK
A: i + 1 B: i C: i + 1 D: i + 1 E: i + 1 F: P M: -n : C_STACK
A: i + 1 B: i C: i D: i + 1 E: i + 1 F: P M: -n : B_C
A: i + 1 B: i C: -n D: i + 1 E: i + 1 F: P M: i : C_EXCHANGE_M
A: (i + 1)**(-n) B: i + 1 C: (i + 1)**(-n) D: -n E: i + 1 F: P M: i : XTY

Code:

print(latex(C))

\left(i + 1\right)^{- n}

\[
Z = \left(i + 1\right)^{- n}
\]

Code:

A = C = Z

DOWN_ROTATE()

C_STACK()

MPY()

B_C()

DOWN_ROTATE()

DOWN_ROTATE()

DIV()

DOWN_ROTATE()

DOWN_ROTATE()

ONE()

SUB()

M_C()

DIV()

DOWN_ROTATE()

DOWN_ROTATE()

SUB()

STACK_A()

C_STACK()

B_EXCHANGE_C()

ADD()

DOWN_ROTATE()

B_EXCHANGE_C()

DOWN_ROTATE()

B_EXCHANGE_C()

DIV()

M_C()

MPY()

C_EXCHANGE_M()

A: Z B: i + 1 C: -n D: i + 1 E: P F: Z M: i : DOWN_ROTATE
A: Z B: i + 1 C: -n D: -n E: i + 1 F: P M: i : C_STACK
A: -Z*n B: Z C: -Z*n D: -n E: i + 1 F: P M: i : MPY
A: -Z*n B: Z C: Z D: -n E: i + 1 F: P M: i : B_C
A: -Z*n B: Z C: -n D: i + 1 E: P F: Z M: i : DOWN_ROTATE
A: -Z*n B: Z C: i + 1 D: P E: Z F: -n M: i : DOWN_ROTATE
A: -Z*n/(i + 1) B: i + 1 C: -Z*n/(i + 1) D: P E: Z F: -n M: i : DIV
A: -Z*n/(i + 1) B: i + 1 C: P D: Z E: -n F: -Z*n/(i + 1) M: i : DOWN_ROTATE
A: -Z*n/(i + 1) B: i + 1 C: Z D: -n E: -Z*n/(i + 1) F: P M: i : DOWN_ROTATE
A: Z B: i + 1 C: 1 D: -n E: -Z*n/(i + 1) F: P M: i : ONE
A: Z - 1 B: Z C: Z - 1 D: -n E: -Z*n/(i + 1) F: P M: i : SUB
A: Z - 1 B: Z C: i D: -n E: -Z*n/(i + 1) F: P M: i : M_C
A: (Z - 1)/i B: i C: (Z - 1)/i D: -n E: -Z*n/(i + 1) F: P M: i : DIV
A: (Z - 1)/i B: i C: -n D: -Z*n/(i + 1) E: P F: (Z - 1)/i M: i : DOWN_ROTATE
A: (Z - 1)/i B: i C: -Z*n/(i + 1) D: P E: (Z - 1)/i F: -n M: i : DOWN_ROTATE
A: Z*n/(i + 1) + (Z - 1)/i B: (Z - 1)/i C: Z*n/(i + 1) + (Z - 1)/i D: P E: (Z - 1)/i F: -n M: i : SUB
A: P B: (Z - 1)/i C: Z*n/(i + 1) + (Z - 1)/i D: (Z - 1)/i E: -n F: -n M: i : STACK_A
A: P B: (Z - 1)/i C: Z*n/(i + 1) + (Z - 1)/i D: Z*n/(i + 1) + (Z - 1)/i E: (Z - 1)/i F: -n M: i : C_STACK
A: P B: Z*n/(i + 1) + (Z - 1)/i C: (Z - 1)/i D: Z*n/(i + 1) + (Z - 1)/i E: (Z - 1)/i F: -n M: i : B_EXCHANGE_C
A: P + (Z - 1)/i B: P C: P + (Z - 1)/i D: Z*n/(i + 1) + (Z - 1)/i E: (Z - 1)/i F: -n M: i : ADD
A: P + (Z - 1)/i B: P C: Z*n/(i + 1) + (Z - 1)/i D: (Z - 1)/i E: -n F: P + (Z - 1)/i M: i : DOWN_ROTATE
A: P + (Z - 1)/i B: Z*n/(i + 1) + (Z - 1)/i C: P D: (Z - 1)/i E: -n F: P + (Z - 1)/i M: i : B_EXCHANGE_C
A: P + (Z - 1)/i B: Z*n/(i + 1) + (Z - 1)/i C: (Z - 1)/i D: -n E: P + (Z - 1)/i F: P M: i : DOWN_ROTATE
A: P + (Z - 1)/i B: (Z - 1)/i C: Z*n/(i + 1) + (Z - 1)/i D: -n E: P + (Z - 1)/i F: P M: i : B_EXCHANGE_C
A: (P + (Z - 1)/i)/(Z*n/(i + 1) + (Z - 1)/i) B: Z*n/(i + 1) + (Z - 1)/i C: (P + (Z - 1)/i)/(Z*n/(i + 1) + (Z - 1)/i) D: -n E: P + (Z - 1)/i F: P M: i : DIV
A: (P + (Z - 1)/i)/(Z*n/(i + 1) + (Z - 1)/i) B: Z*n/(i + 1) + (Z - 1)/i C: i D: -n E: P + (Z - 1)/i F: P M: i : M_C
A: i*(P + (Z - 1)/i)/(Z*n/(i + 1) + (Z - 1)/i) B: (P + (Z - 1)/i)/(Z*n/(i + 1) + (Z - 1)/i) C: i*(P + (Z - 1)/i)/(Z*n/(i + 1) + (Z - 1)/i) D: -n E: P + (Z - 1)/i F: P M: i : MPY
A: i*(P + (Z - 1)/i)/(Z*n/(i + 1) + (Z - 1)/i) B: (P + (Z - 1)/i)/(Z*n/(i + 1) + (Z - 1)/i) C: i D: -n E: P + (Z - 1)/i F: P M: i*(P + (Z - 1)/i)/(Z*n/(i + 1) + (Z - 1)/i) : C_EXCHANGE_M

Code:

print(latex((M)))

\frac{i \left(P + \frac{Z - 1}{i}\right)}{\frac{Z n}{i + 1} + \frac{Z - 1}{i}}

\[
\Delta i = \frac{i \left(P + \frac{Z - 1}{i}\right)}{\frac{Z n}{i + 1} + \frac{Z - 1}{i}}
\]

Code:

A = M = Δi

ADD()

C_EXCHANGE_M()

A: i + Δi B: Δi C: i + Δi D: -n E: P + (Z - 1)/i F: P M: Δi : ADD
A: i + Δi B: Δi C: Δi D: -n E: P + (Z - 1)/i F: P M: i + Δi : C_EXCHANGE_M

Code:

print(latex((M)))

i + Δi

\[
i_\text{next} = i + \Delta i
\]

Thomas Klemm · 04-18-2022, 01:58 PM

(04-11-2022 08:48 PM)Thomas Klemm Wrote: If anyone can shed some light on the subject I would be very grateful.

Meanwhile I found the following comment in asm/80.asm of nonpareil:

Quote:; The code from the patent may not match released code in actual HP-80
; calculators. It does not match the interest rate calculation code
; listed in Figure 3 of the article "A Pocket-Sized Answer Machine for
; Business and Finance" in the May 1973 issue of the Hewlett-Packard
; Journal. The patent application was filed on October 7, 1972, so
; the code in the Journal code might be more recent.

Obviously I'm not the first to notice this discrepancy.

**rprosperi** · 04-18-2022, 06:41 PM

On (far too) many occasions, when I think I may be the first to discover something, upon checking, I have found that Eric had found it, documented it, shared it, and just as likely forgotten it, long before I got there.

Same with Joe Horn...

<sigh>...

Thomas Klemm · 04-18-2022, 07:03 PM

Well, it's not the first time I've taken a close look at the assembly code from this patent.

Interestingly the bug is still there in asm/80.asm:

Code:

dd2:    shift right c[w]

    p - 1 -> p

    if p # 14

         then go to dd2

    if c[x] = 0

         then go to err71

Albert Chan · 05-11-2022, 06:22 PM

(04-11-2022 03:11 PM)Albert Chan Wrote: Plotting errors of both rate estimate, error(n version) ≈ -2 * error(p version)
In other words, this is much better rate estimate.

$\displaystyle I_0 = \frac{2\;(N-P)}{N+1}
\left( \frac{{1 \over N} + {2 \over P}}{3} \right)
$

Examples:

I₀(N= 36, P= 30) = .010210 // true rate = .010207
I₀(N=-36, P=-50) = .018074 // true rate = .017958

There are better ways to compensate the 1:2 error ratio.

I recently updated guess_i(n, pv, pmt, fv) that keep C series terms upto I^2
Converting it with FV=0, I1 = guess_i(N, P, -1, 0):

Rough estimate: $\displaystyle I_0 = \frac{2\;(N-P)}{P(N+1)} $

With correction: $\displaystyle I_1
= I_0 \left( \frac{6 + 1\;(n-1)\;I_0}{6 + 2\;(n-1)\;I_0} \right)
= I_0 - \frac{I_0}{2 + \large \frac{6}{(N-1)\;I_0}} $

Redo examples:

N= 36, P= 30 → I0 = .010811 → I1 = .010205
N=-36, P=-50 → I0 = .016000 → I1 = .017967