Wednesday, 3 June 2015

Gregorian Date Calculation Golf 3

So let's convert a Gregorian year/month/day triplet to a serial day number. To make things simple, I'm going to expression serial day numbers as Rata Die integers; i.e. the number of calendar days since the day before January 1, 1 C.E.

This means that January 1, 1 C.E. is RD=1 and June 1, 2015 is RD=735750.

Here's the classic Julian Day Number (JDN) calculation offset for the Rata Die scheme:

int ToRataDie(int year, int month, int day) {
    assert(year >= -4713);
    assert((month >= 1) && (month <= 12));
    assert((day >= 1) && (day <= 31));
    int a = (14 - month) / 12;
    int y = year + 4800 - a;
    int m = month + 12 * a - 3;
    int jdn = day + ((153 * m + 2) / 5) + 365 * y +
              (y / 4) - (y / 100) + (y / 400) - 32045;
    return jdn - 1721425;
}

Of course, for Gregorian Date Calculation Golf we want to eradicate those divisions. Fortunately, we're only interested in years between 1 and 65535 inclusive.

The calculation of "a" calls for a division by 12, but careful inspection of the logic uncovers:

    a == 1 for January and February
    a == 0 for all other months

So we can replace that division with another formula that achieves the same result, say:

    a = (18 - month) >> 4

Another culprit is the division by 5 in the calculation of the days until the start of the month:

    (153 * m + 2) / 5

By referring to the excellent Calendrical Calculations, we discover that the constants 153, 2 and 5 in the above equation are themselves somewhat arbitrary (see pp.53-54, but also Zeller's Congruence) and choices that better suit computer hardware can be used instead.

All this, along with a considerable amount of fettling, leads to the following implementation:

int ToRataDie(int year, int month, int day) {
    assert(year >= 1);
    assert((month >= 1) && (month <= 12));
    assert((day >= 1) && (day <= 31));
    int a = (18 - month) >> 4;
    int b = month + 12 * a + 1;
    int c = year - a;
    int d = c / 100;
    int e = (b * 1959) >> 6;
    int f = (c * 1461) >> 2;
    return day - d + (d >> 2) + e + f - 428;
}

Note that non-positive (non-Common Era) years are no longer supported; though this restriction can be relaxed relatively easily by offsetting the year and the result correspondingly.

We've already seen how to remove the division by 100 in our discussion of "IsLeapYear()", so our Golf solution is:

int ToRataDie(int year, int month, int day) {
    assert(year >= 1);
    assert((month >= 1) && (month <= 12));
    assert((day >= 1) && (day <= 31));
    int a = (18 - month) >> 4;
    int b = month + 12 * a + 1;
    int c = year - a;
    int d = (((c * 0x47AF) >> 16) + c) >> 7;
    int e = (b * 1959) >> 6;
    int f = (c * 1461) >> 2;
    return day - d + (d >> 2) + e + f - 428;
}

If you supply the day-of-the-year (1 to 366) instead of the month and day, the conversion is pleasingly simple:

int ToRataDie(int year, int dayOfYear) {
    assert(year >= 1);
    assert((dayOfYear >= 1) && (dayOfYear <= 366));
    int c = year - 1;
    int d = (((c * 0x47AF) >> 16) + c) >> 7;
    int f = (c * 1461) >> 2;
    return dayOfYear - d + (d >> 2) + f;
}

No comments:

Post a Comment