Home | History | Annotate | Line # | Download | only in libntp
ntp_calendar.c revision 1.1.1.7
      1 /*
      2  * ntp_calendar.c - calendar and helper functions
      3  *
      4  * Written by Juergen Perlinger (perlinger (at) ntp.org) for the NTP project.
      5  * The contents of 'html/copyright.html' apply.
      6  *
      7  * --------------------------------------------------------------------
      8  * Some notes on the implementation:
      9  *
     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.
     18  *
     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...
     22  *
     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
     28  * purpose.
     29  *
     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.
     36  *
     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.
     41  *
     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.
     48  *
     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.
     52  *
     53  */
     54 
     55 #include <config.h>
     56 #include <sys/types.h>
     57 
     58 #include "ntp_types.h"
     59 #include "ntp_calendar.h"
     60 #include "ntp_stdlib.h"
     61 #include "ntp_fp.h"
     62 #include "ntp_unixtime.h"
     63 
     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.
     68  */
     69 #ifndef TARGET_HAS_2CPL
     70 # if defined(__GNUC__)
     71 #  if defined(__i386__) || defined(__x86_64__) || defined(__arm__)
     72 #   define TARGET_HAS_2CPL 1
     73 #  else
     74 #   define TARGET_HAS_2CPL 0
     75 #  endif
     76 # elif defined(_MSC_VER)
     77 #  if defined(_M_IX86) || defined(_M_X64) || defined(_M_ARM)
     78 #   define TARGET_HAS_2CPL 1
     79 #  else
     80 #   define TARGET_HAS_2CPL 0
     81 #  endif
     82 # else
     83 #  define TARGET_HAS_2CPL 0
     84 # endif
     85 #endif
     86 
     87 #ifndef TARGET_HAS_SAR
     88 # define TARGET_HAS_SAR 0
     89 #endif
     90 
     91 /*
     92  *---------------------------------------------------------------------
     93  * replacing the 'time()' function
     94  * --------------------------------------------------------------------
     95  */
     96 
     97 static systime_func_ptr systime_func = &time;
     98 static inline time_t now(void);
     99 
    100 
    101 systime_func_ptr
    102 ntpcal_set_timefunc(
    103 	systime_func_ptr nfunc
    104 	)
    105 {
    106 	systime_func_ptr res;
    107 
    108 	res = systime_func;
    109 	if (NULL == nfunc)
    110 		nfunc = &time;
    111 	systime_func = nfunc;
    112 
    113 	return res;
    114 }
    115 
    116 
    117 static inline time_t
    118 now(void)
    119 {
    120 	return (*systime_func)(NULL);
    121 }
    122 
    123 /*
    124  *---------------------------------------------------------------------
    125  * Get sign extension mask and unsigned 2cpl rep for a signed integer
    126  *---------------------------------------------------------------------
    127  */
    128 
    129 static inline uint32_t
    130 int32_sflag(
    131 	const int32_t v)
    132 {
    133 #   if TARGET_HAS_2CPL && TARGET_HAS_SAR && SIZEOF_INT >= 4
    134 
    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.
    140 	 */
    141 	return (uint32_t)(v >> 31);
    142 
    143 #   else
    144 
    145 	/* This should be a rather generic approach for getting a sign
    146 	 * extension mask...
    147 	 */
    148 	return UINT32_C(0) - (uint32_t)(v < 0);
    149 
    150 #   endif
    151 }
    152 
    153 static inline uint32_t
    154 int32_to_uint32_2cpl(
    155 	const int32_t v)
    156 {
    157 	uint32_t vu;
    158 
    159 #   if TARGET_HAS_2CPL
    160 
    161 	/* Just copy through the 32 bits from the signed value if we're
    162 	 * on a two's complement target.
    163 	 */
    164 	vu = (uint32_t)v;
    165 
    166 #   else
    167 
    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.
    174 	 */
    175 	if (v < 0)
    176 		vu = ~(uint32_t)(-(v + 1));
    177 	else
    178 		vu = (uint32_t)v;
    179 
    180 #   endif
    181 
    182 	return vu;
    183 }
    184 
    185 static inline int32_t
    186 uint32_2cpl_to_int32(
    187 	const uint32_t vu)
    188 {
    189 	int32_t v;
    190 
    191 #   if TARGET_HAS_2CPL
    192 
    193 	/* Just copy through the 32 bits from the unsigned value if
    194 	 * we're on a two's complement target.
    195 	 */
    196 	v = (int32_t)vu;
    197 
    198 #   else
    199 
    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.
    204 	 */
    205 	if (vu > INT32_MAX)
    206 		v = -(int32_t)(~vu) - 1;
    207 	else
    208 		v = (int32_t)vu;
    209 
    210 #   endif
    211 
    212 	return v;
    213 }
    214 
    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.
    220  */
    221 static inline uint32_t
    222 uint32_saturate(
    223 	uint32_t vu,
    224 	uint32_t mu)
    225 {
    226 	static const uint32_t limit = UINT32_MAX/4u;
    227 	if ((mu ^ vu) > limit) {
    228 		vu    = mu ^ limit;
    229 		errno = EDOM;
    230 	}
    231 	return vu;
    232 }
    233 
    234 /*
    235  *---------------------------------------------------------------------
    236  * Convert between 'time_t' and 'vint64'
    237  *---------------------------------------------------------------------
    238  */
    239 vint64
    240 time_to_vint64(
    241 	const time_t * ptt
    242 	)
    243 {
    244 	vint64 res;
    245 	time_t tt;
    246 
    247 	tt = *ptt;
    248 
    249 #   if SIZEOF_TIME_T <= 4
    250 
    251 	res.D_s.hi = 0;
    252 	if (tt < 0) {
    253 		res.D_s.lo = (uint32_t)-tt;
    254 		M_NEG(res.D_s.hi, res.D_s.lo);
    255 	} else {
    256 		res.D_s.lo = (uint32_t)tt;
    257 	}
    258 
    259 #   elif defined(HAVE_INT64)
    260 
    261 	res.q_s = tt;
    262 
    263 #   else
    264 	/*
    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!
    268 	 */
    269 	if (tt < 0) {
    270 		tt = -tt;
    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);
    274 	} else {
    275 		res.D_s.lo = (uint32_t)tt;
    276 		res.D_s.hi = (uint32_t)(tt >> 32);
    277 	}
    278 
    279 #   endif
    280 
    281 	return res;
    282 }
    283 
    284 
    285 time_t
    286 vint64_to_time(
    287 	const vint64 *tv
    288 	)
    289 {
    290 	time_t res;
    291 
    292 #   if SIZEOF_TIME_T <= 4
    293 
    294 	res = (time_t)tv->D_s.lo;
    295 
    296 #   elif defined(HAVE_INT64)
    297 
    298 	res = (time_t)tv->q_s;
    299 
    300 #   else
    301 
    302 	res = ((time_t)tv->d_s.hi << 32) | tv->D_s.lo;
    303 
    304 #   endif
    305 
    306 	return res;
    307 }
    308 
    309 /*
    310  *---------------------------------------------------------------------
    311  * Get the build date & time
    312  *---------------------------------------------------------------------
    313  */
    314 int
    315 ntpcal_get_build_date(
    316 	struct calendar * jd
    317 	)
    318 {
    319 	/* The C standard tells us the format of '__DATE__':
    320 	 *
    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.
    328 	 *
    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.
    334 	 *
    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
    339 	 * problem.
    340 	 *
    341 	 */
    342 #   ifdef MKREPRO_DATE
    343 	static const char build[] = MKREPRO_TIME "/" MKREPRO_DATE;
    344 #   else
    345 	static const char build[] = __TIME__ "/" __DATE__;
    346 #   endif
    347 	static const char mlist[] = "JanFebMarAprMayJunJulAugSepOctNovDec";
    348 
    349 	char		  monstr[4];
    350 	const char *	  cp;
    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!
    354 	 */
    355 
    356 #   ifdef DEBUG
    357 	static int        ignore  = 0;
    358 #   endif
    359 
    360 	ZERO(*jd);
    361 	jd->year     = 1970;
    362 	jd->month    = 1;
    363 	jd->monthday = 1;
    364 
    365 #   ifdef DEBUG
    366 	/* check environment if build date should be ignored */
    367 	if (0 == ignore) {
    368 	    const char * envstr;
    369 	    envstr = getenv("NTPD_IGNORE_BUILD_DATE");
    370 	    ignore = 1 + (envstr && (!*envstr || !strcasecmp(envstr, "yes")));
    371 	}
    372 	if (ignore > 1)
    373 	    return FALSE;
    374 #   endif
    375 
    376 	if (6 == sscanf(build, "%hu:%hu:%hu/%3s %hu %hu",
    377 			&hour, &minute, &second, monstr, &day, &year)) {
    378 		cp = strstr(mlist, monstr);
    379 		if (NULL != cp) {
    380 			jd->year     = year;
    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;
    386 
    387 			return TRUE;
    388 		}
    389 	}
    390 
    391 	return FALSE;
    392 }
    393 
    394 
    395 /*
    396  *---------------------------------------------------------------------
    397  * basic calendar stuff
    398  * --------------------------------------------------------------------
    399  */
    400 
    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
    404 };
    405 
    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 }
    412 };
    413 
    414 /*
    415  * Some notes on the terminology:
    416  *
    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.
    422  *
    423  * Normally date parts are given as current cycles, while time parts
    424  * are given as elapsed cycles:
    425  *
    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.
    428  *
    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
    432  * that notation.
    433  *
    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.
    438  *
    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).
    443  */
    444 
    445 /*
    446  * ==================================================================
    447  *
    448  * General algorithmic stuff
    449  *
    450  * ==================================================================
    451  */
    452 
    453 /*
    454  *---------------------------------------------------------------------
    455  * Do a periodic extension of 'value' around 'pivot' with a period of
    456  * 'cycle'.
    457  *
    458  * The result 'res' is a number that holds to the following properties:
    459  *
    460  *   1)	 res MOD cycle == value MOD cycle
    461  *   2)	 pivot <= res < pivot + cycle
    462  *	 (replace </<= with >/>= for negative cycles)
    463  *
    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.
    469  *
    470  * This function has some useful applications:
    471  *
    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")
    477  *
    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
    482  *   unfolding")
    483  *
    484  * + Wherever you have a truncated periodic value and a non-truncated
    485  *   base value and you want to match them somehow...
    486  *
    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.
    491  *
    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  * ---------------------------------------------------------------------
    499  */
    500 int32_t
    501 ntpcal_periodic_extend(
    502 	int32_t pivot,
    503 	int32_t value,
    504 	int32_t cycle
    505 	)
    506 {
    507 	uint32_t diff;
    508 	char	 cpl = 0; /* modulo complement flag */
    509 	char	 neg = 0; /* sign change flag	    */
    510 
    511 	/* make the cycle positive and adjust the flags */
    512 	if (cycle < 0) {
    513 		cycle = - cycle;
    514 		neg ^= 1;
    515 		cpl ^= 1;
    516 	}
    517 	/* guard against div by zero or one */
    518 	if (cycle > 1) {
    519 		/*
    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
    523 		 * one.
    524 		 */
    525 		if (value >= pivot) {
    526 			diff = int32_to_uint32_2cpl(value)
    527 			     - int32_to_uint32_2cpl(pivot);
    528 		} else {
    529 			diff = int32_to_uint32_2cpl(pivot)
    530 			     - int32_to_uint32_2cpl(value);
    531 			cpl ^= 1;
    532 		}
    533 		diff %= (uint32_t)cycle;
    534 		if (diff) {
    535 			if (cpl)
    536 				diff = (uint32_t)cycle - diff;
    537 			if (neg)
    538 				diff = ~diff + 1;
    539 			pivot += uint32_2cpl_to_int32(diff);
    540 		}
    541 	}
    542 	return pivot;
    543 }
    544 
    545 /*
    546  *-------------------------------------------------------------------
    547  * Convert a timestamp in NTP scale to a 64bit seconds value in the UN*X
    548  * scale with proper epoch unfolding around a given pivot or the current
    549  * system time. This function happily accepts negative pivot values as
    550  * timestamps befor 1970-01-01, so be aware of possible trouble on
    551  * platforms with 32bit 'time_t'!
    552  *
    553  * This is also a periodic extension, but since the cycle is 2^32 and
    554  * the shift is 2^31, we can do some *very* fast math without explicit
    555  * divisions.
    556  *-------------------------------------------------------------------
    557  */
    558 vint64
    559 ntpcal_ntp_to_time(
    560 	uint32_t	ntp,
    561 	const time_t *	pivot
    562 	)
    563 {
    564 	vint64 res;
    565 
    566 #   if defined(HAVE_INT64)
    567 
    568 	res.q_s = (pivot != NULL)
    569 		      ? *pivot
    570 		      : now();
    571 	res.Q_s -= 0x80000000;		/* unshift of half range */
    572 	ntp	-= (uint32_t)JAN_1970;	/* warp into UN*X domain */
    573 	ntp	-= res.D_s.lo;		/* cycle difference	 */
    574 	res.Q_s += (uint64_t)ntp;	/* get expanded time	 */
    575 
    576 #   else /* no 64bit scalars */
    577 
    578 	time_t tmp;
    579 
    580 	tmp = (pivot != NULL)
    581 		  ? *pivot
    582 		  : now();
    583 	res = time_to_vint64(&tmp);
    584 	M_SUB(res.D_s.hi, res.D_s.lo, 0, 0x80000000);
    585 	ntp -= (uint32_t)JAN_1970;	/* warp into UN*X domain */
    586 	ntp -= res.D_s.lo;		/* cycle difference	 */
    587 	M_ADD(res.D_s.hi, res.D_s.lo, 0, ntp);
    588 
    589 #   endif /* no 64bit scalars */
    590 
    591 	return res;
    592 }
    593 
    594 /*
    595  *-------------------------------------------------------------------
    596  * Convert a timestamp in NTP scale to a 64bit seconds value in the NTP
    597  * scale with proper epoch unfolding around a given pivot or the current
    598  * system time.
    599  *
    600  * Note: The pivot must be given in the UN*X time domain!
    601  *
    602  * This is also a periodic extension, but since the cycle is 2^32 and
    603  * the shift is 2^31, we can do some *very* fast math without explicit
    604  * divisions.
    605  *-------------------------------------------------------------------
    606  */
    607 vint64
    608 ntpcal_ntp_to_ntp(
    609 	uint32_t      ntp,
    610 	const time_t *pivot
    611 	)
    612 {
    613 	vint64 res;
    614 
    615 #   if defined(HAVE_INT64)
    616 
    617 	res.q_s = (pivot)
    618 		      ? *pivot
    619 		      : now();
    620 	res.Q_s -= 0x80000000;		/* unshift of half range */
    621 	res.Q_s += (uint32_t)JAN_1970;	/* warp into NTP domain	 */
    622 	ntp	-= res.D_s.lo;		/* cycle difference	 */
    623 	res.Q_s += (uint64_t)ntp;	/* get expanded time	 */
    624 
    625 #   else /* no 64bit scalars */
    626 
    627 	time_t tmp;
    628 
    629 	tmp = (pivot)
    630 		  ? *pivot
    631 		  : now();
    632 	res = time_to_vint64(&tmp);
    633 	M_SUB(res.D_s.hi, res.D_s.lo, 0, 0x80000000u);
    634 	M_ADD(res.D_s.hi, res.D_s.lo, 0, (uint32_t)JAN_1970);/*into NTP */
    635 	ntp -= res.D_s.lo;		/* cycle difference	 */
    636 	M_ADD(res.D_s.hi, res.D_s.lo, 0, ntp);
    637 
    638 #   endif /* no 64bit scalars */
    639 
    640 	return res;
    641 }
    642 
    643 
    644 /*
    645  * ==================================================================
    646  *
    647  * Splitting values to composite entities
    648  *
    649  * ==================================================================
    650  */
    651 
    652 /*
    653  *-------------------------------------------------------------------
    654  * Split a 64bit seconds value into elapsed days in 'res.hi' and
    655  * elapsed seconds since midnight in 'res.lo' using explicit floor
    656  * division. This function happily accepts negative time values as
    657  * timestamps before the respective epoch start.
    658  * -------------------------------------------------------------------
    659  */
    660 ntpcal_split
    661 ntpcal_daysplit(
    662 	const vint64 *ts
    663 	)
    664 {
    665 	ntpcal_split res;
    666 	uint32_t Q;
    667 
    668 #   if defined(HAVE_INT64)
    669 
    670 	/* Manual floor division by SECSPERDAY. This uses the one's
    671 	 * complement trick, too, but without an extra flag value: The
    672 	 * flag would be 64bit, and that's a bit of overkill on a 32bit
    673 	 * target that has to use a register pair for a 64bit number.
    674 	 */
    675 	if (ts->q_s < 0)
    676 		Q = ~(uint32_t)(~ts->Q_s / SECSPERDAY);
    677 	else
    678 		Q = (uint32_t)(ts->Q_s / SECSPERDAY);
    679 
    680 #   else
    681 
    682 	uint32_t ah, al, sflag, A;
    683 
    684 	/* get operand into ah/al (either ts or ts' one's complement,
    685 	 * for later floor division)
    686 	 */
    687 	sflag = int32_sflag(ts->d_s.hi);
    688 	ah = sflag ^ ts->D_s.hi;
    689 	al = sflag ^ ts->D_s.lo;
    690 
    691 	/* Since 86400 == 128*675 we can drop the least 7 bits and
    692 	 * divide by 675 instead of 86400. Then the maximum remainder
    693 	 * after each devision step is 674, and we need 10 bits for
    694 	 * that. So in the next step we can shift in 22 bits from the
    695 	 * numerator.
    696 	 *
    697 	 * Therefore we load the accu with the top 13 bits (51..63) in
    698 	 * the first shot. We don't have to remember the quotient -- it
    699 	 * would be shifted out anyway.
    700 	 */
    701 	A = ah >> 19;
    702 	if (A >= 675)
    703 		A = (A % 675u);
    704 
    705 	/* Now assemble the remainder with bits 29..50 from the
    706 	 * numerator and divide. This creates the upper ten bits of the
    707 	 * quotient. (Well, the top 22 bits of a 44bit result. But that
    708 	 * will be truncated to 32 bits anyway.)
    709 	 */
    710 	A = (A << 19) | (ah & 0x0007FFFFu);
    711 	A = (A <<  3) | (al >> 29);
    712 	Q = A / 675u;
    713 	A = A % 675u;
    714 
    715 	/* Now assemble the remainder with bits 7..28 from the numerator
    716 	 * and do a final division step.
    717 	 */
    718 	A = (A << 22) | ((al >> 7) & 0x003FFFFFu);
    719 	Q = (Q << 22) | (A / 675u);
    720 
    721 	/* The last 7 bits get simply dropped, as they have no affect on
    722 	 * the quotient when dividing by 86400.
    723 	 */
    724 
    725 	/* apply sign correction and calculate the true floor
    726 	 * remainder.
    727 	 */
    728 	Q ^= sflag;
    729 
    730 #   endif
    731 
    732 	res.hi = uint32_2cpl_to_int32(Q);
    733 	res.lo = ts->D_s.lo - Q * SECSPERDAY;
    734 
    735 	return res;
    736 }
    737 
    738 /*
    739  *-------------------------------------------------------------------
    740  * Split a 32bit seconds value into h/m/s and excessive days.  This
    741  * function happily accepts negative time values as timestamps before
    742  * midnight.
    743  * -------------------------------------------------------------------
    744  */
    745 static int32_t
    746 priv_timesplit(
    747 	int32_t split[3],
    748 	int32_t ts
    749 	)
    750 {
    751 	/* Do 3 chained floor divisions by positive constants, using the
    752 	 * one's complement trick and factoring out the intermediate XOR
    753 	 * ops to reduce the number of operations.
    754 	 */
    755 	uint32_t us, um, uh, ud, sflag;
    756 
    757 	sflag = int32_sflag(ts);
    758 	us    = int32_to_uint32_2cpl(ts);
    759 
    760 	um = (sflag ^ us) / SECSPERMIN;
    761 	uh = um / MINSPERHR;
    762 	ud = uh / HRSPERDAY;
    763 
    764 	um ^= sflag;
    765 	uh ^= sflag;
    766 	ud ^= sflag;
    767 
    768 	split[0] = (int32_t)(uh - ud * HRSPERDAY );
    769 	split[1] = (int32_t)(um - uh * MINSPERHR );
    770 	split[2] = (int32_t)(us - um * SECSPERMIN);
    771 
    772 	return uint32_2cpl_to_int32(ud);
    773 }
    774 
    775 /*
    776  * ---------------------------------------------------------------------
    777  * Given the number of elapsed days in the calendar era, split this
    778  * number into the number of elapsed years in 'res.hi' and the number
    779  * of elapsed days of that year in 'res.lo'.
    780  *
    781  * if 'isleapyear' is not NULL, it will receive an integer that is 0 for
    782  * regular years and a non-zero value for leap years.
    783  *---------------------------------------------------------------------
    784  */
    785 ntpcal_split
    786 ntpcal_split_eradays(
    787 	int32_t days,
    788 	int  *isleapyear
    789 	)
    790 {
    791 	/* Use the fast cyclesplit algorithm here, to calculate the
    792 	 * centuries and years in a century with one division each. This
    793 	 * reduces the number of division operations to two, but is
    794 	 * susceptible to internal range overflow. We make sure the
    795 	 * input operands are in the safe range; this still gives us
    796 	 * approx +/-2.9 million years.
    797 	 */
    798 	ntpcal_split res;
    799 	int32_t	 n100, n001; /* calendar year cycles */
    800 	uint32_t uday, Q, sflag;
    801 
    802 	/* split off centuries first */
    803 	sflag = int32_sflag(days);
    804 	uday  = uint32_saturate(int32_to_uint32_2cpl(days), sflag);
    805 	uday  = (4u * uday) | 3u;
    806 	Q    = sflag ^ ((sflag ^ uday) / GREGORIAN_CYCLE_DAYS);
    807 	uday = uday - Q * GREGORIAN_CYCLE_DAYS;
    808 	n100 = uint32_2cpl_to_int32(Q);
    809 
    810 	/* Split off years in century -- days >= 0 here, and we're far
    811 	 * away from integer overflow trouble now. */
    812 	uday |= 3;
    813 	n001 = uday / GREGORIAN_NORMAL_LEAP_CYCLE_DAYS;
    814 	uday = uday % GREGORIAN_NORMAL_LEAP_CYCLE_DAYS;
    815 
    816 	/* Assemble the year and day in year */
    817 	res.hi = n100 * 100 + n001;
    818 	res.lo = uday / 4u;
    819 
    820 	/* Eventually set the leap year flag. Note: 0 <= n001 <= 99 and
    821 	 * Q is still the two's complement representation of the
    822 	 * centuries: The modulo 4 ops can be done with masking here.
    823 	 * We also shift the year and the century by one, so the tests
    824 	 * can be done against zero instead of 3.
    825 	 */
    826 	if (isleapyear)
    827 		*isleapyear = !((n001+1) & 3)
    828 		    && ((n001 != 99) || !((Q+1) & 3));
    829 
    830 	return res;
    831 }
    832 
    833 /*
    834  *---------------------------------------------------------------------
    835  * Given a number of elapsed days in a year and a leap year indicator,
    836  * split the number of elapsed days into the number of elapsed months in
    837  * 'res.hi' and the number of elapsed days of that month in 'res.lo'.
    838  *
    839  * This function will fail and return {-1,-1} if the number of elapsed
    840  * days is not in the valid range!
    841  *---------------------------------------------------------------------
    842  */
    843 ntpcal_split
    844 ntpcal_split_yeardays(
    845 	int32_t eyd,
    846 	int     isleapyear
    847 	)
    848 {
    849 	ntpcal_split    res;
    850 	const uint16_t *lt;	/* month length table	*/
    851 
    852 	/* check leap year flag and select proper table */
    853 	lt = real_month_table[(isleapyear != 0)];
    854 	if (0 <= eyd && eyd < lt[12]) {
    855 		/* get zero-based month by approximation & correction step */
    856 		res.hi = eyd >> 5;	   /* approx month; might be 1 too low */
    857 		if (lt[res.hi + 1] <= eyd) /* fixup approximative month value  */
    858 			res.hi += 1;
    859 		res.lo = eyd - lt[res.hi];
    860 	} else {
    861 		res.lo = res.hi = -1;
    862 	}
    863 
    864 	return res;
    865 }
    866 
    867 /*
    868  *---------------------------------------------------------------------
    869  * Convert a RD into the date part of a 'struct calendar'.
    870  *---------------------------------------------------------------------
    871  */
    872 int
    873 ntpcal_rd_to_date(
    874 	struct calendar *jd,
    875 	int32_t		 rd
    876 	)
    877 {
    878 	ntpcal_split split;
    879 	int	     leapy;
    880 	u_int	     ymask;
    881 
    882 	/* Get day-of-week first. Since rd is signed, the remainder can
    883 	 * be in the range [-6..+6], but the assignment to an unsigned
    884 	 * variable maps the negative values to positive values >=7.
    885 	 * This makes the sign correction look strange, but adding 7
    886 	 * causes the needed wrap-around into the desired value range of
    887 	 * zero to six, both inclusive.
    888 	 */
    889 	jd->weekday = rd % DAYSPERWEEK;
    890 	if (jd->weekday >= DAYSPERWEEK)	/* weekday is unsigned! */
    891 		jd->weekday += DAYSPERWEEK;
    892 
    893 	split = ntpcal_split_eradays(rd - 1, &leapy);
    894 	/* Get year and day-of-year, with overflow check. If any of the
    895 	 * upper 16 bits is set after shifting to unity-based years, we
    896 	 * will have an overflow when converting to an unsigned 16bit
    897 	 * year. Shifting to the right is OK here, since it does not
    898 	 * matter if the shift is logic or arithmetic.
    899 	 */
    900 	split.hi += 1;
    901 	ymask = 0u - ((split.hi >> 16) == 0);
    902 	jd->year = (uint16_t)(split.hi & ymask);
    903 	jd->yearday = (uint16_t)split.lo + 1;
    904 
    905 	/* convert to month and mday */
    906 	split = ntpcal_split_yeardays(split.lo, leapy);
    907 	jd->month    = (uint8_t)split.hi + 1;
    908 	jd->monthday = (uint8_t)split.lo + 1;
    909 
    910 	return ymask ? leapy : -1;
    911 }
    912 
    913 /*
    914  *---------------------------------------------------------------------
    915  * Convert a RD into the date part of a 'struct tm'.
    916  *---------------------------------------------------------------------
    917  */
    918 int
    919 ntpcal_rd_to_tm(
    920 	struct tm  *utm,
    921 	int32_t	    rd
    922 	)
    923 {
    924 	ntpcal_split split;
    925 	int	     leapy;
    926 
    927 	/* get day-of-week first */
    928 	utm->tm_wday = rd % DAYSPERWEEK;
    929 	if (utm->tm_wday < 0)
    930 		utm->tm_wday += DAYSPERWEEK;
    931 
    932 	/* get year and day-of-year */
    933 	split = ntpcal_split_eradays(rd - 1, &leapy);
    934 	utm->tm_year = split.hi - 1899;
    935 	utm->tm_yday = split.lo;	/* 0-based */
    936 
    937 	/* convert to month and mday */
    938 	split = ntpcal_split_yeardays(split.lo, leapy);
    939 	utm->tm_mon  = split.hi;	/* 0-based */
    940 	utm->tm_mday = split.lo + 1;	/* 1-based */
    941 
    942 	return leapy;
    943 }
    944 
    945 /*
    946  *---------------------------------------------------------------------
    947  * Take a value of seconds since midnight and split it into hhmmss in a
    948  * 'struct calendar'.
    949  *---------------------------------------------------------------------
    950  */
    951 int32_t
    952 ntpcal_daysec_to_date(
    953 	struct calendar *jd,
    954 	int32_t		sec
    955 	)
    956 {
    957 	int32_t days;
    958 	int   ts[3];
    959 
    960 	days = priv_timesplit(ts, sec);
    961 	jd->hour   = (uint8_t)ts[0];
    962 	jd->minute = (uint8_t)ts[1];
    963 	jd->second = (uint8_t)ts[2];
    964 
    965 	return days;
    966 }
    967 
    968 /*
    969  *---------------------------------------------------------------------
    970  * Take a value of seconds since midnight and split it into hhmmss in a
    971  * 'struct tm'.
    972  *---------------------------------------------------------------------
    973  */
    974 int32_t
    975 ntpcal_daysec_to_tm(
    976 	struct tm *utm,
    977 	int32_t	   sec
    978 	)
    979 {
    980 	int32_t days;
    981 	int32_t ts[3];
    982 
    983 	days = priv_timesplit(ts, sec);
    984 	utm->tm_hour = ts[0];
    985 	utm->tm_min  = ts[1];
    986 	utm->tm_sec  = ts[2];
    987 
    988 	return days;
    989 }
    990 
    991 /*
    992  *---------------------------------------------------------------------
    993  * take a split representation for day/second-of-day and day offset
    994  * and convert it to a 'struct calendar'. The seconds will be normalised
    995  * into the range of a day, and the day will be adjusted accordingly.
    996  *
    997  * returns >0 if the result is in a leap year, 0 if in a regular
    998  * year and <0 if the result did not fit into the calendar struct.
    999  *---------------------------------------------------------------------
   1000  */
   1001 int
   1002 ntpcal_daysplit_to_date(
   1003 	struct calendar	   *jd,
   1004 	const ntpcal_split *ds,
   1005 	int32_t		    dof
   1006 	)
   1007 {
   1008 	dof += ntpcal_daysec_to_date(jd, ds->lo);
   1009 	return ntpcal_rd_to_date(jd, ds->hi + dof);
   1010 }
   1011 
   1012 /*
   1013  *---------------------------------------------------------------------
   1014  * take a split representation for day/second-of-day and day offset
   1015  * and convert it to a 'struct tm'. The seconds will be normalised
   1016  * into the range of a day, and the day will be adjusted accordingly.
   1017  *
   1018  * returns 1 if the result is in a leap year and zero if in a regular
   1019  * year.
   1020  *---------------------------------------------------------------------
   1021  */
   1022 int
   1023 ntpcal_daysplit_to_tm(
   1024 	struct tm	   *utm,
   1025 	const ntpcal_split *ds ,
   1026 	int32_t		    dof
   1027 	)
   1028 {
   1029 	dof += ntpcal_daysec_to_tm(utm, ds->lo);
   1030 
   1031 	return ntpcal_rd_to_tm(utm, ds->hi + dof);
   1032 }
   1033 
   1034 /*
   1035  *---------------------------------------------------------------------
   1036  * Take a UN*X time and convert to a calendar structure.
   1037  *---------------------------------------------------------------------
   1038  */
   1039 int
   1040 ntpcal_time_to_date(
   1041 	struct calendar	*jd,
   1042 	const vint64	*ts
   1043 	)
   1044 {
   1045 	ntpcal_split ds;
   1046 
   1047 	ds = ntpcal_daysplit(ts);
   1048 	ds.hi += ntpcal_daysec_to_date(jd, ds.lo);
   1049 	ds.hi += DAY_UNIX_STARTS;
   1050 
   1051 	return ntpcal_rd_to_date(jd, ds.hi);
   1052 }
   1053 
   1054 
   1055 /*
   1056  * ==================================================================
   1057  *
   1058  * merging composite entities
   1059  *
   1060  * ==================================================================
   1061  */
   1062 
   1063 /*
   1064  *---------------------------------------------------------------------
   1065  * Merge a number of days and a number of seconds into seconds,
   1066  * expressed in 64 bits to avoid overflow.
   1067  *---------------------------------------------------------------------
   1068  */
   1069 vint64
   1070 ntpcal_dayjoin(
   1071 	int32_t days,
   1072 	int32_t secs
   1073 	)
   1074 {
   1075 	vint64 res;
   1076 
   1077 #   if defined(HAVE_INT64)
   1078 
   1079 	res.q_s	 = days;
   1080 	res.q_s *= SECSPERDAY;
   1081 	res.q_s += secs;
   1082 
   1083 #   else
   1084 
   1085 	uint32_t p1, p2;
   1086 	int	 isneg;
   1087 
   1088 	/*
   1089 	 * res = days *86400 + secs, using manual 16/32 bit
   1090 	 * multiplications and shifts.
   1091 	 */
   1092 	isneg = (days < 0);
   1093 	if (isneg)
   1094 		days = -days;
   1095 
   1096 	/* assemble days * 675 */
   1097 	res.D_s.lo = (days & 0xFFFF) * 675u;
   1098 	res.D_s.hi = 0;
   1099 	p1 = (days >> 16) * 675u;
   1100 	p2 = p1 >> 16;
   1101 	p1 = p1 << 16;
   1102 	M_ADD(res.D_s.hi, res.D_s.lo, p2, p1);
   1103 
   1104 	/* mul by 128, using shift */
   1105 	res.D_s.hi = (res.D_s.hi << 7) | (res.D_s.lo >> 25);
   1106 	res.D_s.lo = (res.D_s.lo << 7);
   1107 
   1108 	/* fix sign */
   1109 	if (isneg)
   1110 		M_NEG(res.D_s.hi, res.D_s.lo);
   1111 
   1112 	/* properly add seconds */
   1113 	p2 = 0;
   1114 	if (secs < 0) {
   1115 		p1 = (uint32_t)-secs;
   1116 		M_NEG(p2, p1);
   1117 	} else {
   1118 		p1 = (uint32_t)secs;
   1119 	}
   1120 	M_ADD(res.D_s.hi, res.D_s.lo, p2, p1);
   1121 
   1122 #   endif
   1123 
   1124 	return res;
   1125 }
   1126 
   1127 /*
   1128  *---------------------------------------------------------------------
   1129  * get leap years since epoch in elapsed years
   1130  *---------------------------------------------------------------------
   1131  */
   1132 int32_t
   1133 ntpcal_leapyears_in_years(
   1134 	int32_t years
   1135 	)
   1136 {
   1137 	/* We use the in-out-in algorithm here, using the one's
   1138 	 * complement division trick for negative numbers. The chained
   1139 	 * division sequence by 4/25/4 gives the compiler the chance to
   1140 	 * get away with only one true division and doing shifts otherwise.
   1141 	 */
   1142 
   1143 	uint32_t sflag, sum, uyear;
   1144 
   1145 	sflag = int32_sflag(years);
   1146 	uyear = int32_to_uint32_2cpl(years);
   1147 	uyear ^= sflag;
   1148 
   1149 	sum  = (uyear /=  4u);	/*   4yr rule --> IN  */
   1150 	sum -= (uyear /= 25u);	/* 100yr rule --> OUT */
   1151 	sum += (uyear /=  4u);	/* 400yr rule --> IN  */
   1152 
   1153 	/* Thanks to the alternation of IN/OUT/IN we can do the sum
   1154 	 * directly and have a single one's complement operation
   1155 	 * here. (Only if the years are negative, of course.) Otherwise
   1156 	 * the one's complement would have to be done when
   1157 	 * adding/subtracting the terms.
   1158 	 */
   1159 	return uint32_2cpl_to_int32(sflag ^ sum);
   1160 }
   1161 
   1162 /*
   1163  *---------------------------------------------------------------------
   1164  * Convert elapsed years in Era into elapsed days in Era.
   1165  *---------------------------------------------------------------------
   1166  */
   1167 int32_t
   1168 ntpcal_days_in_years(
   1169 	int32_t years
   1170 	)
   1171 {
   1172 	return years * DAYSPERYEAR + ntpcal_leapyears_in_years(years);
   1173 }
   1174 
   1175 /*
   1176  *---------------------------------------------------------------------
   1177  * Convert a number of elapsed month in a year into elapsed days in year.
   1178  *
   1179  * The month will be normalized, and 'res.hi' will contain the
   1180  * excessive years that must be considered when converting the years,
   1181  * while 'res.lo' will contain the number of elapsed days since start
   1182  * of the year.
   1183  *
   1184  * This code uses the shifted-month-approach to convert month to days,
   1185  * because then there is no need to have explicit leap year
   1186  * information.	 The slight disadvantage is that for most month values
   1187  * the result is a negative value, and the year excess is one; the
   1188  * conversion is then simply based on the start of the following year.
   1189  *---------------------------------------------------------------------
   1190  */
   1191 ntpcal_split
   1192 ntpcal_days_in_months(
   1193 	int32_t m
   1194 	)
   1195 {
   1196 	ntpcal_split res;
   1197 
   1198 	/* Add ten months and correct if needed. (It likely is...) */
   1199 	res.lo  = m + 10;
   1200 	res.hi  = (res.lo >= 12);
   1201 	if (res.hi)
   1202 		res.lo -= 12;
   1203 
   1204 	/* if still out of range, normalise by floor division ... */
   1205 	if (res.lo < 0 || res.lo >= 12) {
   1206 		uint32_t mu, Q, sflag;
   1207 		sflag = int32_sflag(res.lo);
   1208 		mu    = int32_to_uint32_2cpl(res.lo);
   1209 		Q     = sflag ^ ((sflag ^ mu) / 12u);
   1210 		res.hi += uint32_2cpl_to_int32(Q);
   1211 		res.lo  = mu - Q * 12u;
   1212 	}
   1213 
   1214 	/* get cummulated days in year with unshift */
   1215 	res.lo = shift_month_table[res.lo] - 306;
   1216 
   1217 	return res;
   1218 }
   1219 
   1220 /*
   1221  *---------------------------------------------------------------------
   1222  * Convert ELAPSED years/months/days of gregorian calendar to elapsed
   1223  * days in Gregorian epoch.
   1224  *
   1225  * If you want to convert years and days-of-year, just give a month of
   1226  * zero.
   1227  *---------------------------------------------------------------------
   1228  */
   1229 int32_t
   1230 ntpcal_edate_to_eradays(
   1231 	int32_t years,
   1232 	int32_t mons,
   1233 	int32_t mdays
   1234 	)
   1235 {
   1236 	ntpcal_split tmp;
   1237 	int32_t	     res;
   1238 
   1239 	if (mons) {
   1240 		tmp = ntpcal_days_in_months(mons);
   1241 		res = ntpcal_days_in_years(years + tmp.hi) + tmp.lo;
   1242 	} else
   1243 		res = ntpcal_days_in_years(years);
   1244 	res += mdays;
   1245 
   1246 	return res;
   1247 }
   1248 
   1249 /*
   1250  *---------------------------------------------------------------------
   1251  * Convert ELAPSED years/months/days of gregorian calendar to elapsed
   1252  * days in year.
   1253  *
   1254  * Note: This will give the true difference to the start of the given year,
   1255  * even if months & days are off-scale.
   1256  *---------------------------------------------------------------------
   1257  */
   1258 int32_t
   1259 ntpcal_edate_to_yeardays(
   1260 	int32_t years,
   1261 	int32_t mons,
   1262 	int32_t mdays
   1263 	)
   1264 {
   1265 	ntpcal_split tmp;
   1266 
   1267 	if (0 <= mons && mons < 12) {
   1268 		years += 1;
   1269 		mdays += real_month_table[is_leapyear(years)][mons];
   1270 	} else {
   1271 		tmp = ntpcal_days_in_months(mons);
   1272 		mdays += tmp.lo
   1273 		       + ntpcal_days_in_years(years + tmp.hi)
   1274 		       - ntpcal_days_in_years(years);
   1275 	}
   1276 
   1277 	return mdays;
   1278 }
   1279 
   1280 /*
   1281  *---------------------------------------------------------------------
   1282  * Convert elapsed days and the hour/minute/second information into
   1283  * total seconds.
   1284  *
   1285  * If 'isvalid' is not NULL, do a range check on the time specification
   1286  * and tell if the time input is in the normal range, permitting for a
   1287  * single leapsecond.
   1288  *---------------------------------------------------------------------
   1289  */
   1290 int32_t
   1291 ntpcal_etime_to_seconds(
   1292 	int32_t hours,
   1293 	int32_t minutes,
   1294 	int32_t seconds
   1295 	)
   1296 {
   1297 	int32_t res;
   1298 
   1299 	res = (hours * MINSPERHR + minutes) * SECSPERMIN + seconds;
   1300 
   1301 	return res;
   1302 }
   1303 
   1304 /*
   1305  *---------------------------------------------------------------------
   1306  * Convert the date part of a 'struct tm' (that is, year, month,
   1307  * day-of-month) into the RD of that day.
   1308  *---------------------------------------------------------------------
   1309  */
   1310 int32_t
   1311 ntpcal_tm_to_rd(
   1312 	const struct tm *utm
   1313 	)
   1314 {
   1315 	return ntpcal_edate_to_eradays(utm->tm_year + 1899,
   1316 				       utm->tm_mon,
   1317 				       utm->tm_mday - 1) + 1;
   1318 }
   1319 
   1320 /*
   1321  *---------------------------------------------------------------------
   1322  * Convert the date part of a 'struct calendar' (that is, year, month,
   1323  * day-of-month) into the RD of that day.
   1324  *---------------------------------------------------------------------
   1325  */
   1326 int32_t
   1327 ntpcal_date_to_rd(
   1328 	const struct calendar *jd
   1329 	)
   1330 {
   1331 	return ntpcal_edate_to_eradays((int32_t)jd->year - 1,
   1332 				       (int32_t)jd->month - 1,
   1333 				       (int32_t)jd->monthday - 1) + 1;
   1334 }
   1335 
   1336 /*
   1337  *---------------------------------------------------------------------
   1338  * convert a year number to rata die of year start
   1339  *---------------------------------------------------------------------
   1340  */
   1341 int32_t
   1342 ntpcal_year_to_ystart(
   1343 	int32_t year
   1344 	)
   1345 {
   1346 	return ntpcal_days_in_years(year - 1) + 1;
   1347 }
   1348 
   1349 /*
   1350  *---------------------------------------------------------------------
   1351  * For a given RD, get the RD of the associated year start,
   1352  * that is, the RD of the last January,1st on or before that day.
   1353  *---------------------------------------------------------------------
   1354  */
   1355 int32_t
   1356 ntpcal_rd_to_ystart(
   1357 	int32_t rd
   1358 	)
   1359 {
   1360 	/*
   1361 	 * Rather simple exercise: split the day number into elapsed
   1362 	 * years and elapsed days, then remove the elapsed days from the
   1363 	 * input value. Nice'n sweet...
   1364 	 */
   1365 	return rd - ntpcal_split_eradays(rd - 1, NULL).lo;
   1366 }
   1367 
   1368 /*
   1369  *---------------------------------------------------------------------
   1370  * For a given RD, get the RD of the associated month start.
   1371  *---------------------------------------------------------------------
   1372  */
   1373 int32_t
   1374 ntpcal_rd_to_mstart(
   1375 	int32_t rd
   1376 	)
   1377 {
   1378 	ntpcal_split split;
   1379 	int	     leaps;
   1380 
   1381 	split = ntpcal_split_eradays(rd - 1, &leaps);
   1382 	split = ntpcal_split_yeardays(split.lo, leaps);
   1383 
   1384 	return rd - split.lo;
   1385 }
   1386 
   1387 /*
   1388  *---------------------------------------------------------------------
   1389  * take a 'struct calendar' and get the seconds-of-day from it.
   1390  *---------------------------------------------------------------------
   1391  */
   1392 int32_t
   1393 ntpcal_date_to_daysec(
   1394 	const struct calendar *jd
   1395 	)
   1396 {
   1397 	return ntpcal_etime_to_seconds(jd->hour, jd->minute,
   1398 				       jd->second);
   1399 }
   1400 
   1401 /*
   1402  *---------------------------------------------------------------------
   1403  * take a 'struct tm' and get the seconds-of-day from it.
   1404  *---------------------------------------------------------------------
   1405  */
   1406 int32_t
   1407 ntpcal_tm_to_daysec(
   1408 	const struct tm *utm
   1409 	)
   1410 {
   1411 	return ntpcal_etime_to_seconds(utm->tm_hour, utm->tm_min,
   1412 				       utm->tm_sec);
   1413 }
   1414 
   1415 /*
   1416  *---------------------------------------------------------------------
   1417  * take a 'struct calendar' and convert it to a 'time_t'
   1418  *---------------------------------------------------------------------
   1419  */
   1420 time_t
   1421 ntpcal_date_to_time(
   1422 	const struct calendar *jd
   1423 	)
   1424 {
   1425 	vint64  join;
   1426 	int32_t days, secs;
   1427 
   1428 	days = ntpcal_date_to_rd(jd) - DAY_UNIX_STARTS;
   1429 	secs = ntpcal_date_to_daysec(jd);
   1430 	join = ntpcal_dayjoin(days, secs);
   1431 
   1432 	return vint64_to_time(&join);
   1433 }
   1434 
   1435 
   1436 /*
   1437  * ==================================================================
   1438  *
   1439  * extended and unchecked variants of caljulian/caltontp
   1440  *
   1441  * ==================================================================
   1442  */
   1443 int
   1444 ntpcal_ntp64_to_date(
   1445 	struct calendar *jd,
   1446 	const vint64    *ntp
   1447 	)
   1448 {
   1449 	ntpcal_split ds;
   1450 
   1451 	ds = ntpcal_daysplit(ntp);
   1452 	ds.hi += ntpcal_daysec_to_date(jd, ds.lo);
   1453 
   1454 	return ntpcal_rd_to_date(jd, ds.hi + DAY_NTP_STARTS);
   1455 }
   1456 
   1457 int
   1458 ntpcal_ntp_to_date(
   1459 	struct calendar *jd,
   1460 	uint32_t	 ntp,
   1461 	const time_t	*piv
   1462 	)
   1463 {
   1464 	vint64	ntp64;
   1465 
   1466 	/*
   1467 	 * Unfold ntp time around current time into NTP domain. Split
   1468 	 * into days and seconds, shift days into CE domain and
   1469 	 * process the parts.
   1470 	 */
   1471 	ntp64 = ntpcal_ntp_to_ntp(ntp, piv);
   1472 	return ntpcal_ntp64_to_date(jd, &ntp64);
   1473 }
   1474 
   1475 
   1476 vint64
   1477 ntpcal_date_to_ntp64(
   1478 	const struct calendar *jd
   1479 	)
   1480 {
   1481 	/*
   1482 	 * Convert date to NTP. Ignore yearday, use d/m/y only.
   1483 	 */
   1484 	return ntpcal_dayjoin(ntpcal_date_to_rd(jd) - DAY_NTP_STARTS,
   1485 			      ntpcal_date_to_daysec(jd));
   1486 }
   1487 
   1488 
   1489 uint32_t
   1490 ntpcal_date_to_ntp(
   1491 	const struct calendar *jd
   1492 	)
   1493 {
   1494 	/*
   1495 	 * Get lower half of 64-bit NTP timestamp from date/time.
   1496 	 */
   1497 	return ntpcal_date_to_ntp64(jd).d_s.lo;
   1498 }
   1499 
   1500 
   1501 
   1502 /*
   1503  * ==================================================================
   1504  *
   1505  * day-of-week calculations
   1506  *
   1507  * ==================================================================
   1508  */
   1509 /*
   1510  * Given a RataDie and a day-of-week, calculate a RDN that is reater-than,
   1511  * greater-or equal, closest, less-or-equal or less-than the given RDN
   1512  * and denotes the given day-of-week
   1513  */
   1514 int32_t
   1515 ntpcal_weekday_gt(
   1516 	int32_t rdn,
   1517 	int32_t dow
   1518 	)
   1519 {
   1520 	return ntpcal_periodic_extend(rdn+1, dow, 7);
   1521 }
   1522 
   1523 int32_t
   1524 ntpcal_weekday_ge(
   1525 	int32_t rdn,
   1526 	int32_t dow
   1527 	)
   1528 {
   1529 	return ntpcal_periodic_extend(rdn, dow, 7);
   1530 }
   1531 
   1532 int32_t
   1533 ntpcal_weekday_close(
   1534 	int32_t rdn,
   1535 	int32_t dow
   1536 	)
   1537 {
   1538 	return ntpcal_periodic_extend(rdn-3, dow, 7);
   1539 }
   1540 
   1541 int32_t
   1542 ntpcal_weekday_le(
   1543 	int32_t rdn,
   1544 	int32_t dow
   1545 	)
   1546 {
   1547 	return ntpcal_periodic_extend(rdn, dow, -7);
   1548 }
   1549 
   1550 int32_t
   1551 ntpcal_weekday_lt(
   1552 	int32_t rdn,
   1553 	int32_t dow
   1554 	)
   1555 {
   1556 	return ntpcal_periodic_extend(rdn-1, dow, -7);
   1557 }
   1558 
   1559 /*
   1560  * ==================================================================
   1561  *
   1562  * ISO week-calendar conversions
   1563  *
   1564  * The ISO8601 calendar defines a calendar of years, weeks and weekdays.
   1565  * It is related to the Gregorian calendar, and a ISO year starts at the
   1566  * Monday closest to Jan,1st of the corresponding Gregorian year.  A ISO
   1567  * calendar year has always 52 or 53 weeks, and like the Grogrian
   1568  * calendar the ISO8601 calendar repeats itself every 400 years, or
   1569  * 146097 days, or 20871 weeks.
   1570  *
   1571  * While it is possible to write ISO calendar functions based on the
   1572  * Gregorian calendar functions, the following implementation takes a
   1573  * different approach, based directly on years and weeks.
   1574  *
   1575  * Analysis of the tabulated data shows that it is not possible to
   1576  * interpolate from years to weeks over a full 400 year range; cyclic
   1577  * shifts over 400 years do not provide a solution here. But it *is*
   1578  * possible to interpolate over every single century of the 400-year
   1579  * cycle. (The centennial leap year rule seems to be the culprit here.)
   1580  *
   1581  * It can be shown that a conversion from years to weeks can be done
   1582  * using a linear transformation of the form
   1583  *
   1584  *   w = floor( y * a + b )
   1585  *
   1586  * where the slope a must hold to
   1587  *
   1588  *  52.1780821918 <= a < 52.1791044776
   1589  *
   1590  * and b must be chosen according to the selected slope and the number
   1591  * of the century in a 400-year period.
   1592  *
   1593  * The inverse calculation can also be done in this way. Careful scaling
   1594  * provides an unlimited set of integer coefficients a,k,b that enable
   1595  * us to write the calulation in the form
   1596  *
   1597  *   w = (y * a	 + b ) / k
   1598  *   y = (w * a' + b') / k'
   1599  *
   1600  * In this implementation the values of k and k' are chosen to be
   1601  * smallest possible powers of two, so the division can be implemented
   1602  * as shifts if the optimiser chooses to do so.
   1603  *
   1604  * ==================================================================
   1605  */
   1606 
   1607 /*
   1608  * Given a number of elapsed (ISO-)years since the begin of the
   1609  * christian era, return the number of elapsed weeks corresponding to
   1610  * the number of years.
   1611  */
   1612 int32_t
   1613 isocal_weeks_in_years(
   1614 	int32_t years
   1615 	)
   1616 {
   1617 	/*
   1618 	 * use: w = (y * 53431 + b[c]) / 1024 as interpolation
   1619 	 */
   1620 	static const uint16_t bctab[4] = { 157, 449, 597, 889 };
   1621 
   1622 	int32_t  cs, cw;
   1623 	uint32_t cc, ci, yu, sflag;
   1624 
   1625 	sflag = int32_sflag(years);
   1626 	yu    = int32_to_uint32_2cpl(years);
   1627 
   1628 	/* split off centuries, using floor division */
   1629 	cc  = sflag ^ ((sflag ^ yu) / 100u);
   1630 	yu -= cc * 100u;
   1631 
   1632 	/* calculate century cycles shift and cycle index:
   1633 	 * Assuming a century is 5217 weeks, we have to add a cycle
   1634 	 * shift that is 3 for every 4 centuries, because 3 of the four
   1635 	 * centuries have 5218 weeks. So '(cc*3 + 1) / 4' is the actual
   1636 	 * correction, and the second century is the defective one.
   1637 	 *
   1638 	 * Needs floor division by 4, which is done with masking and
   1639 	 * shifting.
   1640 	 */
   1641 	ci = cc * 3u + 1;
   1642 	cs = uint32_2cpl_to_int32(sflag ^ ((sflag ^ ci) / 4u));
   1643 	ci = ci % 4u;
   1644 
   1645 	/* Get weeks in century. Can use plain division here as all ops
   1646 	 * are >= 0,  and let the compiler sort out the possible
   1647 	 * optimisations.
   1648 	 */
   1649 	cw = (yu * 53431u + bctab[ci]) / 1024u;
   1650 
   1651 	return uint32_2cpl_to_int32(cc) * 5217 + cs + cw;
   1652 }
   1653 
   1654 /*
   1655  * Given a number of elapsed weeks since the begin of the christian
   1656  * era, split this number into the number of elapsed years in res.hi
   1657  * and the excessive number of weeks in res.lo. (That is, res.lo is
   1658  * the number of elapsed weeks in the remaining partial year.)
   1659  */
   1660 ntpcal_split
   1661 isocal_split_eraweeks(
   1662 	int32_t weeks
   1663 	)
   1664 {
   1665 	/*
   1666 	 * use: y = (w * 157 + b[c]) / 8192 as interpolation
   1667 	 */
   1668 
   1669 	static const uint16_t bctab[4] = { 85, 130, 17, 62 };
   1670 
   1671 	ntpcal_split res;
   1672 	int32_t  cc, ci;
   1673 	uint32_t sw, cy, Q, sflag;
   1674 
   1675 	/* Use two fast cycle-split divisions here. This is again
   1676 	 * susceptible to internal overflow, so we check the range. This
   1677 	 * still permits more than +/-20 million years, so this is
   1678 	 * likely a pure academical problem.
   1679 	 *
   1680 	 * We want to execute '(weeks * 4 + 2) /% 20871' under floor
   1681 	 * division rules in the first step.
   1682 	 */
   1683 	sflag = int32_sflag(weeks);
   1684 	sw  = uint32_saturate(int32_to_uint32_2cpl(weeks), sflag);
   1685 	sw  = 4u * sw + 2;
   1686 	Q   = sflag ^ ((sflag ^ sw) / GREGORIAN_CYCLE_WEEKS);
   1687 	sw -= Q * GREGORIAN_CYCLE_WEEKS;
   1688 	ci  = Q % 4u;
   1689 	cc  = uint32_2cpl_to_int32(Q);
   1690 
   1691 	/* Split off years; sw >= 0 here! The scaled weeks in the years
   1692 	 * are scaled up by 157 afterwards.
   1693 	 */
   1694 	sw  = (sw / 4u) * 157u + bctab[ci];
   1695 	cy  = sw / 8192u;	/* ws >> 13 , let the compiler sort it out */
   1696 	sw  = sw % 8192u;	/* ws & 8191, let the compiler sort it out */
   1697 
   1698 	/* assemble elapsed years and downscale the elapsed weeks in
   1699 	 * the year.
   1700 	 */
   1701 	res.hi = 100*cc + cy;
   1702 	res.lo = sw / 157u;
   1703 
   1704 	return res;
   1705 }
   1706 
   1707 /*
   1708  * Given a second in the NTP time scale and a pivot, expand the NTP
   1709  * time stamp around the pivot and convert into an ISO calendar time
   1710  * stamp.
   1711  */
   1712 int
   1713 isocal_ntp64_to_date(
   1714 	struct isodate *id,
   1715 	const vint64   *ntp
   1716 	)
   1717 {
   1718 	ntpcal_split ds;
   1719 	int32_t      ts[3];
   1720 	uint32_t     uw, ud, sflag;
   1721 
   1722 	/*
   1723 	 * Split NTP time into days and seconds, shift days into CE
   1724 	 * domain and process the parts.
   1725 	 */
   1726 	ds = ntpcal_daysplit(ntp);
   1727 
   1728 	/* split time part */
   1729 	ds.hi += priv_timesplit(ts, ds.lo);
   1730 	id->hour   = (uint8_t)ts[0];
   1731 	id->minute = (uint8_t)ts[1];
   1732 	id->second = (uint8_t)ts[2];
   1733 
   1734 	/* split days into days and weeks, using floor division in unsigned */
   1735 	ds.hi += DAY_NTP_STARTS - 1; /* shift from NTP to RDN */
   1736 	sflag = int32_sflag(ds.hi);
   1737 	ud  = int32_to_uint32_2cpl(ds.hi);
   1738 	uw  = sflag ^ ((sflag ^ ud) / DAYSPERWEEK);
   1739 	ud -= uw * DAYSPERWEEK;
   1740 	ds.hi = uint32_2cpl_to_int32(uw);
   1741 	ds.lo = ud;
   1742 
   1743 	id->weekday = (uint8_t)ds.lo + 1;	/* weekday result    */
   1744 
   1745 	/* get year and week in year */
   1746 	ds = isocal_split_eraweeks(ds.hi);	/* elapsed years&week*/
   1747 	id->year = (uint16_t)ds.hi + 1;		/* shift to current  */
   1748 	id->week = (uint8_t )ds.lo + 1;
   1749 
   1750 	return (ds.hi >= 0 && ds.hi < 0x0000FFFF);
   1751 }
   1752 
   1753 int
   1754 isocal_ntp_to_date(
   1755 	struct isodate *id,
   1756 	uint32_t	ntp,
   1757 	const time_t   *piv
   1758 	)
   1759 {
   1760 	vint64	ntp64;
   1761 
   1762 	/*
   1763 	 * Unfold ntp time around current time into NTP domain, then
   1764 	 * convert the full time stamp.
   1765 	 */
   1766 	ntp64 = ntpcal_ntp_to_ntp(ntp, piv);
   1767 	return isocal_ntp64_to_date(id, &ntp64);
   1768 }
   1769 
   1770 /*
   1771  * Convert a ISO date spec into a second in the NTP time scale,
   1772  * properly truncated to 32 bit.
   1773  */
   1774 vint64
   1775 isocal_date_to_ntp64(
   1776 	const struct isodate *id
   1777 	)
   1778 {
   1779 	int32_t weeks, days, secs;
   1780 
   1781 	weeks = isocal_weeks_in_years((int32_t)id->year - 1)
   1782 	      + (int32_t)id->week - 1;
   1783 	days = weeks * 7 + (int32_t)id->weekday;
   1784 	/* days is RDN of ISO date now */
   1785 	secs = ntpcal_etime_to_seconds(id->hour, id->minute, id->second);
   1786 
   1787 	return ntpcal_dayjoin(days - DAY_NTP_STARTS, secs);
   1788 }
   1789 
   1790 uint32_t
   1791 isocal_date_to_ntp(
   1792 	const struct isodate *id
   1793 	)
   1794 {
   1795 	/*
   1796 	 * Get lower half of 64-bit NTP timestamp from date/time.
   1797 	 */
   1798 	return isocal_date_to_ntp64(id).d_s.lo;
   1799 }
   1800 
   1801 /* -*-EOF-*- */
   1802