2 * ntp_calendar.c - calendar and helper functions
4 * Written by Juergen Perlinger (perlinger@ntp.org) for the NTP project.
5 * The contents of 'html/copyright.html' apply.
7 * --------------------------------------------------------------------
8 * Some notes on the implementation:
10 * Calendar algorithms thrive on the division operation, which is one of
11 * the slowest numerical operations in any CPU. What saves us here from
12 * abysmal performance is the fact that all divisions are divisions by
13 * constant numbers, and most compilers can do this by a multiplication
14 * operation. But this might not work when using the div/ldiv/lldiv
15 * function family, because many compilers are not able to do inline
16 * expansion of the code with following optimisation for the
17 * constant-divider case.
19 * Also div/ldiv/lldiv are defined in terms of int/long/longlong, which
20 * are inherently target dependent. Nothing that could not be cured with
21 * autoconf, but still a mess...
23 * Furthermore, we need floor division in many places. C either leaves
24 * the division behaviour undefined (< C99) or demands truncation to
25 * zero (>= C99), so additional steps are required to make sure the
26 * algorithms work. The {l,ll}div function family is requested to
27 * truncate towards zero, which is also the wrong direction for our
30 * For all this, all divisions by constant are coded manually, even when
31 * there is a joined div/mod operation: The optimiser should sort that
32 * out, if possible. Most of the calculations are done with unsigned
33 * types, explicitely using two's complement arithmetics where
34 * necessary. This minimises the dependecies to compiler and target,
35 * while still giving reasonable to good performance.
37 * The implementation uses a few tricks that exploit properties of the
38 * two's complement: Floor division on negative dividents can be
39 * executed by using the one's complement of the divident. One's
40 * complement can be easily created using XOR and a mask.
42 * Finally, check for overflow conditions is minimal. There are only two
43 * calculation steps in the whole calendar that suffer from an internal
44 * overflow, and these conditions are checked: errno is set to EDOM and
45 * the results are clamped/saturated in this case. All other functions
46 * do not suffer from internal overflow and simply return the result
47 * truncated to 32 bits.
49 * This is a sacrifice made for execution speed. Since a 32-bit day
50 * counter covers +/- 5,879,610 years and the clamp limits the effective
51 * range to +/-2.9 million years, this should not pose a problem here.
56 #include <sys/types.h>
58 #include "ntp_types.h"
59 #include "ntp_calendar.h"
60 #include "ntp_stdlib.h"
62 #include "ntp_unixtime.h"
64 /* For now, let's take the conservative approach: if the target property
65 * macros are not defined, check a few well-known compiler/architecture
66 * settings. Default is to assume that the representation of signed
67 * integers is unknown and shift-arithmetic-right is not available.
69 #ifndef TARGET_HAS_2CPL
70 # if defined(__GNUC__)
71 # if defined(__i386__) || defined(__x86_64__) || defined(__arm__)
72 # define TARGET_HAS_2CPL 1
74 # define TARGET_HAS_2CPL 0
76 # elif defined(_MSC_VER)
77 # if defined(_M_IX86) || defined(_M_X64) || defined(_M_ARM)
78 # define TARGET_HAS_2CPL 1
80 # define TARGET_HAS_2CPL 0
83 # define TARGET_HAS_2CPL 0
87 #ifndef TARGET_HAS_SAR
88 # define TARGET_HAS_SAR 0
92 *---------------------------------------------------------------------
93 * replacing the 'time()' function
94 *---------------------------------------------------------------------
97 static systime_func_ptr systime_func = &time;
98 static inline time_t now(void);
103 systime_func_ptr nfunc
106 systime_func_ptr res;
111 systime_func = nfunc;
120 return (*systime_func)(NULL);
124 *---------------------------------------------------------------------
125 * Get sign extension mask and unsigned 2cpl rep for a signed integer
126 *---------------------------------------------------------------------
129 static inline uint32_t
133 # if TARGET_HAS_2CPL && TARGET_HAS_SAR && SIZEOF_INT >= 4
135 /* Let's assume that shift is the fastest way to get the sign
136 * extension of of a signed integer. This might not always be
137 * true, though -- On 8bit CPUs or machines without barrel
138 * shifter this will kill the performance. So we make sure
139 * we do this only if 'int' has at least 4 bytes.
141 return (uint32_t)(v >> 31);
145 /* This should be a rather generic approach for getting a sign
148 return UINT32_C(0) - (uint32_t)(v < 0);
153 static inline uint32_t
154 int32_to_uint32_2cpl(
161 /* Just copy through the 32 bits from the signed value if we're
162 * on a two's complement target.
168 /* Convert from signed int to unsigned int two's complement. Do
169 * not make any assumptions about the representation of signed
170 * integers, but make sure signed integer overflow cannot happen
171 * here. A compiler on a two's complement target *might* find
172 * out that this is just a complicated cast (as above), but your
173 * mileage might vary.
176 vu = ~(uint32_t)(-(v + 1));
185 static inline int32_t
186 uint32_2cpl_to_int32(
193 /* Just copy through the 32 bits from the unsigned value if
194 * we're on a two's complement target.
200 /* Convert to signed integer, making sure signed integer
201 * overflow cannot happen. Again, the optimiser might or might
202 * not find out that this is just a copy of 32 bits on a target
203 * with two's complement representation for signed integers.
206 v = -(int32_t)(~vu) - 1;
215 /* Some of the calculations need to multiply the input by 4 before doing
216 * a division. This can cause overflow and strange results. Therefore we
217 * clamp / saturate the input operand. And since we do the calculations
218 * in unsigned int with an extra sign flag/mask, we only loose one bit
219 * of the input value range.
221 static inline uint32_t
226 static const uint32_t limit = UINT32_MAX/4u;
227 if ((mu ^ vu) > limit) {
235 *---------------------------------------------------------------------
236 * Convert between 'time_t' and 'vint64'
237 *---------------------------------------------------------------------
249 # if SIZEOF_TIME_T <= 4
253 res.D_s.lo = (uint32_t)-tt;
254 M_NEG(res.D_s.hi, res.D_s.lo);
256 res.D_s.lo = (uint32_t)tt;
259 # elif defined(HAVE_INT64)
265 * shifting negative signed quantities is compiler-dependent, so
266 * we better avoid it and do it all manually. And shifting more
267 * than the width of a quantity is undefined. Also a don't do!
271 res.D_s.lo = (uint32_t)tt;
272 res.D_s.hi = (uint32_t)(tt >> 32);
273 M_NEG(res.D_s.hi, res.D_s.lo);
275 res.D_s.lo = (uint32_t)tt;
276 res.D_s.hi = (uint32_t)(tt >> 32);
292 # if SIZEOF_TIME_T <= 4
294 res = (time_t)tv->D_s.lo;
296 # elif defined(HAVE_INT64)
298 res = (time_t)tv->q_s;
302 res = ((time_t)tv->d_s.hi << 32) | tv->D_s.lo;
310 *---------------------------------------------------------------------
311 * Get the build date & time
312 *---------------------------------------------------------------------
315 ntpcal_get_build_date(
319 /* The C standard tells us the format of '__DATE__':
321 * __DATE__ The date of translation of the preprocessing
322 * translation unit: a character string literal of the form "Mmm
323 * dd yyyy", where the names of the months are the same as those
324 * generated by the asctime function, and the first character of
325 * dd is a space character if the value is less than 10. If the
326 * date of translation is not available, an
327 * implementation-defined valid date shall be supplied.
329 * __TIME__ The time of translation of the preprocessing
330 * translation unit: a character string literal of the form
331 * "hh:mm:ss" as in the time generated by the asctime
332 * function. If the time of translation is not available, an
333 * implementation-defined valid time shall be supplied.
335 * Note that MSVC declares DATE and TIME to be in the local time
336 * zone, while neither the C standard nor the GCC docs make any
337 * statement about this. As a result, we may be +/-12hrs off
338 * UTC. But for practical purposes, this should not be a
343 static const char build[] = MKREPRO_TIME "/" MKREPRO_DATE;
345 static const char build[] = __TIME__ "/" __DATE__;
347 static const char mlist[] = "JanFebMarAprMayJunJulAugSepOctNovDec";
351 unsigned short hour, minute, second, day, year;
352 /* Note: The above quantities are used for sscanf 'hu' format,
353 * so using 'uint16_t' is contra-indicated!
357 static int ignore = 0;
366 /* check environment if build date should be ignored */
369 envstr = getenv("NTPD_IGNORE_BUILD_DATE");
370 ignore = 1 + (envstr && (!*envstr || !strcasecmp(envstr, "yes")));
376 if (6 == sscanf(build, "%hu:%hu:%hu/%3s %hu %hu",
377 &hour, &minute, &second, monstr, &day, &year)) {
378 cp = strstr(mlist, monstr);
381 jd->month = (uint8_t)((cp - mlist) / 3 + 1);
382 jd->monthday = (uint8_t)day;
383 jd->hour = (uint8_t)hour;
384 jd->minute = (uint8_t)minute;
385 jd->second = (uint8_t)second;
396 *---------------------------------------------------------------------
397 * basic calendar stuff
398 *---------------------------------------------------------------------
401 /* month table for a year starting with March,1st */
402 static const uint16_t shift_month_table[13] = {
403 0, 31, 61, 92, 122, 153, 184, 214, 245, 275, 306, 337, 366
406 /* month tables for years starting with January,1st; regular & leap */
407 static const uint16_t real_month_table[2][13] = {
408 /* -*- table for regular years -*- */
409 { 0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334, 365 },
410 /* -*- table for leap years -*- */
411 { 0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335, 366 }
415 * Some notes on the terminology:
417 * We use the proleptic Gregorian calendar, which is the Gregorian
418 * calendar extended in both directions ad infinitum. This totally
419 * disregards the fact that this calendar was invented in 1582, and
420 * was adopted at various dates over the world; sometimes even after
421 * the start of the NTP epoch.
423 * Normally date parts are given as current cycles, while time parts
424 * are given as elapsed cycles:
426 * 1970-01-01/03:04:05 means 'IN the 1970st. year, IN the first month,
427 * ON the first day, with 3hrs, 4minutes and 5 seconds elapsed.
429 * The basic calculations for this calendar implementation deal with
430 * ELAPSED date units, which is the number of full years, full months
431 * and full days before a date: 1970-01-01 would be (1969, 0, 0) in
434 * To ease the numeric computations, month and day values outside the
435 * normal range are acceptable: 2001-03-00 will be treated as the day
436 * before 2001-03-01, 2000-13-32 will give the same result as
437 * 2001-02-01 and so on.
439 * 'rd' or 'RD' is used as an abbreviation for the latin 'rata die'
440 * (day number). This is the number of days elapsed since 0000-12-31
441 * in the proleptic Gregorian calendar. The begin of the Christian Era
442 * (0001-01-01) is RD(1).
446 * ====================================================================
448 * General algorithmic stuff
450 * ====================================================================
454 *---------------------------------------------------------------------
455 * Do a periodic extension of 'value' around 'pivot' with a period of
458 * The result 'res' is a number that holds to the following properties:
460 * 1) res MOD cycle == value MOD cycle
461 * 2) pivot <= res < pivot + cycle
462 * (replace </<= with >/>= for negative cycles)
464 * where 'MOD' denotes the modulo operator for FLOOR DIVISION, which
465 * is not the same as the '%' operator in C: C requires division to be
466 * a truncated division, where remainder and dividend have the same
467 * sign if the remainder is not zero, whereas floor division requires
468 * divider and modulus to have the same sign for a non-zero modulus.
470 * This function has some useful applications:
472 * + let Y be a calendar year and V a truncated 2-digit year: then
473 * periodic_extend(Y-50, V, 100)
474 * is the closest expansion of the truncated year with respect to
475 * the full year, that is a 4-digit year with a difference of less
476 * than 50 years to the year Y. ("century unfolding")
478 * + let T be a UN*X time stamp and V be seconds-of-day: then
479 * perodic_extend(T-43200, V, 86400)
480 * is a time stamp that has the same seconds-of-day as the input
481 * value, with an absolute difference to T of <= 12hrs. ("day
484 * + Wherever you have a truncated periodic value and a non-truncated
485 * base value and you want to match them somehow...
487 * Basically, the function delivers 'pivot + (value - pivot) % cycle',
488 * but the implementation takes some pains to avoid internal signed
489 * integer overflows in the '(value - pivot) % cycle' part and adheres
490 * to the floor division convention.
492 * If 64bit scalars where available on all intended platforms, writing a
493 * version that uses 64 bit ops would be easy; writing a general
494 * division routine for 64bit ops on a platform that can only do
495 * 32/16bit divisions and is still performant is a bit more
496 * difficult. Since most usecases can be coded in a way that does only
497 * require the 32-bit version a 64bit version is NOT provided here.
498 *---------------------------------------------------------------------
501 ntpcal_periodic_extend(
508 char cpl = 0; /* modulo complement flag */
509 char neg = 0; /* sign change flag */
511 /* make the cycle positive and adjust the flags */
517 /* guard against div by zero or one */
520 * Get absolute difference as unsigned quantity and
521 * the complement flag. This is done by always
522 * subtracting the smaller value from the bigger
525 if (value >= pivot) {
526 diff = int32_to_uint32_2cpl(value)
527 - int32_to_uint32_2cpl(pivot);
529 diff = int32_to_uint32_2cpl(pivot)
530 - int32_to_uint32_2cpl(value);
533 diff %= (uint32_t)cycle;
536 diff = (uint32_t)cycle - diff;
539 pivot += uint32_2cpl_to_int32(diff);
545 /*---------------------------------------------------------------------
546 * Note to the casual reader
548 * In the next two functions you will find (or would have found...)
551 * res.Q_s -= 0x80000000;
553 * There was some ruckus about a possible programming error due to
554 * integer overflow and sign propagation.
556 * This assumption is based on a lack of understanding of the C
557 * standard. (Though this is admittedly not one of the most 'natural'
558 * aspects of the 'C' language and easily to get wrong.)
561 * http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1570.pdf
562 * "ISO/IEC 9899:201x Committee Draft — April 12, 2011"
563 * 6.4.4.1 Integer constants, clause 5
565 * why there is no sign extension/overflow problem here.
567 * But to ease the minds of the doubtful, I added back the 'u' qualifiers
568 * that somehow got lost over the last years.
573 *---------------------------------------------------------------------
574 * Convert a timestamp in NTP scale to a 64bit seconds value in the UN*X
575 * scale with proper epoch unfolding around a given pivot or the current
576 * system time. This function happily accepts negative pivot values as
577 * timestamps befor 1970-01-01, so be aware of possible trouble on
578 * platforms with 32bit 'time_t'!
580 * This is also a periodic extension, but since the cycle is 2^32 and
581 * the shift is 2^31, we can do some *very* fast math without explicit
583 *---------------------------------------------------------------------
593 # if defined(HAVE_INT64)
595 res.q_s = (pivot != NULL)
598 res.Q_s -= 0x80000000u; /* unshift of half range */
599 ntp -= (uint32_t)JAN_1970; /* warp into UN*X domain */
600 ntp -= res.D_s.lo; /* cycle difference */
601 res.Q_s += (uint64_t)ntp; /* get expanded time */
603 # else /* no 64bit scalars */
607 tmp = (pivot != NULL)
610 res = time_to_vint64(&tmp);
611 M_SUB(res.D_s.hi, res.D_s.lo, 0, 0x80000000u);
612 ntp -= (uint32_t)JAN_1970; /* warp into UN*X domain */
613 ntp -= res.D_s.lo; /* cycle difference */
614 M_ADD(res.D_s.hi, res.D_s.lo, 0, ntp);
616 # endif /* no 64bit scalars */
622 *---------------------------------------------------------------------
623 * Convert a timestamp in NTP scale to a 64bit seconds value in the NTP
624 * scale with proper epoch unfolding around a given pivot or the current
627 * Note: The pivot must be given in the UN*X time domain!
629 * This is also a periodic extension, but since the cycle is 2^32 and
630 * the shift is 2^31, we can do some *very* fast math without explicit
632 *---------------------------------------------------------------------
642 # if defined(HAVE_INT64)
647 res.Q_s -= 0x80000000u; /* unshift of half range */
648 res.Q_s += (uint32_t)JAN_1970; /* warp into NTP domain */
649 ntp -= res.D_s.lo; /* cycle difference */
650 res.Q_s += (uint64_t)ntp; /* get expanded time */
652 # else /* no 64bit scalars */
659 res = time_to_vint64(&tmp);
660 M_SUB(res.D_s.hi, res.D_s.lo, 0, 0x80000000u);
661 M_ADD(res.D_s.hi, res.D_s.lo, 0, (uint32_t)JAN_1970);/*into NTP */
662 ntp -= res.D_s.lo; /* cycle difference */
663 M_ADD(res.D_s.hi, res.D_s.lo, 0, ntp);
665 # endif /* no 64bit scalars */
672 * ====================================================================
674 * Splitting values to composite entities
676 * ====================================================================
680 *---------------------------------------------------------------------
681 * Split a 64bit seconds value into elapsed days in 'res.hi' and
682 * elapsed seconds since midnight in 'res.lo' using explicit floor
683 * division. This function happily accepts negative time values as
684 * timestamps before the respective epoch start.
685 *---------------------------------------------------------------------
695 # if defined(HAVE_INT64)
697 /* Manual floor division by SECSPERDAY. This uses the one's
698 * complement trick, too, but without an extra flag value: The
699 * flag would be 64bit, and that's a bit of overkill on a 32bit
700 * target that has to use a register pair for a 64bit number.
703 Q = ~(uint32_t)(~ts->Q_s / SECSPERDAY);
705 Q = (uint32_t)(ts->Q_s / SECSPERDAY);
709 uint32_t ah, al, sflag, A;
711 /* get operand into ah/al (either ts or ts' one's complement,
712 * for later floor division)
714 sflag = int32_sflag(ts->d_s.hi);
715 ah = sflag ^ ts->D_s.hi;
716 al = sflag ^ ts->D_s.lo;
718 /* Since 86400 == 128*675 we can drop the least 7 bits and
719 * divide by 675 instead of 86400. Then the maximum remainder
720 * after each devision step is 674, and we need 10 bits for
721 * that. So in the next step we can shift in 22 bits from the
724 * Therefore we load the accu with the top 13 bits (51..63) in
725 * the first shot. We don't have to remember the quotient -- it
726 * would be shifted out anyway.
732 /* Now assemble the remainder with bits 29..50 from the
733 * numerator and divide. This creates the upper ten bits of the
734 * quotient. (Well, the top 22 bits of a 44bit result. But that
735 * will be truncated to 32 bits anyway.)
737 A = (A << 19) | (ah & 0x0007FFFFu);
738 A = (A << 3) | (al >> 29);
742 /* Now assemble the remainder with bits 7..28 from the numerator
743 * and do a final division step.
745 A = (A << 22) | ((al >> 7) & 0x003FFFFFu);
746 Q = (Q << 22) | (A / 675u);
748 /* The last 7 bits get simply dropped, as they have no affect on
749 * the quotient when dividing by 86400.
752 /* apply sign correction and calculate the true floor
759 res.hi = uint32_2cpl_to_int32(Q);
760 res.lo = ts->D_s.lo - Q * SECSPERDAY;
766 *---------------------------------------------------------------------
767 * Split a 32bit seconds value into h/m/s and excessive days. This
768 * function happily accepts negative time values as timestamps before
770 *---------------------------------------------------------------------
778 /* Do 3 chained floor divisions by positive constants, using the
779 * one's complement trick and factoring out the intermediate XOR
780 * ops to reduce the number of operations.
782 uint32_t us, um, uh, ud, sflag;
784 sflag = int32_sflag(ts);
785 us = int32_to_uint32_2cpl(ts);
787 um = (sflag ^ us) / SECSPERMIN;
795 split[0] = (int32_t)(uh - ud * HRSPERDAY );
796 split[1] = (int32_t)(um - uh * MINSPERHR );
797 split[2] = (int32_t)(us - um * SECSPERMIN);
799 return uint32_2cpl_to_int32(ud);
803 *---------------------------------------------------------------------
804 * Given the number of elapsed days in the calendar era, split this
805 * number into the number of elapsed years in 'res.hi' and the number
806 * of elapsed days of that year in 'res.lo'.
808 * if 'isleapyear' is not NULL, it will receive an integer that is 0 for
809 * regular years and a non-zero value for leap years.
810 *---------------------------------------------------------------------
813 ntpcal_split_eradays(
818 /* Use the fast cyclesplit algorithm here, to calculate the
819 * centuries and years in a century with one division each. This
820 * reduces the number of division operations to two, but is
821 * susceptible to internal range overflow. We make sure the
822 * input operands are in the safe range; this still gives us
823 * approx +/-2.9 million years.
826 int32_t n100, n001; /* calendar year cycles */
827 uint32_t uday, Q, sflag;
829 /* split off centuries first */
830 sflag = int32_sflag(days);
831 uday = uint32_saturate(int32_to_uint32_2cpl(days), sflag);
832 uday = (4u * uday) | 3u;
833 Q = sflag ^ ((sflag ^ uday) / GREGORIAN_CYCLE_DAYS);
834 uday = uday - Q * GREGORIAN_CYCLE_DAYS;
835 n100 = uint32_2cpl_to_int32(Q);
837 /* Split off years in century -- days >= 0 here, and we're far
838 * away from integer overflow trouble now. */
840 n001 = uday / GREGORIAN_NORMAL_LEAP_CYCLE_DAYS;
841 uday = uday % GREGORIAN_NORMAL_LEAP_CYCLE_DAYS;
843 /* Assemble the year and day in year */
844 res.hi = n100 * 100 + n001;
847 /* Eventually set the leap year flag. Note: 0 <= n001 <= 99 and
848 * Q is still the two's complement representation of the
849 * centuries: The modulo 4 ops can be done with masking here.
850 * We also shift the year and the century by one, so the tests
851 * can be done against zero instead of 3.
854 *isleapyear = !((n001+1) & 3)
855 && ((n001 != 99) || !((Q+1) & 3));
861 *---------------------------------------------------------------------
862 * Given a number of elapsed days in a year and a leap year indicator,
863 * split the number of elapsed days into the number of elapsed months in
864 * 'res.hi' and the number of elapsed days of that month in 'res.lo'.
866 * This function will fail and return {-1,-1} if the number of elapsed
867 * days is not in the valid range!
868 *---------------------------------------------------------------------
871 ntpcal_split_yeardays(
877 const uint16_t *lt; /* month length table */
879 /* check leap year flag and select proper table */
880 lt = real_month_table[(isleapyear != 0)];
881 if (0 <= eyd && eyd < lt[12]) {
882 /* get zero-based month by approximation & correction step */
883 res.hi = eyd >> 5; /* approx month; might be 1 too low */
884 if (lt[res.hi + 1] <= eyd) /* fixup approximative month value */
886 res.lo = eyd - lt[res.hi];
888 res.lo = res.hi = -1;
895 *---------------------------------------------------------------------
896 * Convert a RD into the date part of a 'struct calendar'.
897 *---------------------------------------------------------------------
909 /* Get day-of-week first. Since rd is signed, the remainder can
910 * be in the range [-6..+6], but the assignment to an unsigned
911 * variable maps the negative values to positive values >=7.
912 * This makes the sign correction look strange, but adding 7
913 * causes the needed wrap-around into the desired value range of
914 * zero to six, both inclusive.
916 jd->weekday = rd % DAYSPERWEEK;
917 if (jd->weekday >= DAYSPERWEEK) /* weekday is unsigned! */
918 jd->weekday += DAYSPERWEEK;
920 split = ntpcal_split_eradays(rd - 1, &leapy);
921 /* Get year and day-of-year, with overflow check. If any of the
922 * upper 16 bits is set after shifting to unity-based years, we
923 * will have an overflow when converting to an unsigned 16bit
924 * year. Shifting to the right is OK here, since it does not
925 * matter if the shift is logic or arithmetic.
928 ymask = 0u - ((split.hi >> 16) == 0);
929 jd->year = (uint16_t)(split.hi & ymask);
930 jd->yearday = (uint16_t)split.lo + 1;
932 /* convert to month and mday */
933 split = ntpcal_split_yeardays(split.lo, leapy);
934 jd->month = (uint8_t)split.hi + 1;
935 jd->monthday = (uint8_t)split.lo + 1;
937 return ymask ? leapy : -1;
941 *---------------------------------------------------------------------
942 * Convert a RD into the date part of a 'struct tm'.
943 *---------------------------------------------------------------------
954 /* get day-of-week first */
955 utm->tm_wday = rd % DAYSPERWEEK;
956 if (utm->tm_wday < 0)
957 utm->tm_wday += DAYSPERWEEK;
959 /* get year and day-of-year */
960 split = ntpcal_split_eradays(rd - 1, &leapy);
961 utm->tm_year = split.hi - 1899;
962 utm->tm_yday = split.lo; /* 0-based */
964 /* convert to month and mday */
965 split = ntpcal_split_yeardays(split.lo, leapy);
966 utm->tm_mon = split.hi; /* 0-based */
967 utm->tm_mday = split.lo + 1; /* 1-based */
973 *---------------------------------------------------------------------
974 * Take a value of seconds since midnight and split it into hhmmss in a
976 *---------------------------------------------------------------------
979 ntpcal_daysec_to_date(
987 days = priv_timesplit(ts, sec);
988 jd->hour = (uint8_t)ts[0];
989 jd->minute = (uint8_t)ts[1];
990 jd->second = (uint8_t)ts[2];
996 *---------------------------------------------------------------------
997 * Take a value of seconds since midnight and split it into hhmmss in a
999 *---------------------------------------------------------------------
1002 ntpcal_daysec_to_tm(
1010 days = priv_timesplit(ts, sec);
1011 utm->tm_hour = ts[0];
1012 utm->tm_min = ts[1];
1013 utm->tm_sec = ts[2];
1019 *---------------------------------------------------------------------
1020 * take a split representation for day/second-of-day and day offset
1021 * and convert it to a 'struct calendar'. The seconds will be normalised
1022 * into the range of a day, and the day will be adjusted accordingly.
1024 * returns >0 if the result is in a leap year, 0 if in a regular
1025 * year and <0 if the result did not fit into the calendar struct.
1026 *---------------------------------------------------------------------
1029 ntpcal_daysplit_to_date(
1030 struct calendar *jd,
1031 const ntpcal_split *ds,
1035 dof += ntpcal_daysec_to_date(jd, ds->lo);
1036 return ntpcal_rd_to_date(jd, ds->hi + dof);
1040 *---------------------------------------------------------------------
1041 * take a split representation for day/second-of-day and day offset
1042 * and convert it to a 'struct tm'. The seconds will be normalised
1043 * into the range of a day, and the day will be adjusted accordingly.
1045 * returns 1 if the result is in a leap year and zero if in a regular
1047 *---------------------------------------------------------------------
1050 ntpcal_daysplit_to_tm(
1052 const ntpcal_split *ds ,
1056 dof += ntpcal_daysec_to_tm(utm, ds->lo);
1058 return ntpcal_rd_to_tm(utm, ds->hi + dof);
1062 *---------------------------------------------------------------------
1063 * Take a UN*X time and convert to a calendar structure.
1064 *---------------------------------------------------------------------
1067 ntpcal_time_to_date(
1068 struct calendar *jd,
1074 ds = ntpcal_daysplit(ts);
1075 ds.hi += ntpcal_daysec_to_date(jd, ds.lo);
1076 ds.hi += DAY_UNIX_STARTS;
1078 return ntpcal_rd_to_date(jd, ds.hi);
1083 * ====================================================================
1085 * merging composite entities
1087 * ====================================================================
1091 *---------------------------------------------------------------------
1092 * Merge a number of days and a number of seconds into seconds,
1093 * expressed in 64 bits to avoid overflow.
1094 *---------------------------------------------------------------------
1104 # if defined(HAVE_INT64)
1107 res.q_s *= SECSPERDAY;
1116 * res = days *86400 + secs, using manual 16/32 bit
1117 * multiplications and shifts.
1123 /* assemble days * 675 */
1124 res.D_s.lo = (days & 0xFFFF) * 675u;
1126 p1 = (days >> 16) * 675u;
1129 M_ADD(res.D_s.hi, res.D_s.lo, p2, p1);
1131 /* mul by 128, using shift */
1132 res.D_s.hi = (res.D_s.hi << 7) | (res.D_s.lo >> 25);
1133 res.D_s.lo = (res.D_s.lo << 7);
1137 M_NEG(res.D_s.hi, res.D_s.lo);
1139 /* properly add seconds */
1142 p1 = (uint32_t)-secs;
1145 p1 = (uint32_t)secs;
1147 M_ADD(res.D_s.hi, res.D_s.lo, p2, p1);
1155 *---------------------------------------------------------------------
1156 * get leap years since epoch in elapsed years
1157 *---------------------------------------------------------------------
1160 ntpcal_leapyears_in_years(
1164 /* We use the in-out-in algorithm here, using the one's
1165 * complement division trick for negative numbers. The chained
1166 * division sequence by 4/25/4 gives the compiler the chance to
1167 * get away with only one true division and doing shifts otherwise.
1170 uint32_t sflag, sum, uyear;
1172 sflag = int32_sflag(years);
1173 uyear = int32_to_uint32_2cpl(years);
1176 sum = (uyear /= 4u); /* 4yr rule --> IN */
1177 sum -= (uyear /= 25u); /* 100yr rule --> OUT */
1178 sum += (uyear /= 4u); /* 400yr rule --> IN */
1180 /* Thanks to the alternation of IN/OUT/IN we can do the sum
1181 * directly and have a single one's complement operation
1182 * here. (Only if the years are negative, of course.) Otherwise
1183 * the one's complement would have to be done when
1184 * adding/subtracting the terms.
1186 return uint32_2cpl_to_int32(sflag ^ sum);
1190 *---------------------------------------------------------------------
1191 * Convert elapsed years in Era into elapsed days in Era.
1192 *---------------------------------------------------------------------
1195 ntpcal_days_in_years(
1199 return years * DAYSPERYEAR + ntpcal_leapyears_in_years(years);
1203 *---------------------------------------------------------------------
1204 * Convert a number of elapsed month in a year into elapsed days in year.
1206 * The month will be normalized, and 'res.hi' will contain the
1207 * excessive years that must be considered when converting the years,
1208 * while 'res.lo' will contain the number of elapsed days since start
1211 * This code uses the shifted-month-approach to convert month to days,
1212 * because then there is no need to have explicit leap year
1213 * information. The slight disadvantage is that for most month values
1214 * the result is a negative value, and the year excess is one; the
1215 * conversion is then simply based on the start of the following year.
1216 *---------------------------------------------------------------------
1219 ntpcal_days_in_months(
1225 /* Add ten months and correct if needed. (It likely is...) */
1227 res.hi = (res.lo >= 12);
1231 /* if still out of range, normalise by floor division ... */
1232 if (res.lo < 0 || res.lo >= 12) {
1233 uint32_t mu, Q, sflag;
1234 sflag = int32_sflag(res.lo);
1235 mu = int32_to_uint32_2cpl(res.lo);
1236 Q = sflag ^ ((sflag ^ mu) / 12u);
1237 res.hi += uint32_2cpl_to_int32(Q);
1238 res.lo = mu - Q * 12u;
1241 /* get cummulated days in year with unshift */
1242 res.lo = shift_month_table[res.lo] - 306;
1248 *---------------------------------------------------------------------
1249 * Convert ELAPSED years/months/days of gregorian calendar to elapsed
1250 * days in Gregorian epoch.
1252 * If you want to convert years and days-of-year, just give a month of
1254 *---------------------------------------------------------------------
1257 ntpcal_edate_to_eradays(
1267 tmp = ntpcal_days_in_months(mons);
1268 res = ntpcal_days_in_years(years + tmp.hi) + tmp.lo;
1270 res = ntpcal_days_in_years(years);
1277 *---------------------------------------------------------------------
1278 * Convert ELAPSED years/months/days of gregorian calendar to elapsed
1281 * Note: This will give the true difference to the start of the given
1282 * year, even if months & days are off-scale.
1283 *---------------------------------------------------------------------
1286 ntpcal_edate_to_yeardays(
1294 if (0 <= mons && mons < 12) {
1296 mdays += real_month_table[is_leapyear(years)][mons];
1298 tmp = ntpcal_days_in_months(mons);
1300 + ntpcal_days_in_years(years + tmp.hi)
1301 - ntpcal_days_in_years(years);
1308 *---------------------------------------------------------------------
1309 * Convert elapsed days and the hour/minute/second information into
1312 * If 'isvalid' is not NULL, do a range check on the time specification
1313 * and tell if the time input is in the normal range, permitting for a
1314 * single leapsecond.
1315 *---------------------------------------------------------------------
1318 ntpcal_etime_to_seconds(
1326 res = (hours * MINSPERHR + minutes) * SECSPERMIN + seconds;
1332 *---------------------------------------------------------------------
1333 * Convert the date part of a 'struct tm' (that is, year, month,
1334 * day-of-month) into the RD of that day.
1335 *---------------------------------------------------------------------
1339 const struct tm *utm
1342 return ntpcal_edate_to_eradays(utm->tm_year + 1899,
1344 utm->tm_mday - 1) + 1;
1348 *---------------------------------------------------------------------
1349 * Convert the date part of a 'struct calendar' (that is, year, month,
1350 * day-of-month) into the RD of that day.
1351 *---------------------------------------------------------------------
1355 const struct calendar *jd
1358 return ntpcal_edate_to_eradays((int32_t)jd->year - 1,
1359 (int32_t)jd->month - 1,
1360 (int32_t)jd->monthday - 1) + 1;
1364 *---------------------------------------------------------------------
1365 * convert a year number to rata die of year start
1366 *---------------------------------------------------------------------
1369 ntpcal_year_to_ystart(
1373 return ntpcal_days_in_years(year - 1) + 1;
1377 *---------------------------------------------------------------------
1378 * For a given RD, get the RD of the associated year start,
1379 * that is, the RD of the last January,1st on or before that day.
1380 *---------------------------------------------------------------------
1383 ntpcal_rd_to_ystart(
1388 * Rather simple exercise: split the day number into elapsed
1389 * years and elapsed days, then remove the elapsed days from the
1390 * input value. Nice'n sweet...
1392 return rd - ntpcal_split_eradays(rd - 1, NULL).lo;
1396 *---------------------------------------------------------------------
1397 * For a given RD, get the RD of the associated month start.
1398 *---------------------------------------------------------------------
1401 ntpcal_rd_to_mstart(
1408 split = ntpcal_split_eradays(rd - 1, &leaps);
1409 split = ntpcal_split_yeardays(split.lo, leaps);
1411 return rd - split.lo;
1415 *---------------------------------------------------------------------
1416 * take a 'struct calendar' and get the seconds-of-day from it.
1417 *---------------------------------------------------------------------
1420 ntpcal_date_to_daysec(
1421 const struct calendar *jd
1424 return ntpcal_etime_to_seconds(jd->hour, jd->minute,
1429 *---------------------------------------------------------------------
1430 * take a 'struct tm' and get the seconds-of-day from it.
1431 *---------------------------------------------------------------------
1434 ntpcal_tm_to_daysec(
1435 const struct tm *utm
1438 return ntpcal_etime_to_seconds(utm->tm_hour, utm->tm_min,
1443 *---------------------------------------------------------------------
1444 * take a 'struct calendar' and convert it to a 'time_t'
1445 *---------------------------------------------------------------------
1448 ntpcal_date_to_time(
1449 const struct calendar *jd
1455 days = ntpcal_date_to_rd(jd) - DAY_UNIX_STARTS;
1456 secs = ntpcal_date_to_daysec(jd);
1457 join = ntpcal_dayjoin(days, secs);
1459 return vint64_to_time(&join);
1464 * ====================================================================
1466 * extended and unchecked variants of caljulian/caltontp
1468 * ====================================================================
1471 ntpcal_ntp64_to_date(
1472 struct calendar *jd,
1478 ds = ntpcal_daysplit(ntp);
1479 ds.hi += ntpcal_daysec_to_date(jd, ds.lo);
1481 return ntpcal_rd_to_date(jd, ds.hi + DAY_NTP_STARTS);
1486 struct calendar *jd,
1494 * Unfold ntp time around current time into NTP domain. Split
1495 * into days and seconds, shift days into CE domain and
1496 * process the parts.
1498 ntp64 = ntpcal_ntp_to_ntp(ntp, piv);
1499 return ntpcal_ntp64_to_date(jd, &ntp64);
1504 ntpcal_date_to_ntp64(
1505 const struct calendar *jd
1509 * Convert date to NTP. Ignore yearday, use d/m/y only.
1511 return ntpcal_dayjoin(ntpcal_date_to_rd(jd) - DAY_NTP_STARTS,
1512 ntpcal_date_to_daysec(jd));
1518 const struct calendar *jd
1522 * Get lower half of 64-bit NTP timestamp from date/time.
1524 return ntpcal_date_to_ntp64(jd).d_s.lo;
1530 * ====================================================================
1532 * day-of-week calculations
1534 * ====================================================================
1537 * Given a RataDie and a day-of-week, calculate a RDN that is reater-than,
1538 * greater-or equal, closest, less-or-equal or less-than the given RDN
1539 * and denotes the given day-of-week
1547 return ntpcal_periodic_extend(rdn+1, dow, 7);
1556 return ntpcal_periodic_extend(rdn, dow, 7);
1560 ntpcal_weekday_close(
1565 return ntpcal_periodic_extend(rdn-3, dow, 7);
1574 return ntpcal_periodic_extend(rdn, dow, -7);
1583 return ntpcal_periodic_extend(rdn-1, dow, -7);
1587 * ====================================================================
1589 * ISO week-calendar conversions
1591 * The ISO8601 calendar defines a calendar of years, weeks and weekdays.
1592 * It is related to the Gregorian calendar, and a ISO year starts at the
1593 * Monday closest to Jan,1st of the corresponding Gregorian year. A ISO
1594 * calendar year has always 52 or 53 weeks, and like the Grogrian
1595 * calendar the ISO8601 calendar repeats itself every 400 years, or
1596 * 146097 days, or 20871 weeks.
1598 * While it is possible to write ISO calendar functions based on the
1599 * Gregorian calendar functions, the following implementation takes a
1600 * different approach, based directly on years and weeks.
1602 * Analysis of the tabulated data shows that it is not possible to
1603 * interpolate from years to weeks over a full 400 year range; cyclic
1604 * shifts over 400 years do not provide a solution here. But it *is*
1605 * possible to interpolate over every single century of the 400-year
1606 * cycle. (The centennial leap year rule seems to be the culprit here.)
1608 * It can be shown that a conversion from years to weeks can be done
1609 * using a linear transformation of the form
1611 * w = floor( y * a + b )
1613 * where the slope a must hold to
1615 * 52.1780821918 <= a < 52.1791044776
1617 * and b must be chosen according to the selected slope and the number
1618 * of the century in a 400-year period.
1620 * The inverse calculation can also be done in this way. Careful scaling
1621 * provides an unlimited set of integer coefficients a,k,b that enable
1622 * us to write the calulation in the form
1624 * w = (y * a + b ) / k
1625 * y = (w * a' + b') / k'
1627 * In this implementation the values of k and k' are chosen to be
1628 * smallest possible powers of two, so the division can be implemented
1629 * as shifts if the optimiser chooses to do so.
1631 * ====================================================================
1635 * Given a number of elapsed (ISO-)years since the begin of the
1636 * christian era, return the number of elapsed weeks corresponding to
1637 * the number of years.
1640 isocal_weeks_in_years(
1645 * use: w = (y * 53431 + b[c]) / 1024 as interpolation
1647 static const uint16_t bctab[4] = { 157, 449, 597, 889 };
1650 uint32_t cc, ci, yu, sflag;
1652 sflag = int32_sflag(years);
1653 yu = int32_to_uint32_2cpl(years);
1655 /* split off centuries, using floor division */
1656 cc = sflag ^ ((sflag ^ yu) / 100u);
1659 /* calculate century cycles shift and cycle index:
1660 * Assuming a century is 5217 weeks, we have to add a cycle
1661 * shift that is 3 for every 4 centuries, because 3 of the four
1662 * centuries have 5218 weeks. So '(cc*3 + 1) / 4' is the actual
1663 * correction, and the second century is the defective one.
1665 * Needs floor division by 4, which is done with masking and
1669 cs = uint32_2cpl_to_int32(sflag ^ ((sflag ^ ci) / 4u));
1672 /* Get weeks in century. Can use plain division here as all ops
1673 * are >= 0, and let the compiler sort out the possible
1676 cw = (yu * 53431u + bctab[ci]) / 1024u;
1678 return uint32_2cpl_to_int32(cc) * 5217 + cs + cw;
1682 * Given a number of elapsed weeks since the begin of the christian
1683 * era, split this number into the number of elapsed years in res.hi
1684 * and the excessive number of weeks in res.lo. (That is, res.lo is
1685 * the number of elapsed weeks in the remaining partial year.)
1688 isocal_split_eraweeks(
1693 * use: y = (w * 157 + b[c]) / 8192 as interpolation
1696 static const uint16_t bctab[4] = { 85, 130, 17, 62 };
1700 uint32_t sw, cy, Q, sflag;
1702 /* Use two fast cycle-split divisions here. This is again
1703 * susceptible to internal overflow, so we check the range. This
1704 * still permits more than +/-20 million years, so this is
1705 * likely a pure academical problem.
1707 * We want to execute '(weeks * 4 + 2) /% 20871' under floor
1708 * division rules in the first step.
1710 sflag = int32_sflag(weeks);
1711 sw = uint32_saturate(int32_to_uint32_2cpl(weeks), sflag);
1713 Q = sflag ^ ((sflag ^ sw) / GREGORIAN_CYCLE_WEEKS);
1714 sw -= Q * GREGORIAN_CYCLE_WEEKS;
1716 cc = uint32_2cpl_to_int32(Q);
1718 /* Split off years; sw >= 0 here! The scaled weeks in the years
1719 * are scaled up by 157 afterwards.
1721 sw = (sw / 4u) * 157u + bctab[ci];
1722 cy = sw / 8192u; /* ws >> 13 , let the compiler sort it out */
1723 sw = sw % 8192u; /* ws & 8191, let the compiler sort it out */
1725 /* assemble elapsed years and downscale the elapsed weeks in
1728 res.hi = 100*cc + cy;
1735 * Given a second in the NTP time scale and a pivot, expand the NTP
1736 * time stamp around the pivot and convert into an ISO calendar time
1740 isocal_ntp64_to_date(
1747 uint32_t uw, ud, sflag;
1750 * Split NTP time into days and seconds, shift days into CE
1751 * domain and process the parts.
1753 ds = ntpcal_daysplit(ntp);
1755 /* split time part */
1756 ds.hi += priv_timesplit(ts, ds.lo);
1757 id->hour = (uint8_t)ts[0];
1758 id->minute = (uint8_t)ts[1];
1759 id->second = (uint8_t)ts[2];
1761 /* split days into days and weeks, using floor division in unsigned */
1762 ds.hi += DAY_NTP_STARTS - 1; /* shift from NTP to RDN */
1763 sflag = int32_sflag(ds.hi);
1764 ud = int32_to_uint32_2cpl(ds.hi);
1765 uw = sflag ^ ((sflag ^ ud) / DAYSPERWEEK);
1766 ud -= uw * DAYSPERWEEK;
1767 ds.hi = uint32_2cpl_to_int32(uw);
1770 id->weekday = (uint8_t)ds.lo + 1; /* weekday result */
1772 /* get year and week in year */
1773 ds = isocal_split_eraweeks(ds.hi); /* elapsed years&week*/
1774 id->year = (uint16_t)ds.hi + 1; /* shift to current */
1775 id->week = (uint8_t )ds.lo + 1;
1777 return (ds.hi >= 0 && ds.hi < 0x0000FFFF);
1790 * Unfold ntp time around current time into NTP domain, then
1791 * convert the full time stamp.
1793 ntp64 = ntpcal_ntp_to_ntp(ntp, piv);
1794 return isocal_ntp64_to_date(id, &ntp64);
1798 * Convert a ISO date spec into a second in the NTP time scale,
1799 * properly truncated to 32 bit.
1802 isocal_date_to_ntp64(
1803 const struct isodate *id
1806 int32_t weeks, days, secs;
1808 weeks = isocal_weeks_in_years((int32_t)id->year - 1)
1809 + (int32_t)id->week - 1;
1810 days = weeks * 7 + (int32_t)id->weekday;
1811 /* days is RDN of ISO date now */
1812 secs = ntpcal_etime_to_seconds(id->hour, id->minute, id->second);
1814 return ntpcal_dayjoin(days - DAY_NTP_STARTS, secs);
1819 const struct isodate *id
1823 * Get lower half of 64-bit NTP timestamp from date/time.
1825 return isocal_date_to_ntp64(id).d_s.lo;