[some cleanup of calendar calculations]

bk: 5cecd12f6fmazsLK8n1OFSass93pzw

diff --git a/include/ntp_calendar.h b/include/ntp_calendar.h
index 85fffbedc9b0515381c0bdec45c9595117c056e1..a13b80306f99f52b26523775054086ce1c9e0047 100644 (file)
--- a/include/ntp_calendar.h
+++ b/include/ntp_calendar.h
@@ -486,6 +486,48 @@ basedate_expand_gpsweek(unsigned short weekno);
  */
 #define        GREGORIAN_CYCLE_WEEKS (GREGORIAN_CYCLE_DAYS / 7)
 
-#define        is_leapyear(y)  (!((y) % 4) && !(!((y) % 100) && (y) % 400))
+/*
+ * Is a Greogorian calendar year a leap year? The obvious solution is to
+ * test the expression
+ *
+ * (y % 4 == 0) && ((y % 100 != 0) || (y % 400 == 0))
+ *
+ * This needs (in theory) 2 true divisions -- most compilers check the
+ * (mod 4) condition by doing a bit test. Some compilers have been
+ * even observed to partially fuse the (mod 100) and (mod 400) test,
+ * but there is an alternative formula that gives the compiler even
+ * better chances:
+ *
+ * (y % 4 == 0) && ((y % 16 == 0) || (y % 25 != 0))
+ *
+ * The order of checks is chosen so that the shorcut evaluation can fix
+ * the result as soon as possible. And the compiler has to do only one
+ * true division here -- the (mod 4) and (mod 16) can be done with
+ * direct bit tests. *If* the compiler chooses to do so.
+ *
+ * The deduction is as follows: rewrite the standard formula as
+ *  (y % 4 == 0) && ((y % 4*25 != 0) || (y % 16*25 == 0))
+ *
+ * then split the congruences:
+ *  (y % 4 == 0) && ((y % 4 != 0 || y % 25 != 0) || (y % 16 == 0 && y % 25 == 0))
+ *
+ * eliminate the 1st inner term, as it is provably false:
+ *  (y % 4 == 0) && (y % 25 != 0 || (y % 16 == 0 && y % 25 == 0))
+ *
+ * Use the distributive laws on the second major group:
+ *  (y % 4 == 0) && ((y % 25 != 0 || y % 16 == 0) && (y % 25 != 0 || y % 25 == 0))
+ *
+ * Eliminate the constant term, reorder, and voila: 
+ */
 
+static inline int
+is_leapyear(int32_t y) {
+       return !(y % 4) && (!(y % 16) || (y % 25));
+}
+/* The (mod 4) test eliminates 3/4 (or 12/16) of all values.
+ * The (mod 16) test eliminates another 1/16 of all values.
+ * 3/16 of all values reach the final division.
+ * Assuming that the true division is the most costly operation, this
+ * sequence should give most bang for the buck.
+ */
 #endif

diff --git a/libntp/ntp_calendar.c b/libntp/ntp_calendar.c
index addd50ab556ed81147e7e81b69886fae4f825cd1..19b23c82679c87bd61d1e2ac1b853cd1cecbb332 100644 (file)
--- a/libntp/ntp_calendar.c
+++ b/libntp/ntp_calendar.c
@@ -1136,6 +1136,13 @@ ntpcal_time_to_date(
  */
 
 #if !defined(HAVE_INT64)
+/* multiplication helper. Seconds in days and weeks are multiples of 128,
+ * and without that factor fit well into 16 bit. So a multiplication
+ * of 32bit by 16bit and some shifting can be used on pure 32bit machines
+ * with compilers that do not support 64bit integers.
+ *
+ * Calculate ( hi * mul * 128 ) + lo
+ */
 static vint64
 _dwjoin(
        uint16_t        mul,
@@ -1144,14 +1151,13 @@ _dwjoin(
        )
 {
        vint64          res;
-       uint32_t        p1, p2;
-       int             sf;
+       uint32_t        p1, p2, sf;
 
-       p1 = (uint32_t)hi;
-       sf = (hi < 0);
-       p1 = (p1 + sf) ^ sf;    /* absolute value if 'hi' */
+       /* get sign flag and absolute value of 'hi' in p1 */
+       sf = (uint32_t)-(hi < 0);
+       p1 = ((uint32_t)hi + sf) ^ sf;
 
-       /* assemble major units */
+       /* assemble major units: res <- |hi| * mul */
        res.D_s.lo = (p1 & 0xFFFF) * mul;
        res.D_s.hi = 0;
        p1 = (p1 >> 16) * mul;
@@ -1159,17 +1165,18 @@ _dwjoin(
        p1 = p1 << 16;
        M_ADD(res.D_s.hi, res.D_s.lo, p2, p1);
 
-       /* mul by 128, using shift */
+       /* mul by 128, using shift: res <-- res << 7 */
        res.D_s.hi = (res.D_s.hi << 7) | (res.D_s.lo >> 25);
        res.D_s.lo = (res.D_s.lo << 7);
 
-       /* fix sign */
-       if (sf)
-               M_NEG(res.D_s.hi, res.D_s.lo);
+       /* fix up sign: res <-- (res + [sf|sf]) ^ [sf|sf] */
+       M_ADD(res.D_s.hi, res.D_s.lo, sf, sf);  
+       res.D_s.lo ^= sf;
+       res.D_s.hi ^= sf;
 
-       /* properly add seconds */
+       /* properly add seconds: res <-- res + [sx(lo)|lo] */
+       p2 = (uint32_t)-(lo < 0);
        p1 = (uint32_t)lo;
-       p2 = UINT32_C(0) - (lo < 0);
        M_ADD(res.D_s.hi, res.D_s.lo, p2, p1);
        return res;
 }

diff --git a/tests/libntp/calendar.c b/tests/libntp/calendar.c
index 97778959fce81783fcf6d7ea80891e37a73da4dd..a54d9d99d6d431d09fc60f30d8a13f9d59b4b3bf 100644 (file)
--- a/tests/libntp/calendar.c
+++ b/tests/libntp/calendar.c
@@ -33,6 +33,7 @@ void  test_SplitYearDays2(void);
 void   test_RataDie1(void);
 void   test_LeapYears1(void);
 void   test_LeapYears2(void);
+void   test_LeapYears3(void);
 void   test_RoundTripDate(void);
 void   test_RoundTripYearStart(void);
 void   test_RoundTripMonthStart(void);
@@ -442,6 +443,21 @@ test_LeapYears2(void)
        return;
 }
 
+/* check the 'is_leapyear()' implementation for 4400 years */
+void
+test_LeapYears3(void)
+{
+       int32_t year;
+       int     l1, l2;
+       
+       for (year = -399; year < 4000; ++year) {
+               l1 = (year % 4 == 0) && ((year % 100 != 0) || (year % 400 == 0));
+               l2 = is_leapyear(year);
+               snprintf(mbuf, sizeof(mbuf), "y=%d", year);
+               TEST_ASSERT_EQUAL_MESSAGE(l1, l2, mbuf);
+       }
+}
+
 /* Full roundtrip from 1601-01-01 to 2400-12-31
  * checks sequence of rata die numbers and validates date output
  * (since the input is all nominal days of the calendar in that range

diff --git a/tests/libntp/run-calendar.c b/tests/libntp/run-calendar.c
index 26089200e7c0634567c156d8a237d664a3205136..31f196952651c85192b3d566adf8c1e8f219cc1e 100644 (file)
--- a/tests/libntp/run-calendar.c
+++ b/tests/libntp/run-calendar.c
@@ -41,6 +41,7 @@ extern void test_SplitYearDays2(void);
 extern void test_RataDie1(void);
 extern void test_LeapYears1(void);
 extern void test_LeapYears2(void);
+extern void test_LeapYears3(void);
 extern void test_RoundTripDate(void);
 extern void test_RoundTripYearStart(void);
 extern void test_RoundTripMonthStart(void);
@@ -89,20 +90,21 @@ int main(int argc, char *argv[])
   RUN_TEST(test_RataDie1, 33);
   RUN_TEST(test_LeapYears1, 34);
   RUN_TEST(test_LeapYears2, 35);
-  RUN_TEST(test_RoundTripDate, 36);
-  RUN_TEST(test_RoundTripYearStart, 37);
-  RUN_TEST(test_RoundTripMonthStart, 38);
-  RUN_TEST(test_RoundTripWeekStart, 39);
-  RUN_TEST(test_RoundTripDayStart, 40);
-  RUN_TEST(test_IsoCalYearsToWeeks, 41);
-  RUN_TEST(test_IsoCalWeeksToYearStart, 42);
-  RUN_TEST(test_IsoCalWeeksToYearEnd, 43);
-  RUN_TEST(test_DaySecToDate, 44);
-  RUN_TEST(test_GpsRollOver, 45);
-  RUN_TEST(test_GpsRemapFunny, 46);
-  RUN_TEST(test_GpsNtpFixpoints, 48);
-  RUN_TEST(test_NtpToNtp, 49);
-  RUN_TEST(test_NtpToTime, 50);
+  RUN_TEST(test_LeapYears3, 36);
+  RUN_TEST(test_RoundTripDate, 37);
+  RUN_TEST(test_RoundTripYearStart, 38);
+  RUN_TEST(test_RoundTripMonthStart, 39);
+  RUN_TEST(test_RoundTripWeekStart, 40);
+  RUN_TEST(test_RoundTripDayStart, 41);
+  RUN_TEST(test_IsoCalYearsToWeeks, 42);
+  RUN_TEST(test_IsoCalWeeksToYearStart, 43);
+  RUN_TEST(test_IsoCalWeeksToYearEnd, 44);
+  RUN_TEST(test_DaySecToDate, 45);
+  RUN_TEST(test_GpsRollOver, 46);
+  RUN_TEST(test_GpsRemapFunny, 47);
+  RUN_TEST(test_GpsNtpFixpoints, 49);
+  RUN_TEST(test_NtpToNtp, 50);
+  RUN_TEST(test_NtpToTime, 51);
 
   return (UnityEnd());
 }