The Calculator

The HP-16C is a versatile and unique calculator, especially designed for the many professionals and students who work with computers and microprocessors—whether as programmers or designers. The HP-16C specialized design provides:

Integer arithmetic in four number bases (hexadecimal, decimal, octal and binary), operating in 1's or 2's Complement or Unsigned mode.
A variable word size, selected by the user, up to a maximum of 64 bits.
Logical operators and bit manipulations.
203 bytes of user memory, providing up to 203 program lines.
Floating-point decimal arithmetic.
Continuous Memory, retaining data and program instructions indefinitely.
Extremely low power consumption and long battery life.

– from the introduction to the HP-16C Computer Scientist Owner's Handbook.

A sketch of the HP-16C

The 16C came in a package that could easily slip into (and out of) a shirt pocket and weighed just under 115g. Not that impressive by today's standards perhaps but this was 1982—a time when a 20M Byte "winchester" hard drive was the size (and weight) of three house bricks. It shared its technology and packaging with the four other calculators that made up the 10C, or Voyager, series.

Architecture

The 16C was both a standard HP BCD-based calculator and an abstraction of a CPU offering a basic instruction set and a small set of flags. Memory was dynamically partitioned into addressable storage and program instructions. The partitioning scheme and other important architectural elements were summarised on this diagram from the Owner's Handbook.

Programs

The net is replete with programs for the HP-48 and the HP-41. Unfortunately there are very few for the HP-16C. This section includes listings of some of my own and any that have been contributed by others. I am particularly interested in useful programs that you may have written for your calculator. If you are willing to share any of your gems, send me a listing via email and I'll make them available here.

NOTE: unless otherwise attributed, all program listings are mine. Feel free to use them but please don't republish them without the permission of the author.

Factorial

Useful for statistics and Taylor series approximations of the trig functions. This implementation uses only 17 steps and preserves both Y and I. It works in both floating-point and integer modes.


; Uses I but preserves its value on entry and retores it
; on exit.

; Uses flag 5--because it will be cleared by the
; multiplications.

; Preserves and restores Y.

+01	LBL FACT

+02	CF 5
+03	X=0?
+04	SF 5
+05	X-I	;	Ie	Ye	Ze	T	N
+06	[1]	;	1	Ie	Ye	Ze	N
+07	F? 5
+08	ISZ

+09	LBL LOOP
+10	RCL I	;	N-1	N	Ie	Ye	N
+11	[X]	;	N=N-1*N	Ie	Ye	Ye	N
+12	DSZ	;	N	Ie	Ye	Ye	N-1
+13	GTO LOOP

+14	X-Y	;	Ie	N	Ye	Ye	0
+15	X-I	;	0	N	Ye	Ye	Ie
+16	Rv	;	N	Ye	Ye	0	Ie
+17	RTN

Basic Trigonometry Library

Here I present a library of the basic trigonometric functions: sin x, cos x, tan x and their inverses.

The obvious question is, "Why bother?" Put it down to a desire to document how things were done in the old days. I find the code to be quite instructive. It gives an insight into how far we've come in the last 20 years. I am particularly proud of the polynomial "engine" in this library. A fine example of 1980's-style code re-use. Not to mention the fact that I need routines like these to exercise my simulator.


Trig in 168 steps. Unfortunately it uses all 16 labels.
It also uses four registers and all three "user" flags.
Sorry.

On the real 16C, that's 168 + 4*7 = 196 bytes. This leaves
you with 1 register to play with while you're using it. If
you don't mind re-keying Pi into R0 you have 5 registers
available while it's parked in your calculator and not being
used.

If you load it into the simulator from my cp16 files, R0
will be primed with a suitable approximation of Pi. If
you're keying it into a real 16C, you'll need to put your
preferred value of Pi into R0 before you use it in "degree"
mode.

Perhaps needless to say, the algorithm I'm using--evaluation
of Taylor polynomials--take a *very* long time on the real
calculator. You can reduce the time, at the expense of
accuracy, by changing the (related) constants at lines 24
and 57. The constant at line 24 controls the number of terms
in the polynomial to evaluate. Line 57 must be one more than
line 24.

The behaviour of these routines is undefined in integer mode.

The accuracy leaves something to be desired. It returns
useful results for the following values of X:

sin	0 <= X <= Pi/2
cos	0 <= X <  Pi/2 (Pi/2 - 5e-8)
tan	0 <= X <  Pi/2

asin	0 <= X < sqrt(2)/2
acos	0 <  X < sqrt(2)/2
atan	0 <= X < sqrt(2)/2

Flag 3 controls interpretation of angles:

	  set = degrees
	clear = radians

Entry Points

A : sin(X)
B : cos(X)
C : tan(X)

D : asin(X)
E : acos(X)
F : atan(X)

Utility Routines

9 : DEG -> RAD
8 : RAD -> DEG
7 : X!
6 : (iterator)
5 : Y^X
4 : (iterator)

Helper Routines

3 : POLY
2 : (iterator)
1 : TERM
0 : DONE

Registers

R0 : Pi - you must prime this before you start

R1 : POLY constant (X)
R2 : POLY accumulator (N)
R3 : TAN scratch


Copyright (C) 1983-1984, 2004 Cameron Paine


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; SIN

LBL SIN [A]	; [.001]

X=0?
; special case: sin(0) = 0
RTN

; Configure POLY to do sin.
SF 0
SF 1

; If we're in degrees mode, convert to radians.
F? 3
GSB RADIANS
GTO POLY


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; COS

LBL COS [B]	; [.009]

CF 1
X=0?
SF 1
[1]
F? 1
; special case: cos(0) = 1
GTO DONE
Rv

; Configure POLY to do cos
SF 1

; If we're in degrees mode, convert to radians.
F? 3
GSB RADIANS
; Fall through to POLY...


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; POLY

; Evaluate the Taylor/Maclaurin polynomials.

; sin(x)  = x - x^3/3! + x^5/5! - x^7/7! + ...
; cos(x)  = 1 - x^2/2! + x^4/4! - x^6/6! + ...
; atan(x) = x - x^3/3  + x^5/5  - x^7/7  + ...

; e^x     = 1 + x^2/2! + x^3/3! + x^4/4! + x^5/5! + ...

; sin and cos are the Taylor series.
; atan is the Maclaurin series.

; Configuration truth table

;	F0	F1
;	1	1	sin
;	0	1	cos
;	1	0	atan

; F2 is a flip-flop that controls addition/subtraction

; R1 : POLY constant (X)
; R2 : POLY accumulator (N)

LBL POLY [3]	; [.020]

; Initialise series (for sin/atan)
STO 2	; N
SF 2	; start with subtraction

; Need the input value as a constant
STO 1	; X

; Prime loop counter with the number of terms to evaluate.
; If you change this you need to change the related
; constant in TERM.
[8]	; [.024]
CHS
STO I

F? 0	; sin?
GTO LOOP

; The cos series starts with 1
[1]
STO 2

LBL LOOP [2]	; [.031]

RCL 1	;	X
GSB TERM;	t	X
GSB POW	;	X^t
GSB TERM;	t	X^t

;;;;; R/S ; [.033] ;;;;;;;;;;;;;;;;;;;;;;;

F? 1
GSB FACT;	t!	X^t
[/]	;	X^t/t!

F? 2
CHS
X<0?
CF 2
X>0?
SF 2
RCL 2	;	N	X^t/t!
[+]	;	N+X^t/t!
STO 2

ISZ
GTO LOOP

LBL DONE [0]	; [.050]

CF 0
CF 1
CF 2

RTN 		; [.054]

LBL TERM [1]	; [.055]

; we need 2, 4, 6... from -7, -6, -5...
;
; 8+I gives: 1, 2, 3...
;
; The constant therefore must be 1 more than the constant
; at [.024].

RCL I	;	I	Ye	Ze	Te
[9]	;	9	I	Ye	Ze	[.057]
[+]	;	t=9+I	Ye	Ze	Ze
[2]	;	2	t	Ye	Ze
[*]	;	t=t*2	Ye	Ze	Ze
F? 0
[1]	;	1	t	Ye	Ze
F? 0
[+]	;	t=t+1	Ye	Ze	Ze

RTN	; [.065]

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; TAN

; Does tan(X), X in degrees or radians

; We do sin(x)/cos(x). We rely on SINCOS to convert
; degrees to radians if necessary.
;
LBL TAN [C]	; [.066]

STO 3
GSB COS
RCL 3
X-Y
STO 3
Rv
GSB SIN
RCL 3
[/]

RTN		; [.076]

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; ARCSINCOS

; arcsin(x) = arctan(x / sqrt(1 - x^2))
; arccos(x) = arctan(sqrt(1 - x^2) / x)

; We use flag 0, which will be re-used/reset by atan.

; This is currently pretty limited. It's reasonably
; correct for 0 < X < sqrt(2)/2
;
; BUGS: asin(1) and acos(0) will return zero. I don't
; have enough program storage to handle these cases
; more gracefully. I return zero so we unwind the stack
; rather than causing Error 0.

LBL ACOS [E]	; [.077]
SF 0

LBL ASIN [D]	; [.079]

ENTER	;	X	X
ENTER	;	X	X	X
[1]	;	1	X	X
X-Y	;	X	1	X
ENTER	;	X	X	1	X
[X]	;	X^2	1	X	X
[-]	;	1-X^2	X	X	X
SQRT
F? 0
X-Y
X=0?
; Special case: asin(1) and acos(0) gives zero divisor
RTN
[/]
GSB ATAN

RTN		; [.094]


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; ATAN

; We select one of two different approximations depending
; on the input value.

; Error computing atan(tan(X))
;
; X		HP-41CX		Approx1		Approx2
; Pi/32		0		0		N/A
; Pi/16		0		0		N/A
; Pi/8		0		0.000000139	0
; Pi/4 atan(1)	0		1.586710373	2.786537-07
; Pi/3		0.000000010	2.783317+04	1.359862-04
; Pi/2.5	0		6.004617+08	9.445856-03
; Pi/2.25	0		2.114579+13	1.316556-01

LBL ATAN [F]	; [.095]
X=0?
; special case: atan(0) == 0
RTN

; Condition POLY for atan
SF 0

; Algorithm selection threshold
[2]
SQRT
[4]
[/]
X-Y

; Select algorithm
X<=Y?
GTO F	; LABEL REUSE >>>>>>>>>>>>>>>>>>>>>>>>>>>

; We use this approximation if X is > the threshold above.
; arctan(x) = 2*arctan((sqrt(1+x^2)-1)/x). We don't recurse
; because: 1) we reused LBL F; 2) we'd potentially have to
; negate the RAD/DEG conversion.

ENTER	;	X	X
ENTER	;	X	X	X
[*]	;	X^2	X
[1]	;	1	X^2	X
[+]	;	s=1+X^2	X
SQRT	;	s	X
[1]	;	1	s	X
[-]	;	s-1	X
X-Y	;	X	s-1
[/]	;	s-1/x
GSB POLY
[2]
[*]

F? 3
GSB DEGREES
RTN	; [.121]

LBL F	; LABEL REUSE <<<<<<<<<<<<<<<<<<<<<<<<<<<

GSB POLY
F? 3
GSB DEGREES
RTN	; [.126]


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; RADIANS

; Converts degrees to radians

; Assumes pi is in R0

; Preserves Y and Z.

; Result in X

LBL RADIANS [9]	; [.127]

RCL 0	;	pi	a	Ye	Ze
[*]	;	a*pi	Ye	Ze	Ze
[1]
[8]
[0]	;	180	a*pi	Ye	Ze
[/]	;	rad	Ye	Ze	Ze

RTN		; [.134]


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; DEGREES

; Converts radians to degrees

; Assumes pi is in R0

; Preserves Y and Z.

; Result in X

LBL DEGREES [8]	; [.135]

[1]
[8]
[0]
[*]
RCL 0
[/]

RTN		; [.142]


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; FACTORIAL

; Does X! for X > 0

; Uses I but preserves its value on entry and retores it
; on exit.

; Preserves and restores Y.

; BUGS:
; This is a modified version of my factorial-17. It does
; not handle zero, which saves us 5 steps. Since it is
; for use in the Taylor series we won't need zero.

+01	LBL FACT [7]	; [.143]

+02	X-I	;	Ie	Ye	Ze	T	N
+03	[1]	;	1	Ie	Ye	Ze	N

+04	LBL LOOP [6]	; [.146]
+05	RCL I	;	N-1	N	Ie	Ye	N
+06	[X]	;	N=N-1*N	Ie	Ye	Ye	N
+07	DSZ	;	N	Ie	Ye	Ye	N-1
+08	GTO LOOP

+09	X-Y	;	Ie	N	Ye	Ye	0
+10	X-I	;	0	N	Ye	Ye	Ie
+11	Rv	;	N	Ye	Ye	0	Ie

+12	RTN		; [.154]


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; POW

; Does Y**X for X > 0

; Preserves I.

; BUGS:
; This is a modified version of my generic power routine.
; We only do X > 0 because that's all we need for the
; Taylor series. This saves us 7 steps and a label.

+01	LBL POW [5]	; [.155]

+02	X-I	;	Ie	Y	Ze	Te
+03	X-Y	;	Y	Ie	Ze	Te

; Generate 1 and put Y into LSTx
+04	ENTER	;	Y	Y	Ie	Ze
+05	[/]	;	1	Ie	Ze	Ze

+06	LBL LOOP [4]	; [.160]
+07	LSTx	;	Y	N	Ie	Ze
+08	[X]	;	N	Ie	Ze	Ze
+09	DSZ
+10	GTO LOOP

+11	X-Y	;	Ie	N	Ze	Ze
+12	X-I	;	0	N	Ze	Ze
+13	Rv	;	N	Ze	Ze	0

+24	RTN		; [.168]

The trigonometry library had its genesis in days before microprocessors had any native real number capabilities. I wrote the code while exploring the speed/space/accuracy trade-offs of various implementations for doing 2D graphics on the Cromemco "Super Dazzler" graphics system.

Internet Protocol (Version 4) Addresses

An IPv4 address is simply a 32-bit unsigned integer. An address is given a sub-structure by conventional rules that divide the address pool into three classes based on numerical range. The sub-structure can be further refined by applying a 32-bit mask to the address.

Conventionally, addresses and masks are written using dotted-decimal notation. In this system, each of the four octets that make up the address (or the mask) are expressed in decimal, with a dot separating each octet. The left-hand position holds the most-significant octet. For example:

192.168.192.42 255.255.192.0

There is an alternative notation for network masks in which a single integer, representing the number of bits that are set in the mask, is appended to a dotted-decimal address, preceded by a slash character. This is sometimes called CIDR notation. Using this scheme, the above address/mask combination would be expressed as:

192.168.192.42/18

This is a useful notation because the actual mask can be conveniently generated using the number after the slash and the 16C's left mask (MASKL) operator.

The Package

The routines in this package manipulate IPv4 addresses and network masks in various useful ways. None of the functions provided are particularly ingenious but collectively they can save quite a few keystrokes.

User Routines

GSB A - JOIN dotted decimal octets to 32-bit integer GSB B - SPLIT 32-bit integer to dotted decimal octets GSB C - NETwork address from CIDR GSB D - BCAST (broadcast) address from CIDR GSB E - Numnber of host ADDResses given a mask

Service Routines

GSB 9 - convert MASK in register 0 to CIDR form GSB 8 - OCTET processor for SPLIT GSB 7 - convert mask in X to CIDR MASK in X

Input and Output

The 16C has no native dotted-decimal representation. Therefore we use the four stack registers as follows:

T	octet 3 (ms)
Z	octet 2
Y	octet 1
X	octet 0 (ls)

For those functions that require CIDR addresses, the mask should be stored in Register 0. If the mask is in binary form it will be converted to CIDR form as required. Thus 192.168.192.42/26 is entered as:

26 STO 0 192 ENTER 168 ENTER 192 ENTER 42

The four octets in a returned address can be easily retrieved by using the roll-down key.

All routines use I for temporary storage. They may be trivially modified to use another register for that purpose. The calculator must be in unsigned decimal mode with a 32-bit word size.


; JOIN - dotted decimal octets to 32-bit integer

+01 - LBL JOIN (A)
+02 - R^
+03 - STO I	; I = o3
+04 - Rv
+05 - X-Y	; x = o1; y = o0
+06 - [8]	; discards T
+07 - RLn	; x = o1 << 8
+08 - LSTx	; bring back 8...
+09 - Rv	; ...and put it in T
+10 - OR	; x = o0 | o1 << 8; y = o2; z = 8; t = 8
+11 - X-Y	; x = o2
+12 - R^	; bring back 8...
+13 - SL	; ...make it 16
+14 - RLn	; x = o2 << 16
+15 - LSTx	; bring back 16...
+16 - R^	; bring back 8...
+17 - OR	; ...make it 24
+18 - Rv	; x = o2 << 16; y = o0 | o1 << 8; t = 24
+19 - OR	; x = o2 << 16 | o0 | o1 << 8
+20 - RCL I	; x = o3
+21 - R^	; x = 24
+22 - RLn	; x = o3 << 24
+23 - OR	; x = o3 << 24 | o2 << 16 | o0 | o1 << 8
+24 - RTN

; SPLIT - 32-bit integer to dotted decimal octets

+25 - LBL SPLIT (B)
+26 - GSB OCTET	; do o0
+27 - Rv
+28 - STO I	; stash in I
+29 - R^
+30 - GSB OCTET	; do o1
+31 - Rv
+32 - X-I	; stash o1, stack o0
+33 - R^
+34 - GSB OCTET	; do o2 and o3
+35 - X-Y	; x = o2; y = o3
+36 - R^
+37 - X-I	; x = o1; y = o2; z = o3
+38 - R^	; x = o0; y = o1; z = o2; t = o3
+39 - RTN

; NET - network address from CIDR host address

+40 - LBL NET (C)
+41 - GSB JOIN
+42 - GSB MASK
+43 - MASKL
+44 - AND
+45 - GTO SPLIT

; BCAST - broadcast address from CIDR network address

+46 - LBL BCAST (D)
+47 - GSB JOIN
+48 - GSB MASK
+49 - MASKL
+50 - NOT
+51 - OR
+52 - GTO SPLIT

; NADDR - number of host addresses on a network

+53 - LBL NADDR (E)
+54 - X=0?
+55 - RCL 0
+56 - X=0?
+57 - RTN
+58 - GSB CIDR_MASK (7)
+59 - MASKL
+60 - NOT
+61 - [1]
+62 - [-]
+63 - RTN

;
; Service Routines
;

; If the high bit of X is set we assume it contains a binary
; mask and we convert it to CIDR form. Otherwise we leave it
; alone.

+64 - LBL CIDR_MASK (7)
+65 - [3]
+66 - [1]
+67 - B?
+68 - B#
+69 - RTN

; Places the low-order octet of X in Y and the remaining
; high-order bits, right justified, in X.

+70 - LBL OCTET (8)
+71 - [8]
+72 - MASKR	; low octet mask
+73 - X-Y	; swap so we stash in LSTx
+74 - AND	; clear high bits
+75 - LSTx	; X = word; Y = octet;
+76 - [8]	; Clear low octet so we can rotate
+77 - MASKR	; low octet mask...
+78 - NOT	; ...becomes remaining octet mask
+79 - AND	; clear low octet
+80 - [8]
+81 - RRn	; rotate cleared octet to MSB
+82 - RTN

+83 - LBL MASK (9)
+84 - RCL 0
+85 - GTO CIDR_MASK (7)

Changing ends

If you spend any time working with network interfaces on Little Endian systems (like those powered by the Intel x86) it is often convenient to be able to re-order the bytes in 16- and 32-bit integers. This handy little program will do that for you. As a bonus, it works correctly with WSIZE set to either 16 or 32 in HEX or DEC mode.


; Swap the endian disposition of a 16- or 32-bit value.
; Works in either HEX or DEC integer mode.
; Uses 3 labels--but they provide useful standalone
; routines.
; Uses no storage registers
; Preserves Y.
; Destroys I.
; Z and T are lost.
; Result returned in X.

; It can be trivially modified to preserve I at the
; expense of Y. This will save 6 steps.

+01	43,22,A	LBL SWEND {A}
+02	1	[1]
+03	36	ENTER
+04	8	[8]
+05	42 A	SL
+06	42 E	RLn
+07	33	Rv
+08	43,6,4	F? 4
+09	22 C	GTO SWAB

; Remove these three lines to preserve I
; at the expense of Y.
+10	34	X-Y
+11	42 22	X-I
+12	34	X-Y

+13	43,22,B	LBL SWEND32 {B}

+14	36	ENTER
+15	36	ENTER
+16	8	[8]
+17	42 A	SL
+18	42 F	RRn
+19	43 36	LSTx
+20	42 44	WSIZE
+21	21 C	GSB SWAB
+22	34	X-Y
+23	21 C	GSB SWAB
+24	43 36	LSTx
+25	42 A	SL
+26	42 A	SL
+27	42 44	WSIZE
+28	43 36	LSTx
+29	42 B	SR
+30	43 E	RLn
+31	42 40	OR

; Remove these three lines to preserve I
; at the expense of Y.
+32	34	X-Y
+33	42 22	X-I
+34	34	X-Y

+35	43 21	RTN

+36	43,22,C	LBL SWAB {C}
+37	8	[8]
+38	42 F	RRn
+39	43 21	RTN

Decimal Shift

Sometimes when in integer mode it's useful to be able to perform decimal shifts. These two routines perform left and right decimal shifts regardless of the currently selected base.


GSB 0 - multiplies X (any base) by 10 (base 10).

001 - 43,22 0	LBL 0
002 - 1		[1]
003 - 42 11	SL
004 - 36	ENTER
005 - 42 15	RLn
006 - 43 36	LSTx
007 - 42 40	OR
008 - 20	MUL
009 - 43 21	RTN

GSB 1 - divides X (any base) by 10 (base 10).

010 - 43,22 1	LBL 1
011 - 1		[1]
012 - 21 0	GSB 0
013 - 10	DIV
014 - 43 21	RTN

NOTE: the line numbers in some of these listings have been "normalised" and are presented only for reference. There are no line-based jumps or sub-routine calls.