ntp_calendar.c revision 289997
1/*
2 * ntp_calendar.c - calendar and helper functions
3 *
4 * Written by Juergen Perlinger (perlinger@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
97static systime_func_ptr systime_func = &time;
98static inline time_t now(void);
99
100
101systime_func_ptr
102ntpcal_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
117static inline time_t
118now(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
129static inline uint32_t
130int32_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
153static inline uint32_t
154int32_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
185static inline int32_t
186uint32_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 */
221static inline uint32_t
222uint32_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 */
239vint64
240time_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
285time_t
286vint64_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 */
314int
315ntpcal_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 */
402static 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 */
407static 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 */
500int32_t
501ntpcal_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 */
558vint64
559ntpcal_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 */
607vint64
608ntpcal_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 */
660ntpcal_split
661ntpcal_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 */
745static int32_t
746priv_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 */
785ntpcal_split
786ntpcal_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 */
843ntpcal_split
844ntpcal_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 */
872int
873ntpcal_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 */
918int
919ntpcal_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 */
951int32_t
952ntpcal_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 */
974int32_t
975ntpcal_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 */
1001int
1002ntpcal_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 */
1022int
1023ntpcal_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 */
1039int
1040ntpcal_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 */
1069vint64
1070ntpcal_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 */
1132int32_t
1133ntpcal_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 */
1167int32_t
1168ntpcal_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 */
1191ntpcal_split
1192ntpcal_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 */
1229int32_t
1230ntpcal_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 */
1258int32_t
1259ntpcal_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 */
1290int32_t
1291ntpcal_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 */
1310int32_t
1311ntpcal_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 */
1326int32_t
1327ntpcal_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 */
1341int32_t
1342ntpcal_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 */
1355int32_t
1356ntpcal_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 */
1373int32_t
1374ntpcal_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 */
1392int32_t
1393ntpcal_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 */
1406int32_t
1407ntpcal_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 */
1420time_t
1421ntpcal_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 */
1443int
1444ntpcal_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
1457int
1458ntpcal_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
1476vint64
1477ntpcal_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
1489uint32_t
1490ntpcal_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 */
1514int32_t
1515ntpcal_weekday_gt(
1516	int32_t rdn,
1517	int32_t dow
1518	)
1519{
1520	return ntpcal_periodic_extend(rdn+1, dow, 7);
1521}
1522
1523int32_t
1524ntpcal_weekday_ge(
1525	int32_t rdn,
1526	int32_t dow
1527	)
1528{
1529	return ntpcal_periodic_extend(rdn, dow, 7);
1530}
1531
1532int32_t
1533ntpcal_weekday_close(
1534	int32_t rdn,
1535	int32_t dow
1536	)
1537{
1538	return ntpcal_periodic_extend(rdn-3, dow, 7);
1539}
1540
1541int32_t
1542ntpcal_weekday_le(
1543	int32_t rdn,
1544	int32_t dow
1545	)
1546{
1547	return ntpcal_periodic_extend(rdn, dow, -7);
1548}
1549
1550int32_t
1551ntpcal_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 */
1612int32_t
1613isocal_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 */
1660ntpcal_split
1661isocal_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 */
1712int
1713isocal_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
1753int
1754isocal_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 */
1774vint64
1775isocal_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
1790uint32_t
1791isocal_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