Wednesday 3 June 2015

Gregorian Date Calculation Golf 4

Calculating a Gregorian year/month/day triplet from a Rata Die integer involves finding our temporal position within a number of cycles:

  • 146097 days for the 400-year Gregorian leap year rule cycle,
  • 1461 days for the 4-year Julian leap year rule cycle,
  • 12 months in a year, and
  • 30.6 days in the days-in-a-month approximation (see here)
This generally means a lot of divisions and modular arithmetic, which causes no end of problems for Gregorian Date Calculation Golf, particularly Rule 8.

Let's start with the canonical Richards' algorithm:


void FromRataDie(int rd, int& year, int& month, int& day) {
    int jdn = rd + 1721425;
    assert(jdn >= 0);
    int f = jdn + 1401 + (((4 * jdn + 274277) / 146097) * 3) / 4 - 38;
    int e = 4 * f + 3;
    int g = (e % 1461) / 4;
    int h = 5 * g + 2;
    day = (h % 153) / 5 + 1;
    month = (h / 153 + 2) % 12 + 1;
    year = (e / 1461) - 4716 + (14 - month) / 12;
    assert(year >= -4713);
    assert((month >= 1) && (month <= 12));
    assert((day >= 1) && (day <= 31));
}

We've already seen how we can get rid of most of these divisions using bitwise shifts or by using different divisors entirely (instead of 153 and 12, for instance):

void FromRataDie(int rd, int& year, int& month, int& day) {
    assert(rd >= 1);
    int a = rd << 2;
    int b = (a + 147321) / 146097;
    int c = ((a + b * 3) | 3) + 1220;
    int d = c / 1461;
    int e = c - d * 1461;
    int f = (e >> 2) * 2141 + 197688;
    int g = f >> 16;
    int h = (g + 3) >> 4;
    day = ((f & 0xFFFF) / 2141) + 1;
    month = g - h * 12;
    year = d + h;
    assert(year >= 1);
    assert((month >= 1) && (month <= 12));
    assert((day >= 1) && (day <= 31));
}

The rather curious bitwise-or with three in the calculation of "c" is due to the observation:

    (x | 3) == (((x / 4) * 4) + 3)


The new constant 2141 is due to the following approximation:

    153 ÷ 5 ~= 65536 ÷ 2141

Alas, our "Division by Integer Constants" trick must resort to 64-bit quantities to obtain a true Golf solution:

void FromRataDie(int rd, int& year, int& month, int& day) {
    assert((rd >= 1) && (rd <= 23936166));
    uint64_t a = uint64_t(rd << 2);
    uint64_t b = (((a + 147321) * 0xE5AC1AF3) >> 49);
    uint64_t c = ((a + b * 3) | 3) + 1220;
    int d = int((c * 0xB36D8398) >> 42);
    int e = int(c - d * 1461);
    int f = (e >> 2) * 2141 + 197688;
    int g = f >> 16;
    int h = (g + 3) >> 4;
    day = (((f & 0xFFFF) * 0x7A71) >> 26) + 1;
    month = g - h * 12;
    year = d + h;
    assert((year >= 1) && (year <= 65535));
    assert((month >= 1) && (month <= 12));
    assert((day >= 1) && (day <= 31));
}

One interesting thing about this algorithm is that none of the constants (with the possible exceptions of 12 and 1461 = 365.25 * 4) look like they've got anything to do with date calculations. Indeed, it looks more like a pseudo-random number generator.

Computing the year and the day-of-year (1 to 366) from the Rata Die integer is left as an exercise for the reader.

No comments:

Post a Comment