Computus 3

After reminiscing with Dr J R Stockton about old computers, it got to converting my 16-bit algorithm for Easter to 8-bit Z80 assembly language. I spent quite some time hand-optimising it (which is an enjoyable past-time in itself, if there's no deadline involved!) and posted the code with a short commentary.

I had thought about writing a bespoke super-optimiser to try to find the "ultimate" Z80 code, but with just under 500 bytes of machine code, that's over 10¹⁰⁰⁰ permutations. I think I might have to divide and conquer!

Modulo Arithmetic in Z80 assembler

I've done a little 8-bit assembly coding in the past, mostly Z80 (on Sinclair machines) and 6502 (on BBC Micro), but I've never really thought about arithmetic optimisations for those platforms as much as I have the last couple of months.

Imagine that you are (in the 21st century) writing Z80 date manipulation functions where the day of the week was important. You'd naturally want to perform "modulo 7" arithmetic on 16-bit quantities. If you had a routine to perform integer division by 7, you might find "z = x mod 7" given "y = floor(x / 7)" for unsigned 16-bit values as follows:

; Z80 assembly to compute 'x mod 7' from 'x' and 'x div 7'
; 'x' in HL
; 'x div 7' in DE
LD B, D       ;  4T 1B -- Take a copy of DE
LD C, E       ;  4T 1B -- BC := x div 7 
EX DE, HL     ;  4T 1B -- DE := x, HL := x div 7 
ADD HL, HL    ; 11T 1B -- HL := (x div 7) * 2
ADD HL, HL    ; 11T 1B -- HL := (x div 7) * 4
ADD HL, HL    ; 11T 1B -- HL := (x div 7) * 8
AND A         ;  4T 1B -- Clear carry flag
SBC HL, BC    ; 15T 2B -- HL := (x div 7) * 7
EX DE, HL     ;  4T 1B -- HL := x, DE := (x div 7) * 7
AND A         ;  4T 1B -- Clear carry flag
SBC HL, BC    ; 15T 2B -- HL := x - (x div 7) * 7
LD A, L       ;  4T 1B -- A:= x mod 7
; 'x mod 7' in A

At first glance, this seems fairly fast, compact code (91 T-cycles and 14 code bytes), but notice that, by the end of it, one side-effect is setting H to zero. Not very surprising as "x mod 7" is clearly an 8-bit quantity. So why go to the trouble of computing full 16-bit intermediates that have only partial bearing on the final result? Providing that the upper eight bits don't "pollute" the lower eight (e.g. during division or change of sign), you need not compute them. Especially as 16-bit arithmetic is disproportionately expensive on Z80:

• Obviously uses more precious registers,
• Addition of 16-bit quantites takes at least 11 T-cycles compared to 4 T-cycles for 8-bit quantities,
• Subtraction via 'SBC' requires worrying about the carry flag as there is no 16-bit 'SUB' instruction,
• Shifting 16-bit quantities (multiplication and division by two) other than 'ADD HL, HL' requires two, long instructions.

So the 8-bit equivalent could be:

 

; Z80 assembly to compute 'x mod 7' from 'x' and 'x div 7'
; 'x' in HL
; 'x div 7' in DE
LD A, E       ;  4T 1B -- A := x div 7 [low bits]
ADD A, A      ;  4T 1B -- A := (x div 7) * 2 [low bits]
ADD A, A      ;  4T 1B -- A := (x div 7) * 4 [low bits]
ADD A, A      ;  4T 1B -- A := (x div 7) * 8 [low bits]
SUB E         ;  4T 1B -- A := (x div 7) * 7 [low bits]
NEG           ;  8T 2B -- A := (x div 7) * -7 [low bits]
ADD A, L      ;  4T 1B -- A := x mod 7
; 'x mod 7' in A 

 

 

This is both faster and smaller (32 T-cycles and 8 code bytes) with no loss of accuracy, so I'm a bit bemused that it's taken me nearly three decades to stumble across it. And a lot of good it will do me too!

P.S. Why is 'NEG' so expensive on Z80? Given that 'CPL' is only 4 T-cycles and 2 code bytes, it seems a waste to dedicate precious silicon on a similar instruction that can be represented with existing instructions ('CPL' then 'INC A') for the same cost. Bah humbug!

Computus 2

Of course, if you're as old as I am, you'll remember programming on machines with no native division instruction (neither integer nor floating-point), so all those divisions and modulos aren't going to get you very far on a "old school" 16-bit computer. Using the information from Hacker's Delight, you can re-cast these operations to multiplications and shifts:

void easter_u16_nodiv(u16 year, u16* month, u16* day)
{
    u16 a = year >> 1;
    u16 b = year >> 2;
    u16 c = year >> 4;
    u16 d = b + c;
    u16 e = (a + d + (d >> 4) + (d >> 5) + (d >> 10)) >> 4;
    u16 f = year - e * 19;
    u16 g = f - ((f + 13) >> 5) * 19;
    d = (b >> 1) + (b >> 3);
    e = (d + (d >> 6) + (d >> 7)) >> 4;
    f = b - e * 25 + 39;
    d = e + (f >> 5);
    u16 h = (d * 3) >> 2;
    d = (d * 8) + 5;
    e = (d >> 1) + (d >> 3);
    e = (e + (e >> 6) + (e >> 7)) >> 4;
    f = d - (e * 25) + 7;
    d = (g * 19) + h - e - (f >> 5) + 15;
    e = ((d >> 1) + (d >> 5)) >> 4;
    f = d - (e * 30);
    d = f - ((f + 2) >> 5) * 30;
    g = d + ((29578 - g - (d * 32)) >> 10);
    f = b - h + g + 2;
    d = a + c + (f >> 1) + (f >> 4);
    e = (d + (d >> 6) + (d >> 12)) >> 2;
    f = year + f - (e * 7);
    e = g - f + ((f * 2341) >> 14) * 7;
    d = e >> 5;
    *day = e - d * 31;
    *month = d + 3;
}

Again, if you want to verify the algorithm (only valid for Gregorian years 0 to 65535), you can do so with Dr J R Stockton's test harness with the following Javascript snippet:

    var a, b, c, d, e, f, g, h;
    a = YR >>> 1;
    b = YR >>> 2;
    c = YR >>> 4;
    d = b + c;
    e = (a + d + (d >>> 4) + (d >>> 5) + (d >>> 10)) >>> 4;
    f = YR - e * 19;
    g = f - ((f + 13) >>> 5) * 19;
    d = (b >>> 1) + (b >>> 3);
    e = (d + (d >>> 6) + (d >>> 7)) >>> 4;
    f = b - e * 25 + 39;
    d = e + (f >>> 5);
    h = (d * 3) >>> 2;
    d = (d * 8) + 5;
    e = (d >>> 1) + (d >>> 3);
    e = (e + (e >>> 6) + (e >>> 7)) >>> 4;
    f = d - (e * 25) + 7;
    d = (g * 19) + h - e - (f >>> 5) + 15;
    e = ((d >>> 1) + (d >>> 5)) >>> 4;
    f = d - (e * 30);
    d = f - ((f + 2) >>> 5) * 30;
    g = d + ((29578 - g - (d * 32)) >>> 10);
    f = b - h + g + 2;
    d = a + c + (f >>> 1) + (f >>> 4);
    e = (d + (d >>> 6) + (d >>> 12)) >>> 2;
    f += YR - (e * 7);
    e = g - f + ((f * 2341) >>> 14) * 7;
    return e;

Computus 1

Computus is the calculation of Easter. Calendrical calculations being a bit of an obsession of mine, I've been looking at this with a view to an implementation using limited-range integer arithmetic. There's a lot of interesting stuff still being researched on this, but I haven't seen a good implementation for 16-bit unsigned arithmetic. So here's my attempt:

void easter_u16(u16 year, u16* month, u16* day)
{
    u16 a = year % 19;
    u16 b = year >> 2;
    u16 c = (b / 25) + 1;
    u16 d = (c * 3) >> 2;
    u16 e = ((a * 19) - ((c * 8 + 5) / 25) + d + 15) % 30;
    e += (29578 - a - e * 32) >> 10;
    e -= ((year % 7) + b - d + e + 2) % 7;
    d = e >> 5;
    *day = e - d * 31;
    *month = d + 3;
}

This is valid for all Gregorian years 0 to 65535. In fact, the core algorithm is valid for all non-negative years: just change the integer type to a wider, signed one.

If you want to test it using Dr J R Stockton's wonderful script tester, use the following Javascript snippet:

    var a = YR % 19;
    var b = YR >>> 2;
    var c = ((b / 25)|0) + 1;
    var d = (c * 3) >>> 2;
    var e = ((a * 19) - (((c * 8 + 5) / 25)|0) + d + 15) % 30;
    e += (29578 - a - (e << 5)) >>> 10;
    e -= ((YR % 7) + b - d + e + 2) % 7;
    return e;