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;

    ReplyDelete