Computus: Calculating Easter Sunday

Computus is the algorithm used to calculate the date of Easter. Easter Sunday falls on the first Sunday after the full moon following the Spring equinox. For the purpose of the calculation the full moon is defined as day 14 of the lunar month and the Spring equinox as 20th March.

Unfortunately Computus defies any attempt to render it with beautiful code! This C function roughly follows the assembly language so is a little uglier than strictly necessary:

easter(year, month, date)
int year, *month, *date;
{
    int gold,cent,sa,la,epact,a18,da;
    gold = year%19;
    cent = year/100;
    sa = cent-cent/4;
    la = (8*cent+13)/25;
    epact = (19*gold-sa-la+15)%30;
    a18 = (gold+11*epact)/319;
    da = ((cent%4+year%100/4)*2+a18+32-year%4-epact)%7;
    *month = (90+epact+da-a18)/25;
    *date = (*month+19+epact+da-a18)%32;
}

The algorithm is similar to MaybeGauss1 found in J R Stockton's collection of algorithms for Easter Sunday and is valid for the Gregorian calendar well into the fourth millenium. The algorithm can be adapted to calculate a number of other dates:

• Shrove Tuesday - 47 days before Easter Sunday
• First Sunday in Lent - 42 days before
• Palm Sunday - 7 days before
• Whit Sunday - 49 days after

Finally here's the same algorithm in 8086 assembly language, length 128 bytes. On entry, AX is the year. On exit AL is the day, AH is the month:

easter:
  push cx
  push dx
  push bx
  push bp
  push si
  push di

  mov bp,ax     ; bp = year (1583:3999) 

  mov cx,100
  cwd
  div cx
  push dx
  xchg si,ax    ; si = century - 1

  mov ax,bp
  mov cl,19
  cwd           
  div cx
  mov bx,dx     ; bx = golden number - 1
  xchg ax,dx

  mul cl
  add ax,15
                ; ax in range (15:357)
  mov dx,si
  add ax,si
  shr dx,1
  shr dx,1
  sub ax,dx
  push ax
  mov ax,8
  mul si
  add ax,13
  mov cl,25
  div cx                      
  xchg dx,ax
  pop ax
  sub ax,dx
  mov cl,30
  cwd
  div cx
  mov di,dx

  mov al,11
  mul dx
  add ax,bx
  mov cl,206    ; multiply by 206 and discard the
  mul cx        ; lower 16 bits of the result.
                ; shorter than dividing by 319

  sub di,dx
  xchg ax,si
  and al,3 
  pop dx
  shr dx,1
  shr dx,1
  add ax,dx
  shl ax,1
  and bp,3
  lea bp,[bp+di-32]
  sub ax,bp
  mov cl,7
  cwd
  div cx
  xchg ax,dx
  add ax,di
  mov bp,ax
  add al,90
  mov cl,25
  div cl
  mov ah,al
  add al,19
  add ax,bp
  and al,31

  pop di
  pop si
  pop bp
  pop bx
  pop dx
  pop cx
  ret