MMIXdb — a Tiptoe through the Tangents

Meet mmixdb, part of checksmix — a Unix-style debugger for MMIX. It lets you approach your program from the point of view of the machine state and that's how MMIX is meant to work.

MMIX

MMIX is Knuth's 64-bit RISC machine, the one that replaced MIX throughout the current edition of The Art of Computer Programming. It has 256 general-purpose registers, a byte-addressed 64-bit space, and an instruction set which is easy to intuit.

checksmix is a Rust implementation that assembles and interprets MMIX.

How it works

Every MMIX instruction is exactly four bytes: an opcode and three operand bytes, OP X,Y,Z. Subroutines are called through PUSHJ, which moves the register window forward so the callee sees its arguments as $0, $1, …, and POP, moves it back and drops the return values into the caller's hole.

p $1 is how you'd normally check a register.

Registers take names through IS:

Sum     IS      $1
Den     IS      $2

mmixdb reads the same symbol table as the assembler. Once Sum IS $1 is declared, p Sum works.

Defining Pi

Pi has no closed form in arithmetic operations that assembly language provides. Pi is computed from a series expansion at the accuracy desired or required. The shortest route to pi from integer arithmetic is based on the tangent. tan(π/4) = 1, and atan(1) = π/4, thus:

π = 4 · atan(1)

Gregory's series expands the arctangent as

atan(x) = x - x³/3 + x⁵/5 - x⁷/7 + …

At x = 1 every power collapses to 1 and the series becomes 1 - 1/3 + 1/5 - 1/7 + …, the alternating reciprocals of the odd numbers. The error at x = 1 after N terms works out to roughly 1/N. The whole computation is a single loop over arithmetic operations. For about four digits of pi we need roughly 100,000 iterations.

Tangents

Tangent and arctangent aren't assembler instructions — MMIX gives you arithmetic. Every transcendental function that is required, you build yourself. In our case, that means a small loop.

The program

pi_series.mms does the arithmetic in fixed point, scaled by 10⁹. It isn't the best way to do this numerically, but it's easy to play with and see which term is which:

Scale   IS      1000000000              % fixed-point scale, 10^9
Terms   IS      100000                  % series terms to sum

Sum     IS      $1                      % running total, scaled
Den     IS      $2                      % denominator: 1, 3, 5, 7, ...
Term    IS      $3                      % Scale/Den
Neg     IS      $4                      % 0 = add this term, 1 = subtract
Cnt     IS      $5                      % terms remaining
One     IS      $6                      % Scale, in a register

Main    SETI    Sum,0
        SETI    Den,1
        SETI    Neg,0
        SETI    One,Scale
        SETI    Cnt,Terms

Loop    DIV     Term,One,Den            % Term = Scale/Den
        BNZ     Neg,Minus
        ADDU    Sum,Sum,Term
        JMP     Step
Minus   SUBU    Sum,Sum,Term
Step    XOR     Neg,Neg,1               % alternate the sign
        ADDU    Den,Den,1               % next denominator
        SUBU    Cnt,Cnt,1
        BNZ     Cnt,Loop

        MULU    Sum,Sum,4               % pi = 4 * atan(1)

The rest of the file is a LOC preamble, an output buffer, and a dozen lines that turn the scaled integer into decimal digits.

$ checksmix pi_series.mms
pi ~= 2.772568160

Definitely not right. Where's the bug?

Is it noise? 2.772568/4=0.693142. It's not a coincidence. ln 2 is not what we want. Mathematicians already know the bug. Full disclosure this bug was made up before I wrote the article. It's a fun quirk!

Finding it

b Loop breaks at the top of the loop body, so run stops on the first term and each c after that advances exactly one term. s steps a single instruction and follows a PUSHJ into the callee; n runs the call to completion and stops after it. Inside this loop there is nothing to step over, so both do the same thing.

Term 1 checks out. p Den is 1, and one s past the DIV gives p Term = 1000000000 — one, scaled.

Term 2 does not. p Den is 2 where it should be 3, and p Term is 500000000: one half.

That narrows it to the increment. b 37 stops on it and l prints the surrounding lines with the offending line:

> 37            ADDU    Den,Den,1               % next denominator

Add one to the denominator, or 1, 2, 3, 4..., but that is not what we wanted. We wanted 1, 3, 5.... This is 1 - 1/2 + 1/3 - 1/4 + … which gives the alternating harmonic series — ln 2.

The fix

We should add 2 here instead.

ADDU    Den,Den,2               % next odd denominator

$ checksmix pi_series.mms
pi ~= 3.141582640

Four correct decimal places for 100,000 terms. Pi is 3.141592654 — fixed. 🕷️