/***********************************************************************

 Julian Day Number: for the Gregorian calendar

 Numbering the days of the year from 1 to 365 (or 366 for leapyears).

 Reference:
 - J. Duglas Robertson;
   Tableless Date Conversion (Remark on Algorithm 398),
   Comminications of the ACM, 15 (Oct. 1972) 918
 - Hatcher, D. A., Simple Formulae for Julian Day Numbers and Calendar Dates,
   Quarterly Journal of the Royal Astronomical Society, 25 (1984) 53-55.
 - Hatcher, D. A., Generalized Equations for Julian Day Numbers and Calendar
   Dates, Quarterly Journal of the Royal Astronomical Society, 26 (1985)
   151-155. Includes coefficients for various calendar systems,
   including the Egyptian, Alexandrian, Roman, Gregorian, and Islamic.
 - M. Keith & T. Craver; 
   The ultimate perpetual calendar?, JoRM 22:4 (1990) 280-282

 Annotation:
   Julian Day Number = day + monthShift(month) - februaryCorrection(month,year)

   The monthShift function:
   month      |   1   2   3   4   5   6   7   8   9  10  11  12
   monthLength|  31  30* 31  30  31  30  31  31  30  31  30  31
   monthShift |   0  31  60  92 122 153 183 214 245 275 306 336

   The februaryCorrection function:
   month      |   1   2   3   4   5   6   7   8   9  10  11  12
   correction |   0   0   1   1   1   1   1   1   1   1   1   1   leapyear
   correction |   0   0   2   2   2   2   2   2   2   2   2   2   normal year

  Method: J. Duglas Robertson
   monthShift = (611*(month+2))/20 - 91;   *** integer arithmetic ***
   The multiplication by 611 = 1001100011_2 is cheap.
   But the division by 20 is expensive.

  Method: Torsten Sillke
   monthShift = 30*month + floor(gamma*month) - 30;
   with 5/9 <= gamma < 4/7.
   A good choise is: gamma = (5+4)/(9+7) = 9/16.
   This is most the most robust using floating point arithmetic.
   The fraction 9/16 is very cheap in integer arithmetic.
   Then we have floor(9/16*month) = floor((month + floor(month/8))/2)
   (M. Keith & T. Craver use gamma=5/9 for the monthShift function.)

   Rounding to zero method.
   monthShift = 31*month - (month-8)/2 - 34;   *** integer arithmetic ***
   Note that the division is not an arithmetical right shift.

  The februaryCorrection function:
   This is the most expansive part in the computation,
   as the leapyear rule needs a real division.
-----------------------------------------------------------------------
   mailto:Torsten.Sillke@uni-bielefeld.de    1999-03-03
***********************************************************************/

int leapyear(int year)
{
   if (year %   4) return 0;
   if (year % 100) return 1;
   if (year % 400) return 0;
   return 1;
}

int leapyear0(int year)
{
   return year%4 == 0 && ( year%100 != 0 || year%400 == 0 );
}

/*** different implementations of the juliandaynumber function ***/

int juliandaynumber(int day, int month, int year)
{
   int j = day + (611*(month+2))/20 - 91;
   /* So far february counts 30 days. Now correct this. */
   if (month >= 3)
      j += leapyear(year) - 2;
   return j;
}


int juliandaynumber0(int day, int month, int year)
{
   int j = day + 30*month + (month + month/8)/2 - 30;
   /* So far february counts 30 days. Now correct this. */
   if (month >= 3)
      j += leapyear(year) - 2;
   return j;
}


int juliandaynumber1(int day, int month, int year)
{
   int j = day + 31*month - (month-8)/2 - 34;
   /* So far february counts 30 days. Now correct this. */
   if (month >= 3)
      j += leapyear(year) - 2;
   return j;
}


/*** T E S T ********************************************************/

int main()
{
   int d, m, y;
#if 0
   scanf("%d%d%d", &d, &m, &y);
   printf("Julian Day Number of (%d,%d,%d) is %d\n", 
	  d, m, y, juliandaynumber(d,m,y));
#else
   d = 0; y = 1;
   for (m=1;m<=12;m++)
     printf("Julian Day Shift %3d %3d \n", 
	  30*m + (m + m/8)/2 - 30,
	  30*m + 5*m/9 - 30
     );
#endif

   return 0;
}