Tuesday 19 May 2009

30.6001

I came across an interesting piece of code today. It calculates the Julian Date from Calendar Date...
int year, month, day;
int century, leapCycle;
double julDate;
if (month <= 2)
{
/* Align us onto a cycle which puts *\
\* February at the end of the year */
year --;
month += 12;
}
century = (int)floor(year/100);
leapCycle = 2 - century + (int)floor(century/4);
julDate = ((floor(365.25*(year + 4716))
+ floor(30.6001*(month+1)) + day) - 2451545)
+ leapCycle - 1524.5;

The idea is nice: align the start of year to 1st March to handle leap days correctly, but the constant 30.6001 just didn't look right: 0.0001 of a day isn't going to effect the answer when multiplied by a few tens of months.

Then I looked harder. The bit about calculating century and leapCycle hardly inspire confidence in the original author: I doubt that year gets converted to a float before being divided by 100 (which, incidentally, isn't a float) then converted to float then back to integer. Ditto for leapCycle - the person who wrote this was clearly having a bad day.

But then there's the precision of 30.6001 - someone who overlooked 100 vs 100.0 is hardly going to pluck this constant out of the air, so it must have a purpose. The math uses single-precision floating point which only copes with half a dozen significant digits so perhaps 30.6 gets converted to something more like 30.5999 in binary?

In fact, googling 30.6001 reveals pages of calculations for Julian Dates: this piece of code has been cut & pasted verbatim since the late 1970's

Digging a little deeper, all modern floating point implementations (software or hardware) use a standard - IEEE754 - published in 1985 which defines the layout and operations of 32-bit and 64-bit data.

Using this online calculator I found that the encoding of 30.6 is accurate to 8 significant digits using single precision 'C' float type with a standard 24-bit mantissa. You only need to worry about rounding 30.6 down to 30.5999 if you're using a much smaller 16-bit mantissa such as may well have been present before IEEE 754 standardisation.

To summarise, the magic value 30.6001 does no harm but hasn't been required for the past 20 years. Any IEEE 754 conformant floating point implementation should work just fine with 30.6 and in my own version I found that I could nicely do away with floating point altogether as 30.6 lends itself much better to fixed point numerical analysis (it's 306/10; fits in 16-bit signed integer arithmetic for any sensible month+1)

1 comment:

  1. Thx for the precision, I want to use the formulea in a script of mine, but I was not confortable to put the 30.6001 without knowing!

    Seb O.

    ReplyDelete