Wednesday, 22 February 2017

Computus 4

I noticed a few days ago that Dr J R Stockton's web pages on Computus are now longer available. The last impression of the site in the Wayback Machine was 2015-09-07. In it he lists a number of JavaScript algorithms for the calculation of the date of Easter Sunday in the Gregorian calendar. The result of the following functions are "the day of March" such that 1 indicates March 1st, 32 indicates April 1st, etc.

The first version is my own humble attempt (here translated to C/C++):

int EasterDayOfMarch(int year) {
// 25 ops = 2 div, 3 rem, 2 mul, 5 shift, 8 add, 5 sub
int a = year % 19;
int b = year >> 2;
int c = (b / 25) + 1;
int d = (c * 3) >> 2;
int e = ((c << 3) + 5) / 25;
int f = (a * 19 + d - e + 15) % 30;
int g = f + ((29578 - a - (f << 5)) >> 10);
int h = g - ((year + b - d + g + 2) % 7);
return h;
}

This is derived from Gauss's algorithm and subsequent revisions. However, the calculation of 'g' is my own work and chosen such that it can be computed on a 16-bit machine efficiently (though it works for all integer sizes).

There are other minor variations outlined on Dr Stockton's page:

// Lichtenberg
int g = 28 + f - ((f + a / 11) / 29);

and

// Hutchins
int g = 28 + f - ((a + f * 11) / 319);

There is also a C implementation which is the equivalent of:

// Pemberton
int g = 28 + f - (int)((f + g + (int)(a > 10)) > 28);

As mentioned previously, Al Petrofsky came up with another variation on the theme:

int EasterDayOfMarchPetrofsky(int year) {
// 21 ops = 4 div, 3 rem, 4 mul, 1 shift, 6 add, 3 sub
int a = (year % 19) * 6060;
int b = year >> 2;
int c = b / 25;
int p = (c * 2267) - ((b / 100) * 6775) + 3411;
int q = ((a + ((p / 25) * 319) - 1) % 9570) / 330;
int r = 28 + q - (year + b + p + q) % 7;
return r;
}

This cleverly uses fewer operations, though requires integer sizes greater than sixteen bits for larger years (not a great imposition, these days!)

To convert a "day of March" 'h' to an actual date, the following code can be used:

void Easter(int year, int* month, int* day) {
int a = year % 19;
int b = year >> 2;
int c = (b / 25) + 1;
int d = (c * 3) >> 2;
int e = ((c << 3) + 5) / 25;
int f = (a * 19 + d - e + 15) % 30;
int g = f + ((29578 - a - (f << 5)) >> 10);
int h = g - ((year + b - d + g + 2) % 7);
int i = h >> 5; // is it April?
*month = i + 3;
*day = (h & 31) + i;
}

1 comment :

1. I've just noticed that you can possibly improve "EasterDayOfMarchPetrofsky" by changing one of the divisions to a shift:

int p = (c * 2267) - ((c >> 2) * 6775) + 3411;