Tony Finch – blog

Turning a linear day number into a date is about twice as difficult as the reverse calculation, both from the computer's point of view and the programmer's. The overview of the function is to (1a) calculate the year (1b) subtract the day count for previous years (2a) calculate the month (2b) subtract the day count for previous months. These pairs of operations are essentially calculating quotients and remainders to produce successive components of the dates.

The trickiest part is dealing with the Gregorian correction. You can work out the Julian year in which a day number occurs using d*4/1461 which is a straight-forward inverse of the corresponding part of the reverse calculation (though you need to fiddle with the epoch to make the leap years fall in the right place). The equivalent for the Gregorian calendar, dividing by the average length of a year using integer arithmetic, is d*400/146097. However this produces leap year intervals of 4 or 5 years, not 4 or 8 years. Fortunately the book describes a neat trick which allows us to use this formula to estimate the day's year, then adjust it to fix any error.

Errors are detected using the forwards conversion from years to day numbers, y*1461/4 - y/100 + y/400. This formula produces the number of days from the start of 1 A.D. up to the end of year y inclusive, according to the conventional numbering. My code internally adjusts the year to start with March instead of January; the formula works for adjustments like this so long as the numbering of February doesn't change. (Note that the formula in my earlier post adjusts the year number so that it refers to the end of the previous February, to which is added a day count corresponding to the fraction of the current year that has passed.)

So what we do is calculate d*400/146097 and add one to estimate the year according to the conventional numbering. The result is always an underestimate (i.e. the formula guesses that years start up to two days late) so we may in fact get the previous year's number. We calculate the count of days up to the end of the estimated year and compare with the original day number. The original day number should normally be before the calculated end of the year. If the estimate was wrong, we need to increment the year number again.

We then need to obtain the remainder of days within this year by subtracting the day count up to the end of the previous year. To get this count we subtract one from the year and do the forward calculation again. If we combine the decrement with the estimate correction, the resulting code looks like this, where dn is day number and dy is day of year.

	y = dn*400/146097 + 1;
	y += (dn >= y*1461/4 - y/100 + y/400) - 1;
	dy = dn - y*1461/4 + y/100 - y/400;

The month calculation is more awkward than I would like because of a mismatch between the forward m*153/5 and reverse dy*5/153 rounding patterns, as I've illustrated below. In each column I've italicised the five-month cycle that corresponds to March - July.

The result is that the code needs a series of different administrative fiddles to work out the month and the day of the month. First add one month of days to line up with the italic part of the reverse pattern, then calculate the month giving March = 1 through to February = 12. Then add three months to line up with the italic part of the forward pattern, March = 4 to February = 15, and calculate the number of days before the month. Subtract this from the day of the year, then add four months of days to compensate for the fiddles, to finally get the day of the month.

	m = (dy + 31)*5/153 + 3;
	d = dy - m*153/5 + 123;

The last step is to restore the year to its usual January - December alignment. This, at least, is a simple inverse of the corresponding code in the forwards direction. There's also a preparatory adjustment of 10 months at the start of the code to set up the March - February alignment in the first place. After a few tweaks, the final code is:

void F(int d, int *py, int *pm, int *pd) {
	int y, m;
	d += 305;
	y = d*400/146097 + 1;
	y -= y*1461/4 - y/100 + y/400 > d;
	d -= y*1461/4 - y/100 + y/400 - 31;
	m = d*5/153 + 3;
	d -= m*153/5 - 92;
	if (m < 14) m -= 1;
	else m -= 13, y += 1;
	*py = y, *pm = m, *pd = d;
}

• 305 is 10 months (March - December) minus the epoch (1st January 1 A.D. is day 1)
• 400 is the number of years in a Gregorian cycle and 146097 is the number of days
• 1461 is the number of days in 4 years, usually
• 153 is the number of days in 5 consecutive months, except February
• 31 is one month of days to align the reverse month calculation, from 31*5/153 == 1
• 3 aligns the forward month calculation
• 92 is three months of days to compensate for the three months added in the previous line, from 3*153/5 plus an extra one because we count days within a month from 1 not 0
• actually 92 = 4*153/5 - 31 + 1 where the extra 31 is the same 31 we added previously

Edited to add...

It's possible to make the forward and reverse patterns line up by adjusting the 153/5 factors, which saves a couple of additions. The smallest denominator for which this works is 17, which produces the following code. Before the final fix-up, months are numbered Mar=1 - Feb=12.

void G(int d, int *py, int *pm, int *pd) {
	int y, m;
	d += 305;
	y = d*400/146097 + 1;
	y -= y*1461/4 - y/100 + y/400 > d;
	d -= y*1461/4 - y/100 + y/400 - 31;
	m = d*17/520;
	d -= m*520/17;
	if (m < 11) m += 2;
	else m -= 10, y += 1;
	*py = y, *pm = m, *pd = d;
}