A programmer's HP-35s review

Some thoughts on the HP-35s after been programming it for a while: Doing the HP-35s Euclid Pack implementation.

Two features I learned to love when programming the 35s: Vector type and equations.

Vector type

The vector type has been critizised for missing key features:

1. No cross product. Must make a program.
   
2. No native way to unpack or address vector components. Must use 'dot product mask vectors'.
3. No operation to tell a vector's dimension (dim or size).

Despite these drawbacks (of which the two first can be worked around), I find the vector type a very valuable feature because of it's data packing capabilities.

Some focus has been on using 3D vectors to augment the machine's 800 indirect registers. For the things I do and the fact that it is a closed system (no I/O), I can't imagine applications using a large fraction of the 800 registers even without packing 3 real numbers in one register!

But, I find the data packing offered by 3D vectors to be very valuable in the 'more valuable real estate territories' of near end user variable registers (A-Z) and working registers (stack and lastx).

For example, for the design of the euclid pack project I ended up with wanting to store definition of two lines (four 2D or 3D points), two planes (two 3D points and two normal vectors) and one triangle (three 3D points) in variables.

If all data entered is in 3D, that would require minimum use 33 registers! And this is data I would like to be close to the end user, so I would like to use variables: Unless had vector type would not have room. Since have the vector type minimum register use is 11.

And implementation actual ended up storing additional data in variables to make certain computations easy and results handy. And I still did not use up all variable registers...

Also now we can have 12 real numbers on the stack! That becomes usefull when using programs that take more than 4 numbers as input. For example to define a 3D triangle (the T program) end user must give in nine numbers (three 3D numbers), without the ability to just put the points on the stack and run the program, the program would had to prompt for the data.

Also the vector type, RPN paradigm and the use of equation is very powerfull when programming:

Equations in programs

I find programming complex calculations in need of much intermediate results to be easy on the 35s by combining the RPN paradigm, vector type and equations:

• The vector type allows 12 values to be stored on the stack.
• Equations in programs can use all these number while being evaluated, only when done is maximum three numbers lost on the stack: A 3D vector in the T stack register.
• The result of evaluating an equation may very well be a 3D vector: A single line in a program may execute three equation which result raises the stack only once.

I think second half of the last program (lines 28 to 35) I made in this series realy utilize this programming style, there is eight lines of equations in sequence!

A very nice feature of the HP 35s that helps programming this way is that when in programming mode and to insert a stack register as a variable in an equation using the RDN key one sees the content of the registers user side. If one has a worked example of the problem being programmed and is programming equation by equation and running the program before programming a new equation then when programming the equation it is easy to see where on the stack needed numbers are.

One drawback of this style of programming is that the evaluation of equation makes for much slower execution then a pure RPN version. For programs that intention is to be called a large (or small...) number of times as subroutines by other programs this style of programming may not be appropiate.

On an improved HP 35s (HP 35s+ :-) perhaps one could have some kind of equation compilation?

No need for me to get into any of the now well known short commings or bugs of the device documented by the community of the MoHPC here, except to say that the reported vector syntax bug is indeed real on some machines and makes this type of programming much harder than needed. I know, I got one... Still I have been able to do whatever I did without ever needing to reset the thing.

HP-35s Euclid Implementation - XV : The inscribed and circumscribed circle of a triangle

The final program to fullfill the requirements! Problem 16 and 17: Finding the inscribed (incircle) and circumscribed (circumcircle) circles of a triangle.

Finding the inscribed circle is the same as solving the problem of finding the circle tangent to three lines (lines defined by the triangle edges) and finding the circumscribed circle is the same as solving the problem of finding the circle that passes through three points (triangle's vertexes).

The triangle in question has been entered using the T program.

                  v3
                 /  \
                /    \
               /      \             
              v1------v2

Refering to the above triangle I give the formulas used in the HP-35s program presented here:

Inscribed circle:

The barycentric coordinates for the incenter are:

(l1/p, l2/p, l3/p)

Where p is the perimeter of the triangle (sum of all edges lengths) and l1 is the length of the edge v3-v2, l2 is the length of the edge v1-v3 and l3 is the length of the edge v2-v1.

The radius of the inscribed circle is:

(2A)/p

Where A is the area of the triangle.

Circumscribed circle:

With the intermediate results

d1=-e2*e3
d2=-e3*e1
d3=-e1*e2

(e1=v3-v2, e2=v1-v3, e3=v2-v1, * is the dot product)

c1=d2d3
c2=d3d1
c3=d1d2

c=c1+c2+c3

the barycentric coordinates for the circumcenter are:

((c2+c3)/2c, (c3+c1)/2c, (c1+c2)/2c)

The radius of the circumscribed circle is:

sqrt((d1+d2)(d2+d3)(d3+d1)/c)/2

These formulas are taken from the excelent book 3D Math Primer for Graphics and Game Development that also provide background for much of the geometric problems in the Euclid Pack specification.

Stack Input/Output:

No input on stack-> XEQ O ->[[bc1, bc2, bc3], cr, [bi1, bi2, bi3], ir | L]

Where (bc1, bc2, bc3) are the barycentric coordinates for the circumscribed circle's center (circumcenter) , cr is its radius, (bi1, bi2, bi3) are the barycentric coordinates for the inscribed circle's center (incenter) and ir is its radius.

Variables:

Reads:

K : First vertex (3D vector). Populated by the T program.
L : Second vertex (3D vector). Populated by the T program.
M : Third vertex (3D vector). Populated by the T program.

Writes:

Z : Scratch data.

Program:

O001 LBL O
O002 5
O003 STO I
O004 STO(I)
O005 4
O006 STO I
O007 LASTX
O008 STO(I)
O009 eq M-L
0O10 eq M-K
O011 eq [ABS(REGY),ABS(REGX),ABS(L-K)]
O012 STO Z
O013 eq REGX/([1,0,0]xREGX+[0,1,0]xREGX+[0,0,1]xREGX)
O014 R^
O015 R^
O016 XEQ X001
O017 R^
O018 X<>Y
O019 ABS
O020 eq REGX/([1,0,0]xZ+[0,1,0]xZ+[0,0,1]xZ)
O021 X<>Y
O022 RDN
O023 X<>Y
O024 RCL(I)
O025 ABS
O026 RDN
O027 XEQ U001
O028 eq M-L
O029 eq K-M
O030 eq L-K
O031 eq [-REGYxREGX,-REGXxREGZ,-REGZxREGY]
O032 eq [([0,1,0]xREGX)x([0,1,0]xREGX),
         ([0,0,1]xREGX)x([1,0,0]xREGX),
         ([1,0,0]xREGX)x([0,1,0]xREGX)]
O033 eq [1,0,0]xREGX+[0,1,0]xREGX+[0,0,1]xREGX
O034 eq 0.5xSQRT((([1,0,0]xREGZ+[0,1,0]xREGZ)x
                  ([0,1,0]xREGZ+[0,0,1]xREGZ)x
                  ([0,0,1]xREGZ+[1,0,0]xREGZ))/REGX)
O035 eq [([0,1,0]xREGZ)+([0,0,1]xREGZ),
         ([0,0,1]xREGZ)+([1,0,0]xREGZ),
         ([1,0,0]xREGZ)+([0,1,0]xREGZ)]/(2xREGY)
O036 XEQ U079
O037 RTN

Comments:

Terms of use.

Mnemonics: O since letter is a circle (or more or less depeding on font, rather circular on the 35s).

Reason for not converting to cartesian coordinates is that since the barycentric coordinates are relative to the vertices it is much easier to 'see' where the center points are relative to the triangle. If one where to be presented the cartesian coordinates one must remember the nine vertex coordinates to know where (relative to triangle) the center points are. Particular one see at once if the circumscribed circle's center is inside or outside the triangle.

Cartesian coordinates are easily obtained using the A program by the following key sequence:

XEQ A
RDN
RDN
XEQ A
RDN
RDN

The last two RDN are only provided to get content of stack in same order as before except that barycentric positional vectors are now replaced with cartesian ones.

HP-35s Euclid Implementation - XIV : Barycentric to Cartesian coordinates

Going from barycentric to cartesian coordinates (problem 15 in the requirements) is the easy part of Cartesian<->Barycentric 2D/3D conversion, this will be a much shorter blog entry than this one!

Stack Input/Output:

3D case: [[b1, b2, b3], Y, Z, T | L]-> XEQ A ->[[x, y, z], Y, Z, T | [b1, b2, b3]]

2D case: [[b1, b2, b3], Y, Z, T | L]-> XEQ A ->[[x, y], Y, Z, T | [b1, b2, b3]]

Variables:

Reads:

K : First vertex (3D vector), in 3D case. Populated by the T program.
L : Second vertex (3D vector), in 3D case. Populated by the T program.
M : Third vertex (3D vector), in 3D case. Populated by the T program.
U : First vertex of triangle projected to 2D (2D vector), in 2D case. Populated by the T program.
V : Second vertex of triangle projected to 2D (2D vector), in 2D case. Populated by the T program.
W : Third vertex of triangle projected to 2D (2D vector), in 2D case. Populated by the T program.

Program:

A001 LBL A
A002 ABS
A003 RDN
A004 FS? 2
A005 GTO A008
A006 eq Kx([1,0,0]xLASTX)+Lx([0,1,0]xLASTX)+Mx([0,0,1]xLASTX)
A007 RTN
A008 eq Ux([1,0,0]xLASTX)+Vx([0,1,0]xLASTX)+Wx([0,0,1]xLASTX)
A009 RTN

Comments:

Terms of use.

Mnemonics: None realy, but letter is left to B, the program label for Barycentric to Cartesian conversion.

If we are working in 2D (flag 2 set) we are applying the weights (barycentric coordinates) to the 2D coordinate's of the projected triangle (line A008) else applying the weights to the vertexes 3D coordinates (line A006).

Fixing the HP-41CV - II : Sending it

Finaly got around to post the HP-41CV to Fix That Calc!

HP-35s Euclid Implementation - XIII : Cartesian to Barycentric coordinates

Arrived at problem 14 of the requirements: Computing barycentric coordinates given cartesian coordinates.

Barycentric coordinates [b1, b2, b3] of a point in a plane are coordinates defined by a triangle in the plane (defines the plane) where the triangle's first vertex's barycentric coordinates are [1,0,0], second vertex's are [0,1,0], third vertex's are [0,0,1] and one have that b1+b2+b3=1.

Coordinates of any point in the plane can be expressed as the weighted average of the triangle's vertices coordinates, these weights are the barycentric coordinates:

 Barycentric     Cartesian
(b1, b2, b3)<=>b1v1+b2v2+b3v3

Where vi is the positional vector for vertex i.

Barycentric coordinates has several uses and are of theoretical importance outside the scope of this blog entry. Here I give two practical uses:

1. Provides a robust test for if a point in triangle's plane with barycentric coordinates [b1, b2, b3] is inside the triangle, on one of the triangle's edge or outside the triangle: If all components are in the range[0,1] point is inside, if bi is 0 and the two others are in the [0,1] range then the point is on the edge oposite of triangle's vertice i, and finaly if at least one component is not in the [0,1] range then point is outside triangle.
2. If a function's values is known at the triangle's vertices one can perform linear interpolation inside the triangle at the point: v1*b1+v2*b2+v3*b3. Where vi is the value at vertex i and [b1, b2, b3] is the barycentric coordinates for point to evaluate at. Note that converting from barycentric coordinates to cartesian coordinates is a special case of this.
   

I choose to convert from cartesian coordinates to barycentric by programming the solution to the equations for the 2D cartesian to barycentric coordinates case:

Q=(y1-y3)*(x2-x3)+(y2-y3)*(x2-x1)

b1=((y-y3)*(x2-x3)+(y2-y3)*(x3-x))/Q
b2=((y-y1)*(x3-x1)+(y3-y1)*(x1-x))/Q
b3=((y-y2)*(x1-x2)+(y1-y2)*(x2-x))/Q

Where (xi, yi) are cartesian coordinates for triangle's i'th vertex and (x, y) is the cartesian coordinates of the point in question. This is the solution to the equation set of three uknown (b1, b2, b3):

x=b1*x1+b2*x2+b3*x3
y=b1*y1+b2*y2+b3*y3
1=b1+b2+b3

To solve 3D cartesian to barycentric I project the triangle and point in question to 2D and solve as a 2D problem.

Implementation is done in two steps:

1. Modify the T program so it: A) Also accept a 2D triangle (flag 2 set='2D mode'. B) Project the triangle to the cardinal plane that gives the largest area, this is the triangle used for computing barycentric coordinates. C) Computes the quotient Q in the above equations.
2. Implement a program labeled B to accepts a 2D or 3D cartesian tuple and use result from B and C above to compute barycentric coordinates.

Step 1 : The modified T program

Stack Input/Output:

2D case: [[x2, y2], [x1, y1], [x0, y0] | L]->XEQ T->[[nx, ny, nz], [x2, y2], [x1, y1], [x0, y0] | L]

3D case: [[x2, y2, z2], [x1, y1, z1], [x0, y0, z0] | L]->XEQ T->[[nx, ny, nz], [x2, y2, z2], [x1, y1, z1], [x0, y0, z0] | L]

Where v0=(x0, y0, z0) is first, v1=(x1, y1, z1) is second, v2=(x2, y2, z2) is the third triangle vertex and [nx, ny, nz] is the cross product (v1-v0)x(v2-v0).

Variables:

Writes:

K : First vertex (3D vector).
L : Second vertex (3D vector).
M : Third vertex (3D vector).
Q : Quotient (y1-y3)*(x2-x3)+(y2-y3)*(x2-x1).
U : First vertex of triangle projected to 2D (2D vector).
V : Second vertex of triangle projected to 2D (2D vector).
W : Third vertex of triangle projected to 2D (2D vector).
X : Index of coordinate discarded in 3D->2D projection.

Program:

T001 LBL T
T002 XEQ U001
T003 FS? 2
T004 XEQ X068
T005 STO M
T006 RDN
T007 FS? 2
T008 XEQ X068
T009 STO L
T010 RDN
T011 FS? 2
T012 XEQ X068
T013 STO K
T014 R^
T015 X<>Y
T016 -
T017 RDN
T018 RDN
T019 LASTX
T020 -
T021 X<>Y
T022 RDN
T023 XEQ X001
T024 ABS
T025 1E-6
T026 X>Y?
T027 GTO T065
T028 LASTX
T029 XEQ X025
T030 STO X
T031 RCL K
T032 X<>Y
T033 XEQ X052
T034 STO U
T035 RDN
T036 RCL L
T037 X<>Y
T038 XEQ X052
T039 STO V
T040 RDN
T041 RCL M
T042 X<>Y
T043 XEQ X052
T044 STO W
T045 eq ([0,1]xU-[0,1]xW)x([1,0]xV-[1,0]xW)+
        ([0,1]xV-[0,1]xW)x([1,0]xW-[1,0]xU)
T046 STO Q
T047 LASTX
T048 4
T049 STO I
T050 RDN
T051 RCL(I)
T052 ABS
TO53 RDN
T054 FS? 2
T055 GTO T060
T056 RCL K
T057 RCL L
T058 RCL M
T059 GTO T063
T060 RCL U
T061 RCL V
T062 RCL W
T063 R^
T064 RTN
T065 SF 10
T066 eq DEGENERATED
T067 PSE
T068 CF 10
T069 XEQ U029
T070 RTN

Comments:

Terms of use.

Mnemonics: T for Triangle.

Uses the X068 subroutine to place 2D triangle in th XZ plane if flag 2 set and the X025 routine to find out which coordinate to discard for the 3D to 2D projection and finaly the X052 routine to actual do the 3D to 2D conversion. In the 2D case this is a rather roundabout way to get end user's input into the variables U, V and W! But I think this worked out quite well in order to have one program for both the 2D and 3D case:

One can compute barycentric coordinates in 3D using four cross product, four dot product and three scalar division operations avoiding projection from 3D to 2D, but involves more scalar math operations than the implemented method. On a vector architecture it may be a better choice, on the HP-35s architecture I think this is a good solution!

The index of the coordinate discarded to do the 3D to 2D projection is remembered in variabe X so the next program to actual compute barycentric coordinates know which coordinate to discard from input.

An unrelated change done to the T program is that it now check if the computed plane normal vector is considered to be of zero length: Then all points are the same or colinear, the message "DEGENERATED" is then shortly displayed and stack registers restored.

Note that one get a cross product calculated normal vector also in the 2D case (direction
[0, 1, 0]), length still twice triangle's area.

Step 2 : Program to convert from Cartesian to Barycentric Coordinates

Stack Input/Output:

[[x, y, z], Y, Z, T | L]->XEQ B->[[b1, b2, b3], Y, Z, T | [x, y, z]]

Where (x, y, z) are cartesian coordinates for point in question and (b1, b2, b3) are computed barycentric coordinates.

Variables:

Reads:

Q : Quotient (y1-y3)*(x2-x3)+(y2-y3)*(x2-x1). Populated by the P program.
U : First vertex of triangle projected to 2D (2D vector). Populated by the P program.
V : Second vertex of triangle projected to 2D (2D vector). Populated by the P program.
W : Third vertex of triangle projected to 2D (2D vector). Populated by the P program.
X : Index of coordinate discarded in 3D->2D projection. Populated by the P program. 

Program:

B001 LBL B
B002 XEQ U001
B003 FS? 2
B004 GTO B007
B005 RCL X
B006 XEQ X052
B007 eq [([0,1]xREGX-[0,1]xW)x([1,0]xV-[1,0]xW)+
         ([0,1]xV-[0,1]xW)x([1,0]xW-[1,0]xREGX),

         ([0,1]xREGX-[0,1]xU)x([1,0]xW-[1,0]xU)+
         ([0,1]xW-[0,1]xU)x([1,0]xU-[1,0]xREGX),

         ([0,1]xREGX-[0,1]xV)x([1,0]xU-[1,0]xV)+
         ([0,1]xU-[0,1]xV)x([1,0]xV-[1,0]xREGX)]/Q
B008 XEQ U050
B009 RTN

Comments:

Terms of use.

Mnemonics: B for Barycentric.

If the input is a 2D (flag 2 set) then goes straight to the equation that calculates baycentric coordinates using the data precomputed by the modified T program, else first project the 3D input to 2D using subroutine X052 and information found in the X variable.

HP-35s 2D to 3D vector conversion

Short program to convert a 2D vector to a 3D by placing it in the z=0 plane. Usefull when want to utilize 3D solutions to 2D data, for example taking cross product of two 2D vectors.

Stack Input/Output:

[[u, v], Y, Z, T | L]-> XEQ 068 ->[[u, v, 0], Y, Z, T | [u, v]]

Program:

X068 ABS
X069 RDN
X070 eq [[1,0]xLASTX,[0,1]xLASTX,0]
X071 RTN

Comments:

Terms of use.

Placed in the X 'library' since a vector utlity, from before have projecting to a cardinal plane, finding the component of a 3D vector with max absolute value, vector unpack routines, cross product and normalization.

HP-35s 3D to 2D projection utility subroutines - II

Subroutine to project a 3D vector to 2D by paralell projection to one of the cardinal planes; discards one coordinate.

Stack Input/Output

[n, [u,v,w], Z, T | L]-> X047 ->[[x,y], n, [u,v,w], Z | L]

Where n tells which coordinate to discard: [x=v, y=w] if n=1, [x=u, y=w] if n=2 and [x=u, y=v] if n=3.

Program:

X052 1
X053 X=Y?
X054 GTO X062
X055 RDN
X056 2
X057 X=Y?
X058 GTO X065
X059 RDN
X060 eq [[1,0,0]xREGY,[0,1,0]xREGY]
X061 RTN
X062 RDN
X063 eq [[0,1,0]xREGY,[0,0,1]xREGY]
X064 RTN
X065 RDN
X066 eq [[1,0,0]xREGY,[0,0,0]xREGY]
X067 RTN

Comments:

Terms of use.

The parameter n is typical found using the X025 subroutine.

Placed in the X 'library' since a vector utlity, from before have finding the component of a 3D vector with max absolute value, vector unpack routines, cross product and normalization.

Change history:

20071014:1723UTC : Because of change in the X025 subroutine line numbers changed here. Entry point used to be X047 but is now X052.