31 December 2018

Alpha stack on the HP-41C and HP-42S

A while back, another member of the MoHPC forum mentioned a program he'd written for the HP-41C that allows the user to manage multiple alpha registers. I responded saying that I had written something similar many moons ago that behaved like a LIFO (last-in-first-out) stack rather than an indexed array of datasets and said that I'd look it up.

I have no idea whatsoever what I did with my little utility so I decided to rewrite it, purely and simply. So here it is.

NB: This program creates and manages a data file in Extended Memory called "ASTACK". If you already have a file of that name, it will be deleted! Note also that running these programs uses Flag 01 and trashes registers R07-R10. Finally, since this uses Extended Memory, you will need to run this on a 41CX or SwissMicros DM41, or on a 41C or 41CV with the "X-Function" module.

The size of the data file created in Extended Memory depends on the depth of the alpha stack that you want to create. Two registers are needed for a header in the file and four registers per stack level are needed. So, if you want a stack that's 6 levels deep, for example, then you'll need room in your Extended Memory for a file that's 2+6*4= 26 registers in size. You'll actually need 28 registers free because the calculator also steals 2 registers for its own internal housekeeping when you create a file in X-Mem.

The three utilities provided are "ASINIT", "APUSH" and "APOP".

Run this with the depth of the desired stack in X, The '41's own error detection will prevent you from creating a file that's too big or from running this on a machine with no "X-Function" module installed (remember, the 41CX and the DM41 have this module baked into their ROM).

This will save the current contents of your alpha register onto the alpha stack and return 0 in X, unless the stack is already full, in which case you'll get -3 back instead. If the alpha stack hasn't been initialized (by running ASINIT) then you'll get -1 back in X.

This takes the string on the top of the alpha stack and transfers it into the '41's alpha register, removing it from the stack. If all went well, X will contain 0 after returning from this program. If the alpha stack has not yet been initialized then you'll get -1 back, or if the stack was already empty (everything already popped off it) when you called APOP then you'll get -2 back.

You can go and grab these utilities here: alpha-stack.zip

Software provided, as usual, as a text listing, a .raw file and a PDF with bar codes.

Edit: I wrote something similar for HP-42S/DM42. This machine has no Extended Memory but can store register data into a matrix. That's exactly what we do with this version.

Grab the listing here: astack.txt and the .raw file for DM42/Free42 here: astack.raw

Listing reproduced here:
00 { 196-Byte Prgm }
02 ABS
03 IP
05 8
06 ×
07 2
08 +
09 1
12 CLX
13 R↓
16 RTN
17▸LBL 10
18 CLX
19 SF 25
21 FC?C 25
22 -1
23 RTN
25 XEQ 10
26 X≠0?
27 RTN
29 I+
31 X≠0?
32 GTO 01
33 -2
34 RTN
35▸LBL 01
36 8
37 ×
38 5
39 -
40 1
42 CLA
43 8
44▸LBL 11
46 I+
48 R↓
50 GTO 11
51 FS?C 01
52 GTO 04
53 2
54 1
57 1
58 -
60▸LBL 04
62 RTN
64 XEQ 10
65 X≠0?
66 RTN
68 I+
70 X=Y?
71 GTO 02
72 1
73 +
75 8
76 ×
77 5
78 -
79 1
81 8
82 R↑
83▸LBL 09
86 I+
89 GTO 09
90 SF 01
92 END

The "N Queens" problem

The "N Queens" problem is a puzzle where you take an NxN chequers or chess board (it's usually 8x8 but this works for the general case) on which you have to place N queens in such a way that none of them threaten or are threatened by any of the others.

Reminder: in the game of chess, a queen can move as many squares as she likes horizontally, vertically or diagonally but cannot jump over any of the other pieces on the board like a knight. This means that any piece in her direct line of sight on the same line, in the same column or on a diagonal to her is fair game.

With this in mind, it is clear that there is only one solution to the problem on a 1x1 board, i.e. a single square (you just plonk a queen on it and you're done), and that there is no way to organise two queens on a 2x2 board or three queens on a 3x3 board without them threatening each other. For a 4x4 or greater board, there is more than one solution.

I thought I'd have a go at solving the general case on an HP-71B and the program listed in its entirety below does just that. There's also a brief video of the 71B in action finding all the solutions to the problem for a 5x5 board in just over 30 seconds.

Here is a line-by-line explanation of the algorithm:
Clear up any variables currently in memory and tell the 71B not to waste time pausing after displaying something on screen. This is going to be slow enough as it is :)
I want the elements in the array that I'm creating later to be numbered from 1 onwards rather than 0, and we need integers to be displayed with no decimals or we're going to get a string overflow (and something far too long to display on screen) in line 200.
OK, let's set up a few variables.

R and C are going to represent the row and column of the square on the board where I want to see if I can place a queen.

S is the size of the board.

T will be used as a counter while I'm seeing if a queen that I place on the board is being threatened by one placed earlier on a row with a lower number.

F is the running count of solutions found so far.

N is the number of nodes of the puzzle examined so far. Each time a queen is tentatively placed on the board, this is a puzzle node.

K will be used as a counter while building the string representing the solution that was found.

W, the only real variable used here, holds the time of day when the run was started for most of the life cycle of this program. At the end it is replaced by the time, in seconds, that it took to find all the solutions to the puzzle.
50 INPUT "Board size 1-9: ";S
60 IF S>9 OR S<1 THEN 50
Ask the user what size board to work with. Anything below 1 makes no sense and anything above 9 will just take far too long to solve. It already takes about 50 minutes to solve for an 8x8 board...
70 INTEGER B(S) @ DISP "Working..."
80 R=0 @ F=0 @ N=0 @ W=TIME
We start by creating an array of integers B() with as many elements as there are rows on the board. This array will hold the number of the column where there is a queen on each row. We then display something on screen to let the user know that we've started crunching numbers. Depending on the size of the board we're working with, it can take a while to find the first solution.

We then initialize the row counter and the running counts and we make a note of the time when we started.
90 R=R+1 @ IF R>S THEN 190
100 C=0
This is the point in the algorithm where we move on to the next row. If we've only just started then the "next" row is actually the first row. Otherwise, we come here when we've successfully placed a queen on the board and we want to try and place one on the next row.

However, if we were already on the last row of the board and we try and count beyond that number, then it means that we just placed a queen on the last row, i.e. we found a solution! In this case we go off to line 190, where this situation is handled. If not, then we reset the column counter so that we can look at all of the squares in this row one after the other.
110 C=C+1 @ IF C>S THEN 170
So, we move on to the next column. If we've only just started looking at the current row then this is the first column. Otherwise it means that we couldn't place a queen in the previous column or that we could, but we're looking for solutions that involve placing her on a different square in the same row.

However, if there are no more columns, i.e. if we attempt to count beyond the size of the board, then it means that we need to backtrack a row and start looking for solutions with a queen in a different column of the previous row. This is handled in line 170.
120 T=0
130 N=N+1
We're still on the board so now we need to see if a queen placed in the square identified by R and C is safe. We're going to use the variable T as a counter iterating through all the rows previous to R to see if a queen in any of those rows would threaten one here. We also increment the variable N to indicate that we've found a new node to examine.
140 T=T+1 @ IF T>=R THEN 160
We're going to look at the next row now, but if the next row happens to be this row (R), then it means we've found no threat coming from previous rows and a queen placed on the square that we're looking at right now is safe. This is handled in line 160.
150 IF B(T)=C OR ABS(B(T)-C)=R-T THEN 110 ELSE 140
This is where we determine if the queen in row T is a threat to us in row R, column C.

Remember, the array B() holds the column numbers of the queens in each row. Thus, B(1) holds the column of the queen in row 1, B(2) that of the queen in row 2 etc. B(T) is the column of the queen in row T. If the queen in that row is in this column (if B(T)=C) then she would be a threat to us. Similarly, if she is removed from us by as many columns as rows (if the absolute value of B(T)-C is the same as R-T) then that queen is on the same diagonal as us and would be a threat. In either case, we can't place a queen on this square, so we go back to line 110 where we look at the next column.

If, however, the queen in line T would not be a threat to us then we can look at the next line by jumping back to line 140.
160 B(R)=C @ GOTO 90
If we make it here then it means that we can place a queen on row R, column C. How do we do this? By setting the relevant element of the B() array, B(R), to the column where we're placing the queen. We then go back to line 90 so that we can move on to the next row.
170 R=R-1 @ IF R=0 THEN 230
180 C=B(R) @ GOTO 110
We came here from line 110 when we tried to look at a column beyond the last one in the current row. This meant that we couldn't place a queen anywhere (else) in this row, so we have to backtrack a row and continue the search. So we start by decrementing the current row counter.

However, if we were on the first row anyway and we can't place a queen in that row now, then we have reached the end of the puzzle. There are no more solutions possible! This is handled in line 230.

If we were able to backtrack, then we set the column counter C to the column of the queen in this earlier row and we go back to line 110 so that we can look at solutions with a queen in the next column.
190 F=F+1
200 F$ = "Found " @ FOR K=1 TO S @ F$ = F$ & STR$(B(K)) @ NEXT K
210 F$ = F$ & " (" & STR$(F) & ")" @ DISP F$
We wind up here if we found a solution.

So, we start by incrementing the running count of solutions found. Lines 200 and 210 construct and display a string indicating the solution just found and the number of solutions found so far. It ends up looking something like:

Found 53142 (10)

The five digits after "Found" indicate that this was a 5x5 puzzle and that the tenth solution found (indicated by "(10)" at the end of the line) was to place a queen in each of columns 5, 3, 1, 4 and 2 of rows 1 to 5 respectively.
220 R=R-1 @ GOTO 110
Remember, we got here by incrementing R beyond the total number of rows in the puzzle (see line 90)? Before we go back into the puzzle to look at the next column of the last row, we need to decrement the row again to bring it back within the bounds of the problem. We then go back to line 110.
230 W=TIME-W @ IF W<0 THEN W=W+86400
This is the end of the puzzle!

We find out how long it took to crunch through all of the possible solutions by subtracting the time of day now from the time of day when we started in line 80. There is, however, one pitfall to this. TIME is a function that reads the clock to give you the time of day in the form of the number of seconds elapsed since midnight, to the nearest 100th of a second. This value can therefore be anything from 0.00 (dead on midnight) to 86399.99 (0.01 second before midnight). Reminder: disregarding leap seconds, there are 60x60x24=86400 seconds in a day.

What happens if you're like me and you sometimes leave a test to run its course overnight after you've gone to bed, resulting in the test starting before midnight and ending after midnight?

Say the test starts at 11:50pm and finishes at 12:10am. 11:50pm is 10 minutes (600 seconds) to midnight, so TIME will return 86400-600=85800. Come 12:10am and the end of the test, TIME will return 600, so TIME-W will be 600-85800=-85200. Congratulations on the negative test duration. You were able to complete the test 23h and 40 minutes before you started...

The HP-71B is a great machine, granted, but it hasn't quite mastered time travel yet. If the duration computed is negative then we compensate for the date flipping over by adding the number of seconds in a day to the result we found. -85200+86400=1200, which looks more like the 20 minutes expected.
240 DISP "SOLUT:" & STR$(F) & ", NODES:" & STR$(N)
We construct and display the string indicating how many solutions were found and how many nodes were examined to find them.

The short video below shows all of this in action, and the complete listing is provided below that.

As shown at the end of the video, the variable W contains the time in seconds taken to complete the puzzle, F is the number of solutions found, N is the number of nodes examined and F$ indicates the last solution found.
50 INPUT "Board size 1-9: ";S
60 IF S>9 OR S<1 THEN 50
70 INTEGER B(S) @ DISP "Working..."
80 R=0 @ F=0 @ N=0 @ W=TIME
90 R=R+1 @ IF R>S THEN 190
100 C=0
110 C=C+1 @ IF C>S THEN 170
120 T=0
130 N=N+1
140 T=T+1 @ IF T>=R THEN 160
150 IF B(T)=C OR ABS(B(T)-C)=R-T THEN 110 ELSE 140
160 B(R)=C @ GOTO 90
170 R=R-1 @ IF R=0 THEN 230
180 C=B(R) @ GOTO 110
190 F=F+1
200 F$ = "Found " @ FOR K=1 TO S @ F$ = F$ & STR$(B(K)) @ NEXT K
210 F$ = F$ & " (" & STR$(F) & ")" @ DISP F$
220 R=R-1 @ GOTO 110
230 W=TIME-W @ IF W<0 THEN W=W+86400
240 DISP "SOLUT:" & STR$(F) & ", NODES:" & STR$(N)

29 December 2018

Numerical solver - Newton-Raphson method

The program downloadable from this page is a numerical solver. It attempts to find the zeroes of a function $f(x)$, i.e. it attempts to find $x$ such that $f(x)=0$

There are several methods for this. One that requires little computing power is Newton's Method, or the Newton-Raphson method, a description of which can be found on Wikipedia.

We have no way of knowing the exact value of $f'(x)$ on the HP-41CX used here because that machine has no symbolic math capabilities, so we calculate an approximation of $f'(x)$ by taking a small value $\epsilon$ and calculating:
\[ f'(x) \approx \frac {f(x+\epsilon)-f(x)} \epsilon \]
This then allows us to calculate the next value of $x$ to use and display so that we see the calculator homing in on the zero that we're looking for:
\[ x_{n+1} = x_n - \frac {f(x_n)} {f'(x_n)} \]
Once we find a value of $x$ such that $|f(x)| \lt \epsilon$, we stop and that value of $x$ is deemed to be a zero of $f(x)$.

The same $\epsilon $ is used here as the small value added to $x$ to calculate $f'(x)$ and as the threshold below which we deem $f(x)$ to be zero. It is calculated according to the display precision of the HP-41CX. The calculator's flags 36-39 are examined to ascertain how many significant figures are displayed and we retrieve the digit after the last used FIX, SCI or ENG command. E.g., if the calculator is in FIX 4 mode then we get a 4 back. Let's call this number $p$ for precision.

We then calculate $\epsilon = 10^{-p} $, so the higher the display precision, the greater precision we seek from the algorithm. Similarly, if less precision is required for the display then less precision is demanded from the algorithm and it can complete after fewer iterations.

It may well be that the algorithm is unable to converge on a zero because of the nature of the function studied. It could be constant or it could send the algorithm off on a wild goose chase. There's some discussion of such cases on the Wikipedia page linked to above.

Finally, the link to download this HP-41 program is here: newton-raphson.zip

It is provided in text listing format, a .raw file and a printable PDF with wand bar codes.

Here is a video of a HP-41CX finding various zeros of the function:
\[ f(x) = x^3-2x^2-11x+12 = (x+3)(x-1)(x-4) \]
The three zeroes are therefore $x=-3, x=1$ and $x=4$.

And here is the equivalent in BASIC for HP-71:
30 P=FLAG(-17)+2*FLAG(-18)+4*FLAG(-19)+8*FLAG(-20) @ E=10^(-P) @ M=2*P+8
50 Y=FNF(X) @ IF ABS(Y)<E THEN 120
60 "X=" & STR$(X)
70 Y1=(FNF(E+X)-FNF(X))/E
90 X=X-Y/Y1
100 IF M>0 THEN M=M-1 @ GOTO 50
120 "ZERO: " & STR$(X) @ END
130 DEF FNF(X) = (X+3)*(X-1)*(X-4)
After running, the variable X contains the zero found and Y is the value of $f(x)$ calculated for that value of $x$ - this gives you the error, the absolute value of which should be lower than our $\epsilon $, which is stored in the variable E.