The HP-41 calculators can be programmed with microcode, aka MCODE. This is how their normal instructions are implemented. The emulator even has an MCODE console where  you can step though its instructions. On the physical machine running your own MCODE requires making ROMs (or something that pretends to be a ROM…). I am not familiar with how this works with the emulator, but I have found an SDK.

Microcode is of course very low level so it could be interesting to compile e.g. C to MCODE. This is also so utterly useless that it would be insanely fun to at least make a tiny prototype or a starting point for such a compiler. Imagine creating e.g. MCODE assembly with clang/llvm:

./clang.exe --target=hp41 -S hello.c

If you have the PPC ROM synthesizing instructions is much easier, and you don’t have to “jailbreak” the calculator. I believe assemblers/compilers that also work with the emulator makes this even easier. I have actually not tested this, but I do think they support synthetic instructions.

One of the most fun things to do on an HP-41 is to start utilizing instructions you cannot key in but that still are actually functioning codes. Slicing bytecodes together to make these instructions is called “Synthetic Programming”, as they will have to be “synthesized” in some more or less artificial manner than using the keyboard. I have a couple of books on this, these and others are now freely available for download at HP Calculator Literature.

There is actually an online emulator for my very first programmable calculator, the TI-51 III, aka TI-55. Just start pushing buttons! Even the on/off switch works.

The TI screenshots below were downloaded from the “TI-SmartView CE-T Emulator Software”, as getting screenshots from the calculator itself does not work regardless of connection method.

The HP screenshots were taken in the “V41” virtual HP-41. Programs were also downloaded from the emulator and decompiled using “HP41UC“, an ancient 16-bit HP-41 compiler/decompiler that today must be run in DOSBox.

It was much more fun to program this on my old HP calculator (not necessarily an optimal solution, but I wrote this around 40 years ago, my very first version…):

LBL "HANOI"
"N?"
PROMPT
STO 01
"A"
ASTO 02
"B"
ASTO 03
"C"
ASTO 04
10
STO 00
XEQ "INFIX"
STOP
LBL "INFIX"
1
RCL 01
X=Y?
GTO 01
XEQ "PUSH"
RCL 02
XEQ "PUSH"
RCL 03
XEQ "PUSH"
RCL 04
XEQ "PUSH"
1
ST- 01
RCL 04
X<> 03
STO 04
XEQ "INFIX"
XEQ "POP"
STO 04
XEQ "POP"
STO 03
XEQ "POP"
STO 02
XEQ "POP"
STO 01
LBL 01
CLA
ARCL 02
>"-"
ARCL 03
AVIEW
1
RCL 01
X=Y?
RTN
XEQ "PUSH"
RCL 02
XEQ "PUSH"
RCL 03
XEQ "PUSH"
RCL 04
XEQ "PUSH"
1
ST- 01
RCL 04
X<> 02
STO 04
XEQ "INFIX"
XEQ "POP"
STO 04
XEQ "POP"
STO 03
XEQ "POP"
STO 02
XEQ "POP"
STO 01
END

LBL "PUSH"
STO IND 00
1
ST+ 00
RTN
LBL "POP"
1
ST- 00
RCL IND 00
END

I called the pegs “A”, “B” and “C” this time:

The Tower of Hanoi is something I “have to” do on any calculator I encounter, but using Python feels almost like cheating, it is just too easy. Unfortunately there is little else of interest for this task on this calculator.

def hanoi(n,a,b,c):
  if n==1:
    print(a," >> ",b)
  else:
   hanoi(n-1,a,c,b)
   hanoi(1,a,b,c)
   hanoi(n-1,c,b,a)

hanoi(3,1,2,3)

To fit it on one screen we are just moving 3 disks from peg 1 to peg 2, using peg 3 as a helper:

Here you can see the program running:

The numbers can of course be complex, so let us try a different input:

Believe it or not, even with their strict limitations, the fairly simplistic functions from the last post are sufficient building blocks to demonstrate a simple quantum circuit. Now we just have to declare the necessary quantum gates and initialize the three input qubits and we are ready to actually “build” the circuit.

This is the matrix definition of the quantum gates needed (see the drawing in the first post, “I” is needed to expand the circuit to more bits as all the gates have either a one or two bit input):

X=[[0,1],[1,0]]
Z=[[1,0],[0,-1]]
I=[[1,0],[0,1]]
H=[[1/sqrt(2),1/sqrt(2)],[1/sqrt(2),-1/sqrt(2)]]
CNOT=[[1,0,0,0],[0,1,0,0],[0,0,0,1],[0,0,1,0]]

The usual suspects (aka “Alice” and “Bob”) can now start doing their qubit dance routine. Let us give Alice a nice little qubit, then initialize the inputs with it (other inputs are zeroed):

#Alice qubit
q=[[-1/sqrt(2)],[-1/sqrt(2)]]
print("Alice qubit:",q)

#Initialize inputs
q=[q[0],[0],[0],[0],q[1],[0],[0],[0]]

Creating the circuit:

#Build circuit
U=kron(I,kron(H,I))
U=mult(kron(I,CNOT),U)
U=mult(kron(CNOT,I),U)
U=mult(kron(H,kron(I,I)),U)

Running the circuit on the three inputs:

#Apply circuit
q=mult(U,q)

Alice now measures her qubits, the result will be made available to Bob:

#Alice measurements
b0=qbit(0,M(q))
q=qbitset(0,b0,q)
b1=qbit(1,M(q))
q=qbitset(1,b1,q)
print("Alice measurements:",b0,b1)

Bob now, based on the measurements from Alice, can do some operations on his entangled qubit (using a small one qubit input circuit, finally turning his qubit into something matching the one Alice started with):

#Bob qubit
r=[[0],[0]]
for i in range(len(q)):
  if i&1==1:
    r[1][0]+=q[i][0]
  else:
    r[0][0]+=q[i][0]

#Conditionally apply X and/or Z
if b1==1:
  r=mult(X,r)
if b0==1:
  r=mult(Z,r)
print("Bob qubit:",r)

We will need a few functions to do necessary mathematical operations on qubits (or more specifically on the state vector or probability vector of a quantum system). Note that I follow what I believe is the most common numbering of inputs with qubit 0 being the first input from above in the circuit and the first bit from the left in the state vector index, i.e the opposite of how we usually number bits in binary, from right to left (e.g if probability[1] in a three qubit system is 100% then the bits in the binary representation are 001 with a 1 in bit number 0, which is qubit number 2). You can see this reversal when we are doing actual bit operations (bit shifts) in the following functions.

We need a matrix multiplication routine (dot product):

def mult(a,b):
  c=[[0 for x in range(len(b[0]))] for x in range(len(a))]
  for i in range(len(a)):
    for j in range(len(b[i])):
      for k in range(len(b)):
        c[i][j]+=a[i][k]*b[k][j]
  return c

We need a matrix tensor product routine (cross product or more specifically the Kronecker product):

def kron(a,b):
  c=[[0 for x in range(len(a[0])*len(b[0]))]\
     for x in range(len(a)*len(b))]
  for i in range(len(a)):
    for k in range(len(b)):
      for j in range(len(a[i])):
        for l in range(len(b[k])):
          c[i*len(b)+k][j*len(b[k])+l]=a[i][j]*b[k][l]
  return c

We need a routine for measuring qubits (or more specifically making a measurement on the probability vector, also note that this function for simplicity measures ALL qubits in the system, but that is ok as the state is not stored back, instead we can set wanted measurements back with a qbitset function, see later):

def M(a):
  c=[[0] for x in range(len(a))]
  r=random()
  j=len(a)-1
  s=0.0
  for i in range(j):
    s+=abs(a[i][0])**2
    if s>=r:
      j=i
      break
  c[j][0]=1
  return c

We need a routine for storing the measured state of a qubit (0 or 1) back into the probability vector (by zeroing out states we now know have probability zero, after which a normalization of the probabilities are needed, see later):

def qbitset(p,v,a):
  for i in range(len(a)):
    if (i>>round(log(len(a),2)-p-1))&1!=v:
      a[i]=[0]
  return norm(a)

We need a routine for normalizing a probability vector (the sum of the probabilities must be 1):

def norm(a):
  s=0.0
  for i in range(len(a)):
    s+=abs(a[i][0])**2
  t=sqrt(s)
  for i in range(len(a)):
    a[i][0]/=t
  return a

We need a function for returning the value of a qubit (0 or 1) after a measurement (by finding the state that now has probability 1, otherwise if the qubit is still in superposition just return -1) :

def qbit(p,a):
  for i in range(len(a)):
    if a[i][0]==1:
      if (i>>round(log(len(a),2)-p-1))&1==1:
        return 1
      else:
        return 0
  return -1