ntp_calendar.c revision 290001
1290001Sglebius/*
2290001Sglebius * ntp_calendar.c - calendar and helper functions
3290001Sglebius *
4290001Sglebius * Written by Juergen Perlinger (perlinger@ntp.org) for the NTP project.
5290001Sglebius * The contents of 'html/copyright.html' apply.
6290001Sglebius *
7290001Sglebius * --------------------------------------------------------------------
8290001Sglebius * Some notes on the implementation:
9290001Sglebius *
10290001Sglebius * Calendar algorithms thrive on the division operation, which is one of
11290001Sglebius * the slowest numerical operations in any CPU. What saves us here from
12290001Sglebius * abysmal performance is the fact that all divisions are divisions by
13290001Sglebius * constant numbers, and most compilers can do this by a multiplication
14290001Sglebius * operation.  But this might not work when using the div/ldiv/lldiv
15290001Sglebius * function family, because many compilers are not able to do inline
16290001Sglebius * expansion of the code with following optimisation for the
17290001Sglebius * constant-divider case.
18290001Sglebius *
19290001Sglebius * Also div/ldiv/lldiv are defined in terms of int/long/longlong, which
20290001Sglebius * are inherently target dependent. Nothing that could not be cured with
21290001Sglebius * autoconf, but still a mess...
22290001Sglebius *
23290001Sglebius * Furthermore, we need floor division in many places. C either leaves
24290001Sglebius * the division behaviour undefined (< C99) or demands truncation to
25290001Sglebius * zero (>= C99), so additional steps are required to make sure the
26290001Sglebius * algorithms work. The {l,ll}div function family is requested to
27290001Sglebius * truncate towards zero, which is also the wrong direction for our
28290001Sglebius * purpose.
29290001Sglebius *
30290001Sglebius * For all this, all divisions by constant are coded manually, even when
31290001Sglebius * there is a joined div/mod operation: The optimiser should sort that
32290001Sglebius * out, if possible. Most of the calculations are done with unsigned
33290001Sglebius * types, explicitely using two's complement arithmetics where
34290001Sglebius * necessary. This minimises the dependecies to compiler and target,
35290001Sglebius * while still giving reasonable to good performance.
36290001Sglebius *
37290001Sglebius * The implementation uses a few tricks that exploit properties of the
38290001Sglebius * two's complement: Floor division on negative dividents can be
39290001Sglebius * executed by using the one's complement of the divident. One's
40290001Sglebius * complement can be easily created using XOR and a mask.
41290001Sglebius *
42290001Sglebius * Finally, check for overflow conditions is minimal. There are only two
43290001Sglebius * calculation steps in the whole calendar that suffer from an internal
44290001Sglebius * overflow, and these conditions are checked: errno is set to EDOM and
45290001Sglebius * the results are clamped/saturated in this case.  All other functions
46290001Sglebius * do not suffer from internal overflow and simply return the result
47290001Sglebius * truncated to 32 bits.
48290001Sglebius *
49290001Sglebius * This is a sacrifice made for execution speed.  Since a 32-bit day
50290001Sglebius * counter covers +/- 5,879,610 years and the clamp limits the effective
51290001Sglebius * range to +/-2.9 million years, this should not pose a problem here.
52290001Sglebius *
53290001Sglebius */
54290001Sglebius
55290001Sglebius#include <config.h>
56290001Sglebius#include <sys/types.h>
57290001Sglebius
58290001Sglebius#include "ntp_types.h"
59290001Sglebius#include "ntp_calendar.h"
60290001Sglebius#include "ntp_stdlib.h"
61290001Sglebius#include "ntp_fp.h"
62290001Sglebius#include "ntp_unixtime.h"
63290001Sglebius
64290001Sglebius/* For now, let's take the conservative approach: if the target property
65290001Sglebius * macros are not defined, check a few well-known compiler/architecture
66290001Sglebius * settings. Default is to assume that the representation of signed
67290001Sglebius * integers is unknown and shift-arithmetic-right is not available.
68290001Sglebius */
69290001Sglebius#ifndef TARGET_HAS_2CPL
70290001Sglebius# if defined(__GNUC__)
71290001Sglebius#  if defined(__i386__) || defined(__x86_64__) || defined(__arm__)
72290001Sglebius#   define TARGET_HAS_2CPL 1
73290001Sglebius#  else
74290001Sglebius#   define TARGET_HAS_2CPL 0
75290001Sglebius#  endif
76290001Sglebius# elif defined(_MSC_VER)
77290001Sglebius#  if defined(_M_IX86) || defined(_M_X64) || defined(_M_ARM)
78290001Sglebius#   define TARGET_HAS_2CPL 1
79290001Sglebius#  else
80290001Sglebius#   define TARGET_HAS_2CPL 0
81290001Sglebius#  endif
82290001Sglebius# else
83290001Sglebius#  define TARGET_HAS_2CPL 0
84290001Sglebius# endif
85290001Sglebius#endif
86290001Sglebius
87290001Sglebius#ifndef TARGET_HAS_SAR
88290001Sglebius# define TARGET_HAS_SAR 0
89290001Sglebius#endif
90290001Sglebius
91290001Sglebius/*
92290001Sglebius *---------------------------------------------------------------------
93290001Sglebius * replacing the 'time()' function
94290001Sglebius * --------------------------------------------------------------------
95290001Sglebius */
96290001Sglebius
97290001Sglebiusstatic systime_func_ptr systime_func = &time;
98290001Sglebiusstatic inline time_t now(void);
99290001Sglebius
100290001Sglebius
101290001Sglebiussystime_func_ptr
102290001Sglebiusntpcal_set_timefunc(
103290001Sglebius	systime_func_ptr nfunc
104290001Sglebius	)
105290001Sglebius{
106290001Sglebius	systime_func_ptr res;
107290001Sglebius
108290001Sglebius	res = systime_func;
109290001Sglebius	if (NULL == nfunc)
110290001Sglebius		nfunc = &time;
111290001Sglebius	systime_func = nfunc;
112290001Sglebius
113290001Sglebius	return res;
114290001Sglebius}
115290001Sglebius
116290001Sglebius
117290001Sglebiusstatic inline time_t
118290001Sglebiusnow(void)
119290001Sglebius{
120290001Sglebius	return (*systime_func)(NULL);
121290001Sglebius}
122290001Sglebius
123290001Sglebius/*
124290001Sglebius *---------------------------------------------------------------------
125290001Sglebius * Get sign extension mask and unsigned 2cpl rep for a signed integer
126290001Sglebius *---------------------------------------------------------------------
127290001Sglebius */
128290001Sglebius
129290001Sglebiusstatic inline uint32_t
130290001Sglebiusint32_sflag(
131290001Sglebius	const int32_t v)
132290001Sglebius{
133290001Sglebius#   if TARGET_HAS_2CPL && TARGET_HAS_SAR && SIZEOF_INT >= 4
134290001Sglebius
135290001Sglebius	/* Let's assume that shift is the fastest way to get the sign
136290001Sglebius	 * extension of of a signed integer. This might not always be
137290001Sglebius	 * true, though -- On 8bit CPUs or machines without barrel
138290001Sglebius	 * shifter this will kill the performance. So we make sure
139290001Sglebius	 * we do this only if 'int' has at least 4 bytes.
140290001Sglebius	 */
141290001Sglebius	return (uint32_t)(v >> 31);
142290001Sglebius
143290001Sglebius#   else
144290001Sglebius
145290001Sglebius	/* This should be a rather generic approach for getting a sign
146290001Sglebius	 * extension mask...
147290001Sglebius	 */
148290001Sglebius	return UINT32_C(0) - (uint32_t)(v < 0);
149290001Sglebius
150290001Sglebius#   endif
151290001Sglebius}
152290001Sglebius
153290001Sglebiusstatic inline uint32_t
154290001Sglebiusint32_to_uint32_2cpl(
155290001Sglebius	const int32_t v)
156290001Sglebius{
157290001Sglebius	uint32_t vu;
158290001Sglebius
159290001Sglebius#   if TARGET_HAS_2CPL
160290001Sglebius
161290001Sglebius	/* Just copy through the 32 bits from the signed value if we're
162290001Sglebius	 * on a two's complement target.
163290001Sglebius	 */
164290001Sglebius	vu = (uint32_t)v;
165290001Sglebius
166290001Sglebius#   else
167290001Sglebius
168290001Sglebius	/* Convert from signed int to unsigned int two's complement. Do
169290001Sglebius	 * not make any assumptions about the representation of signed
170290001Sglebius	 * integers, but make sure signed integer overflow cannot happen
171290001Sglebius	 * here. A compiler on a two's complement target *might* find
172290001Sglebius	 * out that this is just a complicated cast (as above), but your
173290001Sglebius	 * mileage might vary.
174290001Sglebius	 */
175290001Sglebius	if (v < 0)
176290001Sglebius		vu = ~(uint32_t)(-(v + 1));
177290001Sglebius	else
178290001Sglebius		vu = (uint32_t)v;
179290001Sglebius
180290001Sglebius#   endif
181290001Sglebius
182290001Sglebius	return vu;
183290001Sglebius}
184290001Sglebius
185290001Sglebiusstatic inline int32_t
186290001Sglebiusuint32_2cpl_to_int32(
187290001Sglebius	const uint32_t vu)
188290001Sglebius{
189290001Sglebius	int32_t v;
190290001Sglebius
191290001Sglebius#   if TARGET_HAS_2CPL
192290001Sglebius
193290001Sglebius	/* Just copy through the 32 bits from the unsigned value if
194290001Sglebius	 * we're on a two's complement target.
195290001Sglebius	 */
196290001Sglebius	v = (int32_t)vu;
197290001Sglebius
198290001Sglebius#   else
199290001Sglebius
200290001Sglebius	/* Convert to signed integer, making sure signed integer
201290001Sglebius	 * overflow cannot happen. Again, the optimiser might or might
202290001Sglebius	 * not find out that this is just a copy of 32 bits on a target
203290001Sglebius	 * with two's complement representation for signed integers.
204290001Sglebius	 */
205290001Sglebius	if (vu > INT32_MAX)
206290001Sglebius		v = -(int32_t)(~vu) - 1;
207290001Sglebius	else
208290001Sglebius		v = (int32_t)vu;
209290001Sglebius
210290001Sglebius#   endif
211290001Sglebius
212290001Sglebius	return v;
213290001Sglebius}
214290001Sglebius
215290001Sglebius/* Some of the calculations need to multiply the input by 4 before doing
216290001Sglebius * a division. This can cause overflow and strange results. Therefore we
217290001Sglebius * clamp / saturate the input operand. And since we do the calculations
218290001Sglebius * in unsigned int with an extra sign flag/mask, we only loose one bit
219290001Sglebius * of the input value range.
220290001Sglebius */
221290001Sglebiusstatic inline uint32_t
222290001Sglebiusuint32_saturate(
223290001Sglebius	uint32_t vu,
224290001Sglebius	uint32_t mu)
225290001Sglebius{
226290001Sglebius	static const uint32_t limit = UINT32_MAX/4u;
227290001Sglebius	if ((mu ^ vu) > limit) {
228290001Sglebius		vu    = mu ^ limit;
229290001Sglebius		errno = EDOM;
230290001Sglebius	}
231290001Sglebius	return vu;
232290001Sglebius}
233290001Sglebius
234290001Sglebius/*
235290001Sglebius *---------------------------------------------------------------------
236290001Sglebius * Convert between 'time_t' and 'vint64'
237290001Sglebius *---------------------------------------------------------------------
238290001Sglebius */
239290001Sglebiusvint64
240290001Sglebiustime_to_vint64(
241290001Sglebius	const time_t * ptt
242290001Sglebius	)
243290001Sglebius{
244290001Sglebius	vint64 res;
245290001Sglebius	time_t tt;
246290001Sglebius
247290001Sglebius	tt = *ptt;
248290001Sglebius
249290001Sglebius#   if SIZEOF_TIME_T <= 4
250290001Sglebius
251290001Sglebius	res.D_s.hi = 0;
252290001Sglebius	if (tt < 0) {
253290001Sglebius		res.D_s.lo = (uint32_t)-tt;
254290001Sglebius		M_NEG(res.D_s.hi, res.D_s.lo);
255290001Sglebius	} else {
256290001Sglebius		res.D_s.lo = (uint32_t)tt;
257290001Sglebius	}
258290001Sglebius
259290001Sglebius#   elif defined(HAVE_INT64)
260290001Sglebius
261290001Sglebius	res.q_s = tt;
262290001Sglebius
263290001Sglebius#   else
264290001Sglebius	/*
265290001Sglebius	 * shifting negative signed quantities is compiler-dependent, so
266290001Sglebius	 * we better avoid it and do it all manually. And shifting more
267290001Sglebius	 * than the width of a quantity is undefined. Also a don't do!
268290001Sglebius	 */
269290001Sglebius	if (tt < 0) {
270290001Sglebius		tt = -tt;
271290001Sglebius		res.D_s.lo = (uint32_t)tt;
272290001Sglebius		res.D_s.hi = (uint32_t)(tt >> 32);
273290001Sglebius		M_NEG(res.D_s.hi, res.D_s.lo);
274290001Sglebius	} else {
275290001Sglebius		res.D_s.lo = (uint32_t)tt;
276290001Sglebius		res.D_s.hi = (uint32_t)(tt >> 32);
277290001Sglebius	}
278290001Sglebius
279290001Sglebius#   endif
280290001Sglebius
281290001Sglebius	return res;
282290001Sglebius}
283290001Sglebius
284290001Sglebius
285290001Sglebiustime_t
286290001Sglebiusvint64_to_time(
287290001Sglebius	const vint64 *tv
288290001Sglebius	)
289290001Sglebius{
290290001Sglebius	time_t res;
291290001Sglebius
292290001Sglebius#   if SIZEOF_TIME_T <= 4
293290001Sglebius
294290001Sglebius	res = (time_t)tv->D_s.lo;
295290001Sglebius
296290001Sglebius#   elif defined(HAVE_INT64)
297290001Sglebius
298290001Sglebius	res = (time_t)tv->q_s;
299290001Sglebius
300290001Sglebius#   else
301290001Sglebius
302290001Sglebius	res = ((time_t)tv->d_s.hi << 32) | tv->D_s.lo;
303290001Sglebius
304290001Sglebius#   endif
305290001Sglebius
306290001Sglebius	return res;
307290001Sglebius}
308290001Sglebius
309290001Sglebius/*
310290001Sglebius *---------------------------------------------------------------------
311290001Sglebius * Get the build date & time
312290001Sglebius *---------------------------------------------------------------------
313290001Sglebius */
314290001Sglebiusint
315290001Sglebiusntpcal_get_build_date(
316290001Sglebius	struct calendar * jd
317290001Sglebius	)
318290001Sglebius{
319290001Sglebius	/* The C standard tells us the format of '__DATE__':
320290001Sglebius	 *
321290001Sglebius	 * __DATE__ The date of translation of the preprocessing
322290001Sglebius	 * translation unit: a character string literal of the form "Mmm
323290001Sglebius	 * dd yyyy", where the names of the months are the same as those
324290001Sglebius	 * generated by the asctime function, and the first character of
325290001Sglebius	 * dd is a space character if the value is less than 10. If the
326290001Sglebius	 * date of translation is not available, an
327290001Sglebius	 * implementation-defined valid date shall be supplied.
328290001Sglebius	 *
329290001Sglebius	 * __TIME__ The time of translation of the preprocessing
330290001Sglebius	 * translation unit: a character string literal of the form
331290001Sglebius	 * "hh:mm:ss" as in the time generated by the asctime
332290001Sglebius	 * function. If the time of translation is not available, an
333290001Sglebius	 * implementation-defined valid time shall be supplied.
334290001Sglebius	 *
335290001Sglebius	 * Note that MSVC declares DATE and TIME to be in the local time
336290001Sglebius	 * zone, while neither the C standard nor the GCC docs make any
337290001Sglebius	 * statement about this. As a result, we may be +/-12hrs off
338290001Sglebius	 * UTC.  But for practical purposes, this should not be a
339290001Sglebius	 * problem.
340290001Sglebius	 *
341290001Sglebius	 */
342290001Sglebius#   ifdef MKREPRO_DATE
343290001Sglebius	static const char build[] = MKREPRO_TIME "/" MKREPRO_DATE;
344290001Sglebius#   else
345290001Sglebius	static const char build[] = __TIME__ "/" __DATE__;
346290001Sglebius#   endif
347290001Sglebius	static const char mlist[] = "JanFebMarAprMayJunJulAugSepOctNovDec";
348290001Sglebius
349290001Sglebius	char		  monstr[4];
350290001Sglebius	const char *	  cp;
351290001Sglebius	unsigned short	  hour, minute, second, day, year;
352290001Sglebius 	/* Note: The above quantities are used for sscanf 'hu' format,
353290001Sglebius	 * so using 'uint16_t' is contra-indicated!
354290001Sglebius	 */
355290001Sglebius
356290001Sglebius#   ifdef DEBUG
357290001Sglebius	static int        ignore  = 0;
358290001Sglebius#   endif
359290001Sglebius
360290001Sglebius	ZERO(*jd);
361290001Sglebius	jd->year     = 1970;
362290001Sglebius	jd->month    = 1;
363290001Sglebius	jd->monthday = 1;
364290001Sglebius
365290001Sglebius#   ifdef DEBUG
366290001Sglebius	/* check environment if build date should be ignored */
367290001Sglebius	if (0 == ignore) {
368290001Sglebius	    const char * envstr;
369290001Sglebius	    envstr = getenv("NTPD_IGNORE_BUILD_DATE");
370290001Sglebius	    ignore = 1 + (envstr && (!*envstr || !strcasecmp(envstr, "yes")));
371290001Sglebius	}
372290001Sglebius	if (ignore > 1)
373290001Sglebius	    return FALSE;
374290001Sglebius#   endif
375290001Sglebius
376290001Sglebius	if (6 == sscanf(build, "%hu:%hu:%hu/%3s %hu %hu",
377290001Sglebius			&hour, &minute, &second, monstr, &day, &year)) {
378290001Sglebius		cp = strstr(mlist, monstr);
379290001Sglebius		if (NULL != cp) {
380290001Sglebius			jd->year     = year;
381290001Sglebius			jd->month    = (uint8_t)((cp - mlist) / 3 + 1);
382290001Sglebius			jd->monthday = (uint8_t)day;
383290001Sglebius			jd->hour     = (uint8_t)hour;
384290001Sglebius			jd->minute   = (uint8_t)minute;
385290001Sglebius			jd->second   = (uint8_t)second;
386290001Sglebius
387290001Sglebius			return TRUE;
388290001Sglebius		}
389290001Sglebius	}
390290001Sglebius
391290001Sglebius	return FALSE;
392290001Sglebius}
393290001Sglebius
394290001Sglebius
395290001Sglebius/*
396290001Sglebius *---------------------------------------------------------------------
397290001Sglebius * basic calendar stuff
398290001Sglebius * --------------------------------------------------------------------
399290001Sglebius */
400290001Sglebius
401290001Sglebius/* month table for a year starting with March,1st */
402290001Sglebiusstatic const uint16_t shift_month_table[13] = {
403290001Sglebius	0, 31, 61, 92, 122, 153, 184, 214, 245, 275, 306, 337, 366
404290001Sglebius};
405290001Sglebius
406290001Sglebius/* month tables for years starting with January,1st; regular & leap */
407290001Sglebiusstatic const uint16_t real_month_table[2][13] = {
408290001Sglebius	/* -*- table for regular years -*- */
409290001Sglebius	{ 0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334, 365 },
410290001Sglebius	/* -*- table for leap years -*- */
411290001Sglebius	{ 0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335, 366 }
412290001Sglebius};
413290001Sglebius
414290001Sglebius/*
415290001Sglebius * Some notes on the terminology:
416290001Sglebius *
417290001Sglebius * We use the proleptic Gregorian calendar, which is the Gregorian
418290001Sglebius * calendar extended in both directions ad infinitum. This totally
419290001Sglebius * disregards the fact that this calendar was invented in 1582, and
420290001Sglebius * was adopted at various dates over the world; sometimes even after
421290001Sglebius * the start of the NTP epoch.
422290001Sglebius *
423290001Sglebius * Normally date parts are given as current cycles, while time parts
424290001Sglebius * are given as elapsed cycles:
425290001Sglebius *
426290001Sglebius * 1970-01-01/03:04:05 means 'IN the 1970st. year, IN the first month,
427290001Sglebius * ON the first day, with 3hrs, 4minutes and 5 seconds elapsed.
428290001Sglebius *
429290001Sglebius * The basic calculations for this calendar implementation deal with
430290001Sglebius * ELAPSED date units, which is the number of full years, full months
431290001Sglebius * and full days before a date: 1970-01-01 would be (1969, 0, 0) in
432290001Sglebius * that notation.
433290001Sglebius *
434290001Sglebius * To ease the numeric computations, month and day values outside the
435290001Sglebius * normal range are acceptable: 2001-03-00 will be treated as the day
436290001Sglebius * before 2001-03-01, 2000-13-32 will give the same result as
437290001Sglebius * 2001-02-01 and so on.
438290001Sglebius *
439290001Sglebius * 'rd' or 'RD' is used as an abbreviation for the latin 'rata die'
440290001Sglebius * (day number).  This is the number of days elapsed since 0000-12-31
441290001Sglebius * in the proleptic Gregorian calendar. The begin of the Christian Era
442290001Sglebius * (0001-01-01) is RD(1).
443290001Sglebius */
444290001Sglebius
445290001Sglebius/*
446290001Sglebius * ==================================================================
447290001Sglebius *
448290001Sglebius * General algorithmic stuff
449290001Sglebius *
450290001Sglebius * ==================================================================
451290001Sglebius */
452290001Sglebius
453290001Sglebius/*
454290001Sglebius *---------------------------------------------------------------------
455290001Sglebius * Do a periodic extension of 'value' around 'pivot' with a period of
456290001Sglebius * 'cycle'.
457290001Sglebius *
458290001Sglebius * The result 'res' is a number that holds to the following properties:
459290001Sglebius *
460290001Sglebius *   1)	 res MOD cycle == value MOD cycle
461290001Sglebius *   2)	 pivot <= res < pivot + cycle
462290001Sglebius *	 (replace </<= with >/>= for negative cycles)
463290001Sglebius *
464290001Sglebius * where 'MOD' denotes the modulo operator for FLOOR DIVISION, which
465290001Sglebius * is not the same as the '%' operator in C: C requires division to be
466290001Sglebius * a truncated division, where remainder and dividend have the same
467290001Sglebius * sign if the remainder is not zero, whereas floor division requires
468290001Sglebius * divider and modulus to have the same sign for a non-zero modulus.
469290001Sglebius *
470290001Sglebius * This function has some useful applications:
471290001Sglebius *
472290001Sglebius * + let Y be a calendar year and V a truncated 2-digit year: then
473290001Sglebius *	periodic_extend(Y-50, V, 100)
474290001Sglebius *   is the closest expansion of the truncated year with respect to
475290001Sglebius *   the full year, that is a 4-digit year with a difference of less
476290001Sglebius *   than 50 years to the year Y. ("century unfolding")
477290001Sglebius *
478290001Sglebius * + let T be a UN*X time stamp and V be seconds-of-day: then
479290001Sglebius *	perodic_extend(T-43200, V, 86400)
480290001Sglebius *   is a time stamp that has the same seconds-of-day as the input
481290001Sglebius *   value, with an absolute difference to T of <= 12hrs.  ("day
482290001Sglebius *   unfolding")
483290001Sglebius *
484290001Sglebius * + Wherever you have a truncated periodic value and a non-truncated
485290001Sglebius *   base value and you want to match them somehow...
486290001Sglebius *
487290001Sglebius * Basically, the function delivers 'pivot + (value - pivot) % cycle',
488290001Sglebius * but the implementation takes some pains to avoid internal signed
489290001Sglebius * integer overflows in the '(value - pivot) % cycle' part and adheres
490290001Sglebius * to the floor division convention.
491290001Sglebius *
492290001Sglebius * If 64bit scalars where available on all intended platforms, writing a
493290001Sglebius * version that uses 64 bit ops would be easy; writing a general
494290001Sglebius * division routine for 64bit ops on a platform that can only do
495290001Sglebius * 32/16bit divisions and is still performant is a bit more
496290001Sglebius * difficult. Since most usecases can be coded in a way that does only
497290001Sglebius * require the 32-bit version a 64bit version is NOT provided here.
498290001Sglebius * ---------------------------------------------------------------------
499290001Sglebius */
500290001Sglebiusint32_t
501290001Sglebiusntpcal_periodic_extend(
502290001Sglebius	int32_t pivot,
503290001Sglebius	int32_t value,
504290001Sglebius	int32_t cycle
505290001Sglebius	)
506290001Sglebius{
507290001Sglebius	uint32_t diff;
508290001Sglebius	char	 cpl = 0; /* modulo complement flag */
509290001Sglebius	char	 neg = 0; /* sign change flag	    */
510290001Sglebius
511290001Sglebius	/* make the cycle positive and adjust the flags */
512290001Sglebius	if (cycle < 0) {
513290001Sglebius		cycle = - cycle;
514290001Sglebius		neg ^= 1;
515290001Sglebius		cpl ^= 1;
516290001Sglebius	}
517290001Sglebius	/* guard against div by zero or one */
518290001Sglebius	if (cycle > 1) {
519290001Sglebius		/*
520290001Sglebius		 * Get absolute difference as unsigned quantity and
521290001Sglebius		 * the complement flag. This is done by always
522290001Sglebius		 * subtracting the smaller value from the bigger
523290001Sglebius		 * one.
524290001Sglebius		 */
525290001Sglebius		if (value >= pivot) {
526290001Sglebius			diff = int32_to_uint32_2cpl(value)
527290001Sglebius			     - int32_to_uint32_2cpl(pivot);
528290001Sglebius		} else {
529290001Sglebius			diff = int32_to_uint32_2cpl(pivot)
530290001Sglebius			     - int32_to_uint32_2cpl(value);
531290001Sglebius			cpl ^= 1;
532290001Sglebius		}
533290001Sglebius		diff %= (uint32_t)cycle;
534290001Sglebius		if (diff) {
535290001Sglebius			if (cpl)
536290001Sglebius				diff = (uint32_t)cycle - diff;
537290001Sglebius			if (neg)
538290001Sglebius				diff = ~diff + 1;
539290001Sglebius			pivot += uint32_2cpl_to_int32(diff);
540290001Sglebius		}
541290001Sglebius	}
542290001Sglebius	return pivot;
543290001Sglebius}
544290001Sglebius
545290001Sglebius/*
546290001Sglebius *-------------------------------------------------------------------
547290001Sglebius * Convert a timestamp in NTP scale to a 64bit seconds value in the UN*X
548290001Sglebius * scale with proper epoch unfolding around a given pivot or the current
549290001Sglebius * system time. This function happily accepts negative pivot values as
550290001Sglebius * timestamps befor 1970-01-01, so be aware of possible trouble on
551290001Sglebius * platforms with 32bit 'time_t'!
552290001Sglebius *
553290001Sglebius * This is also a periodic extension, but since the cycle is 2^32 and
554290001Sglebius * the shift is 2^31, we can do some *very* fast math without explicit
555290001Sglebius * divisions.
556290001Sglebius *-------------------------------------------------------------------
557290001Sglebius */
558290001Sglebiusvint64
559290001Sglebiusntpcal_ntp_to_time(
560290001Sglebius	uint32_t	ntp,
561290001Sglebius	const time_t *	pivot
562290001Sglebius	)
563290001Sglebius{
564290001Sglebius	vint64 res;
565290001Sglebius
566290001Sglebius#   if defined(HAVE_INT64)
567290001Sglebius
568290001Sglebius	res.q_s = (pivot != NULL)
569290001Sglebius		      ? *pivot
570290001Sglebius		      : now();
571290001Sglebius	res.Q_s -= 0x80000000;		/* unshift of half range */
572290001Sglebius	ntp	-= (uint32_t)JAN_1970;	/* warp into UN*X domain */
573290001Sglebius	ntp	-= res.D_s.lo;		/* cycle difference	 */
574290001Sglebius	res.Q_s += (uint64_t)ntp;	/* get expanded time	 */
575290001Sglebius
576290001Sglebius#   else /* no 64bit scalars */
577290001Sglebius
578290001Sglebius	time_t tmp;
579290001Sglebius
580290001Sglebius	tmp = (pivot != NULL)
581290001Sglebius		  ? *pivot
582290001Sglebius		  : now();
583290001Sglebius	res = time_to_vint64(&tmp);
584290001Sglebius	M_SUB(res.D_s.hi, res.D_s.lo, 0, 0x80000000);
585290001Sglebius	ntp -= (uint32_t)JAN_1970;	/* warp into UN*X domain */
586290001Sglebius	ntp -= res.D_s.lo;		/* cycle difference	 */
587290001Sglebius	M_ADD(res.D_s.hi, res.D_s.lo, 0, ntp);
588290001Sglebius
589290001Sglebius#   endif /* no 64bit scalars */
590290001Sglebius
591290001Sglebius	return res;
592290001Sglebius}
593290001Sglebius
594290001Sglebius/*
595290001Sglebius *-------------------------------------------------------------------
596290001Sglebius * Convert a timestamp in NTP scale to a 64bit seconds value in the NTP
597290001Sglebius * scale with proper epoch unfolding around a given pivot or the current
598290001Sglebius * system time.
599290001Sglebius *
600290001Sglebius * Note: The pivot must be given in the UN*X time domain!
601290001Sglebius *
602290001Sglebius * This is also a periodic extension, but since the cycle is 2^32 and
603290001Sglebius * the shift is 2^31, we can do some *very* fast math without explicit
604290001Sglebius * divisions.
605290001Sglebius *-------------------------------------------------------------------
606290001Sglebius */
607290001Sglebiusvint64
608290001Sglebiusntpcal_ntp_to_ntp(
609290001Sglebius	uint32_t      ntp,
610290001Sglebius	const time_t *pivot
611290001Sglebius	)
612290001Sglebius{
613290001Sglebius	vint64 res;
614290001Sglebius
615290001Sglebius#   if defined(HAVE_INT64)
616290001Sglebius
617290001Sglebius	res.q_s = (pivot)
618290001Sglebius		      ? *pivot
619290001Sglebius		      : now();
620290001Sglebius	res.Q_s -= 0x80000000;		/* unshift of half range */
621290001Sglebius	res.Q_s += (uint32_t)JAN_1970;	/* warp into NTP domain	 */
622290001Sglebius	ntp	-= res.D_s.lo;		/* cycle difference	 */
623290001Sglebius	res.Q_s += (uint64_t)ntp;	/* get expanded time	 */
624290001Sglebius
625290001Sglebius#   else /* no 64bit scalars */
626290001Sglebius
627290001Sglebius	time_t tmp;
628290001Sglebius
629290001Sglebius	tmp = (pivot)
630290001Sglebius		  ? *pivot
631290001Sglebius		  : now();
632290001Sglebius	res = time_to_vint64(&tmp);
633290001Sglebius	M_SUB(res.D_s.hi, res.D_s.lo, 0, 0x80000000u);
634290001Sglebius	M_ADD(res.D_s.hi, res.D_s.lo, 0, (uint32_t)JAN_1970);/*into NTP */
635290001Sglebius	ntp -= res.D_s.lo;		/* cycle difference	 */
636290001Sglebius	M_ADD(res.D_s.hi, res.D_s.lo, 0, ntp);
637290001Sglebius
638290001Sglebius#   endif /* no 64bit scalars */
639290001Sglebius
640290001Sglebius	return res;
641290001Sglebius}
642290001Sglebius
643290001Sglebius
644290001Sglebius/*
645290001Sglebius * ==================================================================
646290001Sglebius *
647290001Sglebius * Splitting values to composite entities
648290001Sglebius *
649290001Sglebius * ==================================================================
650290001Sglebius */
651290001Sglebius
652290001Sglebius/*
653290001Sglebius *-------------------------------------------------------------------
654290001Sglebius * Split a 64bit seconds value into elapsed days in 'res.hi' and
655290001Sglebius * elapsed seconds since midnight in 'res.lo' using explicit floor
656290001Sglebius * division. This function happily accepts negative time values as
657290001Sglebius * timestamps before the respective epoch start.
658290001Sglebius * -------------------------------------------------------------------
659290001Sglebius */
660290001Sglebiusntpcal_split
661290001Sglebiusntpcal_daysplit(
662290001Sglebius	const vint64 *ts
663290001Sglebius	)
664290001Sglebius{
665290001Sglebius	ntpcal_split res;
666290001Sglebius	uint32_t Q;
667290001Sglebius
668290001Sglebius#   if defined(HAVE_INT64)
669290001Sglebius
670290001Sglebius	/* Manual floor division by SECSPERDAY. This uses the one's
671290001Sglebius	 * complement trick, too, but without an extra flag value: The
672290001Sglebius	 * flag would be 64bit, and that's a bit of overkill on a 32bit
673290001Sglebius	 * target that has to use a register pair for a 64bit number.
674290001Sglebius	 */
675290001Sglebius	if (ts->q_s < 0)
676290001Sglebius		Q = ~(uint32_t)(~ts->Q_s / SECSPERDAY);
677290001Sglebius	else
678290001Sglebius		Q = (uint32_t)(ts->Q_s / SECSPERDAY);
679290001Sglebius
680290001Sglebius#   else
681290001Sglebius
682290001Sglebius	uint32_t ah, al, sflag, A;
683290001Sglebius
684290001Sglebius	/* get operand into ah/al (either ts or ts' one's complement,
685290001Sglebius	 * for later floor division)
686290001Sglebius	 */
687290001Sglebius	sflag = int32_sflag(ts->d_s.hi);
688290001Sglebius	ah = sflag ^ ts->D_s.hi;
689290001Sglebius	al = sflag ^ ts->D_s.lo;
690290001Sglebius
691290001Sglebius	/* Since 86400 == 128*675 we can drop the least 7 bits and
692290001Sglebius	 * divide by 675 instead of 86400. Then the maximum remainder
693290001Sglebius	 * after each devision step is 674, and we need 10 bits for
694290001Sglebius	 * that. So in the next step we can shift in 22 bits from the
695290001Sglebius	 * numerator.
696290001Sglebius	 *
697290001Sglebius	 * Therefore we load the accu with the top 13 bits (51..63) in
698290001Sglebius	 * the first shot. We don't have to remember the quotient -- it
699290001Sglebius	 * would be shifted out anyway.
700290001Sglebius	 */
701290001Sglebius	A = ah >> 19;
702290001Sglebius	if (A >= 675)
703290001Sglebius		A = (A % 675u);
704290001Sglebius
705290001Sglebius	/* Now assemble the remainder with bits 29..50 from the
706290001Sglebius	 * numerator and divide. This creates the upper ten bits of the
707290001Sglebius	 * quotient. (Well, the top 22 bits of a 44bit result. But that
708290001Sglebius	 * will be truncated to 32 bits anyway.)
709290001Sglebius	 */
710290001Sglebius	A = (A << 19) | (ah & 0x0007FFFFu);
711290001Sglebius	A = (A <<  3) | (al >> 29);
712290001Sglebius	Q = A / 675u;
713290001Sglebius	A = A % 675u;
714290001Sglebius
715290001Sglebius	/* Now assemble the remainder with bits 7..28 from the numerator
716290001Sglebius	 * and do a final division step.
717290001Sglebius	 */
718290001Sglebius	A = (A << 22) | ((al >> 7) & 0x003FFFFFu);
719290001Sglebius	Q = (Q << 22) | (A / 675u);
720290001Sglebius
721290001Sglebius	/* The last 7 bits get simply dropped, as they have no affect on
722290001Sglebius	 * the quotient when dividing by 86400.
723290001Sglebius	 */
724290001Sglebius
725290001Sglebius	/* apply sign correction and calculate the true floor
726290001Sglebius	 * remainder.
727290001Sglebius	 */
728290001Sglebius	Q ^= sflag;
729290001Sglebius
730290001Sglebius#   endif
731290001Sglebius
732290001Sglebius	res.hi = uint32_2cpl_to_int32(Q);
733290001Sglebius	res.lo = ts->D_s.lo - Q * SECSPERDAY;
734290001Sglebius
735290001Sglebius	return res;
736290001Sglebius}
737290001Sglebius
738290001Sglebius/*
739290001Sglebius *-------------------------------------------------------------------
740290001Sglebius * Split a 32bit seconds value into h/m/s and excessive days.  This
741290001Sglebius * function happily accepts negative time values as timestamps before
742290001Sglebius * midnight.
743290001Sglebius * -------------------------------------------------------------------
744290001Sglebius */
745290001Sglebiusstatic int32_t
746290001Sglebiuspriv_timesplit(
747290001Sglebius	int32_t split[3],
748290001Sglebius	int32_t ts
749290001Sglebius	)
750290001Sglebius{
751290001Sglebius	/* Do 3 chained floor divisions by positive constants, using the
752290001Sglebius	 * one's complement trick and factoring out the intermediate XOR
753290001Sglebius	 * ops to reduce the number of operations.
754290001Sglebius	 */
755290001Sglebius	uint32_t us, um, uh, ud, sflag;
756290001Sglebius
757290001Sglebius	sflag = int32_sflag(ts);
758290001Sglebius	us    = int32_to_uint32_2cpl(ts);
759290001Sglebius
760290001Sglebius	um = (sflag ^ us) / SECSPERMIN;
761290001Sglebius	uh = um / MINSPERHR;
762290001Sglebius	ud = uh / HRSPERDAY;
763290001Sglebius
764290001Sglebius	um ^= sflag;
765290001Sglebius	uh ^= sflag;
766290001Sglebius	ud ^= sflag;
767290001Sglebius
768290001Sglebius	split[0] = (int32_t)(uh - ud * HRSPERDAY );
769290001Sglebius	split[1] = (int32_t)(um - uh * MINSPERHR );
770290001Sglebius	split[2] = (int32_t)(us - um * SECSPERMIN);
771290001Sglebius
772290001Sglebius	return uint32_2cpl_to_int32(ud);
773290001Sglebius}
774290001Sglebius
775290001Sglebius/*
776290001Sglebius * ---------------------------------------------------------------------
777290001Sglebius * Given the number of elapsed days in the calendar era, split this
778290001Sglebius * number into the number of elapsed years in 'res.hi' and the number
779290001Sglebius * of elapsed days of that year in 'res.lo'.
780290001Sglebius *
781290001Sglebius * if 'isleapyear' is not NULL, it will receive an integer that is 0 for
782290001Sglebius * regular years and a non-zero value for leap years.
783290001Sglebius *---------------------------------------------------------------------
784290001Sglebius */
785290001Sglebiusntpcal_split
786290001Sglebiusntpcal_split_eradays(
787290001Sglebius	int32_t days,
788290001Sglebius	int  *isleapyear
789290001Sglebius	)
790290001Sglebius{
791290001Sglebius	/* Use the fast cyclesplit algorithm here, to calculate the
792290001Sglebius	 * centuries and years in a century with one division each. This
793290001Sglebius	 * reduces the number of division operations to two, but is
794290001Sglebius	 * susceptible to internal range overflow. We make sure the
795290001Sglebius	 * input operands are in the safe range; this still gives us
796290001Sglebius	 * approx +/-2.9 million years.
797290001Sglebius	 */
798290001Sglebius	ntpcal_split res;
799290001Sglebius	int32_t	 n100, n001; /* calendar year cycles */
800290001Sglebius	uint32_t uday, Q, sflag;
801290001Sglebius
802290001Sglebius	/* split off centuries first */
803290001Sglebius	sflag = int32_sflag(days);
804290001Sglebius	uday  = uint32_saturate(int32_to_uint32_2cpl(days), sflag);
805290001Sglebius	uday  = (4u * uday) | 3u;
806290001Sglebius	Q    = sflag ^ ((sflag ^ uday) / GREGORIAN_CYCLE_DAYS);
807290001Sglebius	uday = uday - Q * GREGORIAN_CYCLE_DAYS;
808290001Sglebius	n100 = uint32_2cpl_to_int32(Q);
809290001Sglebius
810290001Sglebius	/* Split off years in century -- days >= 0 here, and we're far
811290001Sglebius	 * away from integer overflow trouble now. */
812290001Sglebius	uday |= 3;
813290001Sglebius	n001 = uday / GREGORIAN_NORMAL_LEAP_CYCLE_DAYS;
814290001Sglebius	uday = uday % GREGORIAN_NORMAL_LEAP_CYCLE_DAYS;
815290001Sglebius
816290001Sglebius	/* Assemble the year and day in year */
817290001Sglebius	res.hi = n100 * 100 + n001;
818290001Sglebius	res.lo = uday / 4u;
819290001Sglebius
820290001Sglebius	/* Eventually set the leap year flag. Note: 0 <= n001 <= 99 and
821290001Sglebius	 * Q is still the two's complement representation of the
822290001Sglebius	 * centuries: The modulo 4 ops can be done with masking here.
823290001Sglebius	 * We also shift the year and the century by one, so the tests
824290001Sglebius	 * can be done against zero instead of 3.
825290001Sglebius	 */
826290001Sglebius	if (isleapyear)
827290001Sglebius		*isleapyear = !((n001+1) & 3)
828290001Sglebius		    && ((n001 != 99) || !((Q+1) & 3));
829290001Sglebius
830290001Sglebius	return res;
831290001Sglebius}
832290001Sglebius
833290001Sglebius/*
834290001Sglebius *---------------------------------------------------------------------
835290001Sglebius * Given a number of elapsed days in a year and a leap year indicator,
836290001Sglebius * split the number of elapsed days into the number of elapsed months in
837290001Sglebius * 'res.hi' and the number of elapsed days of that month in 'res.lo'.
838290001Sglebius *
839290001Sglebius * This function will fail and return {-1,-1} if the number of elapsed
840290001Sglebius * days is not in the valid range!
841290001Sglebius *---------------------------------------------------------------------
842290001Sglebius */
843290001Sglebiusntpcal_split
844290001Sglebiusntpcal_split_yeardays(
845290001Sglebius	int32_t eyd,
846290001Sglebius	int     isleapyear
847290001Sglebius	)
848290001Sglebius{
849290001Sglebius	ntpcal_split    res;
850290001Sglebius	const uint16_t *lt;	/* month length table	*/
851290001Sglebius
852290001Sglebius	/* check leap year flag and select proper table */
853290001Sglebius	lt = real_month_table[(isleapyear != 0)];
854290001Sglebius	if (0 <= eyd && eyd < lt[12]) {
855290001Sglebius		/* get zero-based month by approximation & correction step */
856290001Sglebius		res.hi = eyd >> 5;	   /* approx month; might be 1 too low */
857290001Sglebius		if (lt[res.hi + 1] <= eyd) /* fixup approximative month value  */
858290001Sglebius			res.hi += 1;
859290001Sglebius		res.lo = eyd - lt[res.hi];
860290001Sglebius	} else {
861290001Sglebius		res.lo = res.hi = -1;
862290001Sglebius	}
863290001Sglebius
864290001Sglebius	return res;
865290001Sglebius}
866290001Sglebius
867290001Sglebius/*
868290001Sglebius *---------------------------------------------------------------------
869290001Sglebius * Convert a RD into the date part of a 'struct calendar'.
870290001Sglebius *---------------------------------------------------------------------
871290001Sglebius */
872290001Sglebiusint
873290001Sglebiusntpcal_rd_to_date(
874290001Sglebius	struct calendar *jd,
875290001Sglebius	int32_t		 rd
876290001Sglebius	)
877290001Sglebius{
878290001Sglebius	ntpcal_split split;
879290001Sglebius	int	     leapy;
880290001Sglebius	u_int	     ymask;
881290001Sglebius
882290001Sglebius	/* Get day-of-week first. Since rd is signed, the remainder can
883290001Sglebius	 * be in the range [-6..+6], but the assignment to an unsigned
884290001Sglebius	 * variable maps the negative values to positive values >=7.
885290001Sglebius	 * This makes the sign correction look strange, but adding 7
886290001Sglebius	 * causes the needed wrap-around into the desired value range of
887290001Sglebius	 * zero to six, both inclusive.
888290001Sglebius	 */
889290001Sglebius	jd->weekday = rd % DAYSPERWEEK;
890290001Sglebius	if (jd->weekday >= DAYSPERWEEK)	/* weekday is unsigned! */
891290001Sglebius		jd->weekday += DAYSPERWEEK;
892290001Sglebius
893290001Sglebius	split = ntpcal_split_eradays(rd - 1, &leapy);
894290001Sglebius	/* Get year and day-of-year, with overflow check. If any of the
895290001Sglebius	 * upper 16 bits is set after shifting to unity-based years, we
896290001Sglebius	 * will have an overflow when converting to an unsigned 16bit
897290001Sglebius	 * year. Shifting to the right is OK here, since it does not
898290001Sglebius	 * matter if the shift is logic or arithmetic.
899290001Sglebius	 */
900290001Sglebius	split.hi += 1;
901290001Sglebius	ymask = 0u - ((split.hi >> 16) == 0);
902290001Sglebius	jd->year = (uint16_t)(split.hi & ymask);
903290001Sglebius	jd->yearday = (uint16_t)split.lo + 1;
904290001Sglebius
905290001Sglebius	/* convert to month and mday */
906290001Sglebius	split = ntpcal_split_yeardays(split.lo, leapy);
907290001Sglebius	jd->month    = (uint8_t)split.hi + 1;
908290001Sglebius	jd->monthday = (uint8_t)split.lo + 1;
909290001Sglebius
910290001Sglebius	return ymask ? leapy : -1;
911290001Sglebius}
912290001Sglebius
913290001Sglebius/*
914290001Sglebius *---------------------------------------------------------------------
915290001Sglebius * Convert a RD into the date part of a 'struct tm'.
916290001Sglebius *---------------------------------------------------------------------
917290001Sglebius */
918290001Sglebiusint
919290001Sglebiusntpcal_rd_to_tm(
920290001Sglebius	struct tm  *utm,
921290001Sglebius	int32_t	    rd
922290001Sglebius	)
923290001Sglebius{
924290001Sglebius	ntpcal_split split;
925290001Sglebius	int	     leapy;
926290001Sglebius
927290001Sglebius	/* get day-of-week first */
928290001Sglebius	utm->tm_wday = rd % DAYSPERWEEK;
929290001Sglebius	if (utm->tm_wday < 0)
930290001Sglebius		utm->tm_wday += DAYSPERWEEK;
931290001Sglebius
932290001Sglebius	/* get year and day-of-year */
933290001Sglebius	split = ntpcal_split_eradays(rd - 1, &leapy);
934290001Sglebius	utm->tm_year = split.hi - 1899;
935290001Sglebius	utm->tm_yday = split.lo;	/* 0-based */
936290001Sglebius
937290001Sglebius	/* convert to month and mday */
938290001Sglebius	split = ntpcal_split_yeardays(split.lo, leapy);
939290001Sglebius	utm->tm_mon  = split.hi;	/* 0-based */
940290001Sglebius	utm->tm_mday = split.lo + 1;	/* 1-based */
941290001Sglebius
942290001Sglebius	return leapy;
943290001Sglebius}
944290001Sglebius
945290001Sglebius/*
946290001Sglebius *---------------------------------------------------------------------
947290001Sglebius * Take a value of seconds since midnight and split it into hhmmss in a
948290001Sglebius * 'struct calendar'.
949290001Sglebius *---------------------------------------------------------------------
950290001Sglebius */
951290001Sglebiusint32_t
952290001Sglebiusntpcal_daysec_to_date(
953290001Sglebius	struct calendar *jd,
954290001Sglebius	int32_t		sec
955290001Sglebius	)
956290001Sglebius{
957290001Sglebius	int32_t days;
958290001Sglebius	int   ts[3];
959290001Sglebius
960290001Sglebius	days = priv_timesplit(ts, sec);
961290001Sglebius	jd->hour   = (uint8_t)ts[0];
962290001Sglebius	jd->minute = (uint8_t)ts[1];
963290001Sglebius	jd->second = (uint8_t)ts[2];
964290001Sglebius
965290001Sglebius	return days;
966290001Sglebius}
967290001Sglebius
968290001Sglebius/*
969290001Sglebius *---------------------------------------------------------------------
970290001Sglebius * Take a value of seconds since midnight and split it into hhmmss in a
971290001Sglebius * 'struct tm'.
972290001Sglebius *---------------------------------------------------------------------
973290001Sglebius */
974290001Sglebiusint32_t
975290001Sglebiusntpcal_daysec_to_tm(
976290001Sglebius	struct tm *utm,
977290001Sglebius	int32_t	   sec
978290001Sglebius	)
979290001Sglebius{
980290001Sglebius	int32_t days;
981290001Sglebius	int32_t ts[3];
982290001Sglebius
983290001Sglebius	days = priv_timesplit(ts, sec);
984290001Sglebius	utm->tm_hour = ts[0];
985290001Sglebius	utm->tm_min  = ts[1];
986290001Sglebius	utm->tm_sec  = ts[2];
987290001Sglebius
988290001Sglebius	return days;
989290001Sglebius}
990290001Sglebius
991290001Sglebius/*
992290001Sglebius *---------------------------------------------------------------------
993290001Sglebius * take a split representation for day/second-of-day and day offset
994290001Sglebius * and convert it to a 'struct calendar'. The seconds will be normalised
995290001Sglebius * into the range of a day, and the day will be adjusted accordingly.
996290001Sglebius *
997290001Sglebius * returns >0 if the result is in a leap year, 0 if in a regular
998290001Sglebius * year and <0 if the result did not fit into the calendar struct.
999290001Sglebius *---------------------------------------------------------------------
1000290001Sglebius */
1001290001Sglebiusint
1002290001Sglebiusntpcal_daysplit_to_date(
1003290001Sglebius	struct calendar	   *jd,
1004290001Sglebius	const ntpcal_split *ds,
1005290001Sglebius	int32_t		    dof
1006290001Sglebius	)
1007290001Sglebius{
1008290001Sglebius	dof += ntpcal_daysec_to_date(jd, ds->lo);
1009290001Sglebius	return ntpcal_rd_to_date(jd, ds->hi + dof);
1010290001Sglebius}
1011290001Sglebius
1012290001Sglebius/*
1013290001Sglebius *---------------------------------------------------------------------
1014290001Sglebius * take a split representation for day/second-of-day and day offset
1015290001Sglebius * and convert it to a 'struct tm'. The seconds will be normalised
1016290001Sglebius * into the range of a day, and the day will be adjusted accordingly.
1017290001Sglebius *
1018290001Sglebius * returns 1 if the result is in a leap year and zero if in a regular
1019290001Sglebius * year.
1020290001Sglebius *---------------------------------------------------------------------
1021290001Sglebius */
1022290001Sglebiusint
1023290001Sglebiusntpcal_daysplit_to_tm(
1024290001Sglebius	struct tm	   *utm,
1025290001Sglebius	const ntpcal_split *ds ,
1026290001Sglebius	int32_t		    dof
1027290001Sglebius	)
1028290001Sglebius{
1029290001Sglebius	dof += ntpcal_daysec_to_tm(utm, ds->lo);
1030290001Sglebius
1031290001Sglebius	return ntpcal_rd_to_tm(utm, ds->hi + dof);
1032290001Sglebius}
1033290001Sglebius
1034290001Sglebius/*
1035290001Sglebius *---------------------------------------------------------------------
1036290001Sglebius * Take a UN*X time and convert to a calendar structure.
1037290001Sglebius *---------------------------------------------------------------------
1038290001Sglebius */
1039290001Sglebiusint
1040290001Sglebiusntpcal_time_to_date(
1041290001Sglebius	struct calendar	*jd,
1042290001Sglebius	const vint64	*ts
1043290001Sglebius	)
1044290001Sglebius{
1045290001Sglebius	ntpcal_split ds;
1046290001Sglebius
1047290001Sglebius	ds = ntpcal_daysplit(ts);
1048290001Sglebius	ds.hi += ntpcal_daysec_to_date(jd, ds.lo);
1049290001Sglebius	ds.hi += DAY_UNIX_STARTS;
1050290001Sglebius
1051290001Sglebius	return ntpcal_rd_to_date(jd, ds.hi);
1052290001Sglebius}
1053290001Sglebius
1054290001Sglebius
1055290001Sglebius/*
1056290001Sglebius * ==================================================================
1057290001Sglebius *
1058290001Sglebius * merging composite entities
1059290001Sglebius *
1060290001Sglebius * ==================================================================
1061290001Sglebius */
1062290001Sglebius
1063290001Sglebius/*
1064290001Sglebius *---------------------------------------------------------------------
1065290001Sglebius * Merge a number of days and a number of seconds into seconds,
1066290001Sglebius * expressed in 64 bits to avoid overflow.
1067290001Sglebius *---------------------------------------------------------------------
1068290001Sglebius */
1069290001Sglebiusvint64
1070290001Sglebiusntpcal_dayjoin(
1071290001Sglebius	int32_t days,
1072290001Sglebius	int32_t secs
1073290001Sglebius	)
1074290001Sglebius{
1075290001Sglebius	vint64 res;
1076290001Sglebius
1077290001Sglebius#   if defined(HAVE_INT64)
1078290001Sglebius
1079290001Sglebius	res.q_s	 = days;
1080290001Sglebius	res.q_s *= SECSPERDAY;
1081290001Sglebius	res.q_s += secs;
1082290001Sglebius
1083290001Sglebius#   else
1084290001Sglebius
1085290001Sglebius	uint32_t p1, p2;
1086290001Sglebius	int	 isneg;
1087290001Sglebius
1088290001Sglebius	/*
1089290001Sglebius	 * res = days *86400 + secs, using manual 16/32 bit
1090290001Sglebius	 * multiplications and shifts.
1091290001Sglebius	 */
1092290001Sglebius	isneg = (days < 0);
1093290001Sglebius	if (isneg)
1094290001Sglebius		days = -days;
1095290001Sglebius
1096290001Sglebius	/* assemble days * 675 */
1097290001Sglebius	res.D_s.lo = (days & 0xFFFF) * 675u;
1098290001Sglebius	res.D_s.hi = 0;
1099290001Sglebius	p1 = (days >> 16) * 675u;
1100290001Sglebius	p2 = p1 >> 16;
1101290001Sglebius	p1 = p1 << 16;
1102290001Sglebius	M_ADD(res.D_s.hi, res.D_s.lo, p2, p1);
1103290001Sglebius
1104290001Sglebius	/* mul by 128, using shift */
1105290001Sglebius	res.D_s.hi = (res.D_s.hi << 7) | (res.D_s.lo >> 25);
1106290001Sglebius	res.D_s.lo = (res.D_s.lo << 7);
1107290001Sglebius
1108290001Sglebius	/* fix sign */
1109290001Sglebius	if (isneg)
1110290001Sglebius		M_NEG(res.D_s.hi, res.D_s.lo);
1111290001Sglebius
1112290001Sglebius	/* properly add seconds */
1113290001Sglebius	p2 = 0;
1114290001Sglebius	if (secs < 0) {
1115290001Sglebius		p1 = (uint32_t)-secs;
1116290001Sglebius		M_NEG(p2, p1);
1117290001Sglebius	} else {
1118290001Sglebius		p1 = (uint32_t)secs;
1119290001Sglebius	}
1120290001Sglebius	M_ADD(res.D_s.hi, res.D_s.lo, p2, p1);
1121290001Sglebius
1122290001Sglebius#   endif
1123290001Sglebius
1124290001Sglebius	return res;
1125290001Sglebius}
1126290001Sglebius
1127290001Sglebius/*
1128290001Sglebius *---------------------------------------------------------------------
1129290001Sglebius * get leap years since epoch in elapsed years
1130290001Sglebius *---------------------------------------------------------------------
1131290001Sglebius */
1132290001Sglebiusint32_t
1133290001Sglebiusntpcal_leapyears_in_years(
1134290001Sglebius	int32_t years
1135290001Sglebius	)
1136290001Sglebius{
1137290001Sglebius	/* We use the in-out-in algorithm here, using the one's
1138290001Sglebius	 * complement division trick for negative numbers. The chained
1139290001Sglebius	 * division sequence by 4/25/4 gives the compiler the chance to
1140290001Sglebius	 * get away with only one true division and doing shifts otherwise.
1141290001Sglebius	 */
1142290001Sglebius
1143290001Sglebius	uint32_t sflag, sum, uyear;
1144290001Sglebius
1145290001Sglebius	sflag = int32_sflag(years);
1146290001Sglebius	uyear = int32_to_uint32_2cpl(years);
1147290001Sglebius	uyear ^= sflag;
1148290001Sglebius
1149290001Sglebius	sum  = (uyear /=  4u);	/*   4yr rule --> IN  */
1150290001Sglebius	sum -= (uyear /= 25u);	/* 100yr rule --> OUT */
1151290001Sglebius	sum += (uyear /=  4u);	/* 400yr rule --> IN  */
1152290001Sglebius
1153290001Sglebius	/* Thanks to the alternation of IN/OUT/IN we can do the sum
1154290001Sglebius	 * directly and have a single one's complement operation
1155290001Sglebius	 * here. (Only if the years are negative, of course.) Otherwise
1156290001Sglebius	 * the one's complement would have to be done when
1157290001Sglebius	 * adding/subtracting the terms.
1158290001Sglebius	 */
1159290001Sglebius	return uint32_2cpl_to_int32(sflag ^ sum);
1160290001Sglebius}
1161290001Sglebius
1162290001Sglebius/*
1163290001Sglebius *---------------------------------------------------------------------
1164290001Sglebius * Convert elapsed years in Era into elapsed days in Era.
1165290001Sglebius *---------------------------------------------------------------------
1166290001Sglebius */
1167290001Sglebiusint32_t
1168290001Sglebiusntpcal_days_in_years(
1169290001Sglebius	int32_t years
1170290001Sglebius	)
1171290001Sglebius{
1172290001Sglebius	return years * DAYSPERYEAR + ntpcal_leapyears_in_years(years);
1173290001Sglebius}
1174290001Sglebius
1175290001Sglebius/*
1176290001Sglebius *---------------------------------------------------------------------
1177290001Sglebius * Convert a number of elapsed month in a year into elapsed days in year.
1178290001Sglebius *
1179290001Sglebius * The month will be normalized, and 'res.hi' will contain the
1180290001Sglebius * excessive years that must be considered when converting the years,
1181290001Sglebius * while 'res.lo' will contain the number of elapsed days since start
1182290001Sglebius * of the year.
1183290001Sglebius *
1184290001Sglebius * This code uses the shifted-month-approach to convert month to days,
1185290001Sglebius * because then there is no need to have explicit leap year
1186290001Sglebius * information.	 The slight disadvantage is that for most month values
1187290001Sglebius * the result is a negative value, and the year excess is one; the
1188290001Sglebius * conversion is then simply based on the start of the following year.
1189290001Sglebius *---------------------------------------------------------------------
1190290001Sglebius */
1191290001Sglebiusntpcal_split
1192290001Sglebiusntpcal_days_in_months(
1193290001Sglebius	int32_t m
1194290001Sglebius	)
1195290001Sglebius{
1196290001Sglebius	ntpcal_split res;
1197290001Sglebius
1198290001Sglebius	/* Add ten months and correct if needed. (It likely is...) */
1199290001Sglebius	res.lo  = m + 10;
1200290001Sglebius	res.hi  = (res.lo >= 12);
1201290001Sglebius	if (res.hi)
1202290001Sglebius		res.lo -= 12;
1203290001Sglebius
1204290001Sglebius	/* if still out of range, normalise by floor division ... */
1205290001Sglebius	if (res.lo < 0 || res.lo >= 12) {
1206290001Sglebius		uint32_t mu, Q, sflag;
1207290001Sglebius		sflag = int32_sflag(res.lo);
1208290001Sglebius		mu    = int32_to_uint32_2cpl(res.lo);
1209290001Sglebius		Q     = sflag ^ ((sflag ^ mu) / 12u);
1210290001Sglebius		res.hi += uint32_2cpl_to_int32(Q);
1211290001Sglebius		res.lo  = mu - Q * 12u;
1212290001Sglebius	}
1213290001Sglebius
1214290001Sglebius	/* get cummulated days in year with unshift */
1215290001Sglebius	res.lo = shift_month_table[res.lo] - 306;
1216290001Sglebius
1217290001Sglebius	return res;
1218290001Sglebius}
1219290001Sglebius
1220290001Sglebius/*
1221290001Sglebius *---------------------------------------------------------------------
1222290001Sglebius * Convert ELAPSED years/months/days of gregorian calendar to elapsed
1223290001Sglebius * days in Gregorian epoch.
1224290001Sglebius *
1225290001Sglebius * If you want to convert years and days-of-year, just give a month of
1226290001Sglebius * zero.
1227290001Sglebius *---------------------------------------------------------------------
1228290001Sglebius */
1229290001Sglebiusint32_t
1230290001Sglebiusntpcal_edate_to_eradays(
1231290001Sglebius	int32_t years,
1232290001Sglebius	int32_t mons,
1233290001Sglebius	int32_t mdays
1234290001Sglebius	)
1235290001Sglebius{
1236290001Sglebius	ntpcal_split tmp;
1237290001Sglebius	int32_t	     res;
1238290001Sglebius
1239290001Sglebius	if (mons) {
1240290001Sglebius		tmp = ntpcal_days_in_months(mons);
1241290001Sglebius		res = ntpcal_days_in_years(years + tmp.hi) + tmp.lo;
1242290001Sglebius	} else
1243290001Sglebius		res = ntpcal_days_in_years(years);
1244290001Sglebius	res += mdays;
1245290001Sglebius
1246290001Sglebius	return res;
1247290001Sglebius}
1248290001Sglebius
1249290001Sglebius/*
1250290001Sglebius *---------------------------------------------------------------------
1251290001Sglebius * Convert ELAPSED years/months/days of gregorian calendar to elapsed
1252290001Sglebius * days in year.
1253290001Sglebius *
1254290001Sglebius * Note: This will give the true difference to the start of the given year,
1255290001Sglebius * even if months & days are off-scale.
1256290001Sglebius *---------------------------------------------------------------------
1257290001Sglebius */
1258290001Sglebiusint32_t
1259290001Sglebiusntpcal_edate_to_yeardays(
1260290001Sglebius	int32_t years,
1261290001Sglebius	int32_t mons,
1262290001Sglebius	int32_t mdays
1263290001Sglebius	)
1264290001Sglebius{
1265290001Sglebius	ntpcal_split tmp;
1266290001Sglebius
1267290001Sglebius	if (0 <= mons && mons < 12) {
1268290001Sglebius		years += 1;
1269290001Sglebius		mdays += real_month_table[is_leapyear(years)][mons];
1270290001Sglebius	} else {
1271290001Sglebius		tmp = ntpcal_days_in_months(mons);
1272290001Sglebius		mdays += tmp.lo
1273290001Sglebius		       + ntpcal_days_in_years(years + tmp.hi)
1274290001Sglebius		       - ntpcal_days_in_years(years);
1275290001Sglebius	}
1276290001Sglebius
1277290001Sglebius	return mdays;
1278290001Sglebius}
1279290001Sglebius
1280290001Sglebius/*
1281290001Sglebius *---------------------------------------------------------------------
1282290001Sglebius * Convert elapsed days and the hour/minute/second information into
1283290001Sglebius * total seconds.
1284290001Sglebius *
1285290001Sglebius * If 'isvalid' is not NULL, do a range check on the time specification
1286290001Sglebius * and tell if the time input is in the normal range, permitting for a
1287290001Sglebius * single leapsecond.
1288290001Sglebius *---------------------------------------------------------------------
1289290001Sglebius */
1290290001Sglebiusint32_t
1291290001Sglebiusntpcal_etime_to_seconds(
1292290001Sglebius	int32_t hours,
1293290001Sglebius	int32_t minutes,
1294290001Sglebius	int32_t seconds
1295290001Sglebius	)
1296290001Sglebius{
1297290001Sglebius	int32_t res;
1298290001Sglebius
1299290001Sglebius	res = (hours * MINSPERHR + minutes) * SECSPERMIN + seconds;
1300290001Sglebius
1301290001Sglebius	return res;
1302290001Sglebius}
1303290001Sglebius
1304290001Sglebius/*
1305290001Sglebius *---------------------------------------------------------------------
1306290001Sglebius * Convert the date part of a 'struct tm' (that is, year, month,
1307290001Sglebius * day-of-month) into the RD of that day.
1308290001Sglebius *---------------------------------------------------------------------
1309290001Sglebius */
1310290001Sglebiusint32_t
1311290001Sglebiusntpcal_tm_to_rd(
1312290001Sglebius	const struct tm *utm
1313290001Sglebius	)
1314290001Sglebius{
1315290001Sglebius	return ntpcal_edate_to_eradays(utm->tm_year + 1899,
1316290001Sglebius				       utm->tm_mon,
1317290001Sglebius				       utm->tm_mday - 1) + 1;
1318290001Sglebius}
1319290001Sglebius
1320290001Sglebius/*
1321290001Sglebius *---------------------------------------------------------------------
1322290001Sglebius * Convert the date part of a 'struct calendar' (that is, year, month,
1323290001Sglebius * day-of-month) into the RD of that day.
1324290001Sglebius *---------------------------------------------------------------------
1325290001Sglebius */
1326290001Sglebiusint32_t
1327290001Sglebiusntpcal_date_to_rd(
1328290001Sglebius	const struct calendar *jd
1329290001Sglebius	)
1330290001Sglebius{
1331290001Sglebius	return ntpcal_edate_to_eradays((int32_t)jd->year - 1,
1332290001Sglebius				       (int32_t)jd->month - 1,
1333290001Sglebius				       (int32_t)jd->monthday - 1) + 1;
1334290001Sglebius}
1335290001Sglebius
1336290001Sglebius/*
1337290001Sglebius *---------------------------------------------------------------------
1338290001Sglebius * convert a year number to rata die of year start
1339290001Sglebius *---------------------------------------------------------------------
1340290001Sglebius */
1341290001Sglebiusint32_t
1342290001Sglebiusntpcal_year_to_ystart(
1343290001Sglebius	int32_t year
1344290001Sglebius	)
1345290001Sglebius{
1346290001Sglebius	return ntpcal_days_in_years(year - 1) + 1;
1347290001Sglebius}
1348290001Sglebius
1349290001Sglebius/*
1350290001Sglebius *---------------------------------------------------------------------
1351290001Sglebius * For a given RD, get the RD of the associated year start,
1352290001Sglebius * that is, the RD of the last January,1st on or before that day.
1353290001Sglebius *---------------------------------------------------------------------
1354290001Sglebius */
1355290001Sglebiusint32_t
1356290001Sglebiusntpcal_rd_to_ystart(
1357290001Sglebius	int32_t rd
1358290001Sglebius	)
1359290001Sglebius{
1360290001Sglebius	/*
1361290001Sglebius	 * Rather simple exercise: split the day number into elapsed
1362290001Sglebius	 * years and elapsed days, then remove the elapsed days from the
1363290001Sglebius	 * input value. Nice'n sweet...
1364290001Sglebius	 */
1365290001Sglebius	return rd - ntpcal_split_eradays(rd - 1, NULL).lo;
1366290001Sglebius}
1367290001Sglebius
1368290001Sglebius/*
1369290001Sglebius *---------------------------------------------------------------------
1370290001Sglebius * For a given RD, get the RD of the associated month start.
1371290001Sglebius *---------------------------------------------------------------------
1372290001Sglebius */
1373290001Sglebiusint32_t
1374290001Sglebiusntpcal_rd_to_mstart(
1375290001Sglebius	int32_t rd
1376290001Sglebius	)
1377290001Sglebius{
1378290001Sglebius	ntpcal_split split;
1379290001Sglebius	int	     leaps;
1380290001Sglebius
1381290001Sglebius	split = ntpcal_split_eradays(rd - 1, &leaps);
1382290001Sglebius	split = ntpcal_split_yeardays(split.lo, leaps);
1383290001Sglebius
1384290001Sglebius	return rd - split.lo;
1385290001Sglebius}
1386290001Sglebius
1387290001Sglebius/*
1388290001Sglebius *---------------------------------------------------------------------
1389290001Sglebius * take a 'struct calendar' and get the seconds-of-day from it.
1390290001Sglebius *---------------------------------------------------------------------
1391290001Sglebius */
1392290001Sglebiusint32_t
1393290001Sglebiusntpcal_date_to_daysec(
1394290001Sglebius	const struct calendar *jd
1395290001Sglebius	)
1396290001Sglebius{
1397290001Sglebius	return ntpcal_etime_to_seconds(jd->hour, jd->minute,
1398290001Sglebius				       jd->second);
1399290001Sglebius}
1400290001Sglebius
1401290001Sglebius/*
1402290001Sglebius *---------------------------------------------------------------------
1403290001Sglebius * take a 'struct tm' and get the seconds-of-day from it.
1404290001Sglebius *---------------------------------------------------------------------
1405290001Sglebius */
1406290001Sglebiusint32_t
1407290001Sglebiusntpcal_tm_to_daysec(
1408290001Sglebius	const struct tm *utm
1409290001Sglebius	)
1410290001Sglebius{
1411290001Sglebius	return ntpcal_etime_to_seconds(utm->tm_hour, utm->tm_min,
1412290001Sglebius				       utm->tm_sec);
1413290001Sglebius}
1414290001Sglebius
1415290001Sglebius/*
1416290001Sglebius *---------------------------------------------------------------------
1417290001Sglebius * take a 'struct calendar' and convert it to a 'time_t'
1418290001Sglebius *---------------------------------------------------------------------
1419290001Sglebius */
1420290001Sglebiustime_t
1421290001Sglebiusntpcal_date_to_time(
1422290001Sglebius	const struct calendar *jd
1423290001Sglebius	)
1424290001Sglebius{
1425290001Sglebius	vint64  join;
1426290001Sglebius	int32_t days, secs;
1427290001Sglebius
1428290001Sglebius	days = ntpcal_date_to_rd(jd) - DAY_UNIX_STARTS;
1429290001Sglebius	secs = ntpcal_date_to_daysec(jd);
1430290001Sglebius	join = ntpcal_dayjoin(days, secs);
1431290001Sglebius
1432290001Sglebius	return vint64_to_time(&join);
1433290001Sglebius}
1434290001Sglebius
1435290001Sglebius
1436290001Sglebius/*
1437290001Sglebius * ==================================================================
1438290001Sglebius *
1439290001Sglebius * extended and unchecked variants of caljulian/caltontp
1440290001Sglebius *
1441290001Sglebius * ==================================================================
1442290001Sglebius */
1443290001Sglebiusint
1444290001Sglebiusntpcal_ntp64_to_date(
1445290001Sglebius	struct calendar *jd,
1446290001Sglebius	const vint64    *ntp
1447290001Sglebius	)
1448290001Sglebius{
1449290001Sglebius	ntpcal_split ds;
1450290001Sglebius
1451290001Sglebius	ds = ntpcal_daysplit(ntp);
1452290001Sglebius	ds.hi += ntpcal_daysec_to_date(jd, ds.lo);
1453290001Sglebius
1454290001Sglebius	return ntpcal_rd_to_date(jd, ds.hi + DAY_NTP_STARTS);
1455290001Sglebius}
1456290001Sglebius
1457290001Sglebiusint
1458290001Sglebiusntpcal_ntp_to_date(
1459290001Sglebius	struct calendar *jd,
1460290001Sglebius	uint32_t	 ntp,
1461290001Sglebius	const time_t	*piv
1462290001Sglebius	)
1463290001Sglebius{
1464290001Sglebius	vint64	ntp64;
1465290001Sglebius
1466290001Sglebius	/*
1467290001Sglebius	 * Unfold ntp time around current time into NTP domain. Split
1468290001Sglebius	 * into days and seconds, shift days into CE domain and
1469290001Sglebius	 * process the parts.
1470290001Sglebius	 */
1471290001Sglebius	ntp64 = ntpcal_ntp_to_ntp(ntp, piv);
1472290001Sglebius	return ntpcal_ntp64_to_date(jd, &ntp64);
1473290001Sglebius}
1474290001Sglebius
1475290001Sglebius
1476290001Sglebiusvint64
1477290001Sglebiusntpcal_date_to_ntp64(
1478290001Sglebius	const struct calendar *jd
1479290001Sglebius	)
1480290001Sglebius{
1481290001Sglebius	/*
1482290001Sglebius	 * Convert date to NTP. Ignore yearday, use d/m/y only.
1483290001Sglebius	 */
1484290001Sglebius	return ntpcal_dayjoin(ntpcal_date_to_rd(jd) - DAY_NTP_STARTS,
1485290001Sglebius			      ntpcal_date_to_daysec(jd));
1486290001Sglebius}
1487290001Sglebius
1488290001Sglebius
1489290001Sglebiusuint32_t
1490290001Sglebiusntpcal_date_to_ntp(
1491290001Sglebius	const struct calendar *jd
1492290001Sglebius	)
1493290001Sglebius{
1494290001Sglebius	/*
1495290001Sglebius	 * Get lower half of 64-bit NTP timestamp from date/time.
1496290001Sglebius	 */
1497290001Sglebius	return ntpcal_date_to_ntp64(jd).d_s.lo;
1498290001Sglebius}
1499290001Sglebius
1500290001Sglebius
1501290001Sglebius
1502290001Sglebius/*
1503290001Sglebius * ==================================================================
1504290001Sglebius *
1505290001Sglebius * day-of-week calculations
1506290001Sglebius *
1507290001Sglebius * ==================================================================
1508290001Sglebius */
1509290001Sglebius/*
1510290001Sglebius * Given a RataDie and a day-of-week, calculate a RDN that is reater-than,
1511290001Sglebius * greater-or equal, closest, less-or-equal or less-than the given RDN
1512290001Sglebius * and denotes the given day-of-week
1513290001Sglebius */
1514290001Sglebiusint32_t
1515290001Sglebiusntpcal_weekday_gt(
1516290001Sglebius	int32_t rdn,
1517290001Sglebius	int32_t dow
1518290001Sglebius	)
1519290001Sglebius{
1520290001Sglebius	return ntpcal_periodic_extend(rdn+1, dow, 7);
1521290001Sglebius}
1522290001Sglebius
1523290001Sglebiusint32_t
1524290001Sglebiusntpcal_weekday_ge(
1525290001Sglebius	int32_t rdn,
1526290001Sglebius	int32_t dow
1527290001Sglebius	)
1528290001Sglebius{
1529290001Sglebius	return ntpcal_periodic_extend(rdn, dow, 7);
1530290001Sglebius}
1531290001Sglebius
1532290001Sglebiusint32_t
1533290001Sglebiusntpcal_weekday_close(
1534290001Sglebius	int32_t rdn,
1535290001Sglebius	int32_t dow
1536290001Sglebius	)
1537290001Sglebius{
1538290001Sglebius	return ntpcal_periodic_extend(rdn-3, dow, 7);
1539290001Sglebius}
1540290001Sglebius
1541290001Sglebiusint32_t
1542290001Sglebiusntpcal_weekday_le(
1543290001Sglebius	int32_t rdn,
1544290001Sglebius	int32_t dow
1545290001Sglebius	)
1546290001Sglebius{
1547290001Sglebius	return ntpcal_periodic_extend(rdn, dow, -7);
1548290001Sglebius}
1549290001Sglebius
1550290001Sglebiusint32_t
1551290001Sglebiusntpcal_weekday_lt(
1552290001Sglebius	int32_t rdn,
1553290001Sglebius	int32_t dow
1554290001Sglebius	)
1555290001Sglebius{
1556290001Sglebius	return ntpcal_periodic_extend(rdn-1, dow, -7);
1557290001Sglebius}
1558290001Sglebius
1559290001Sglebius/*
1560290001Sglebius * ==================================================================
1561290001Sglebius *
1562290001Sglebius * ISO week-calendar conversions
1563290001Sglebius *
1564290001Sglebius * The ISO8601 calendar defines a calendar of years, weeks and weekdays.
1565290001Sglebius * It is related to the Gregorian calendar, and a ISO year starts at the
1566290001Sglebius * Monday closest to Jan,1st of the corresponding Gregorian year.  A ISO
1567290001Sglebius * calendar year has always 52 or 53 weeks, and like the Grogrian
1568290001Sglebius * calendar the ISO8601 calendar repeats itself every 400 years, or
1569290001Sglebius * 146097 days, or 20871 weeks.
1570290001Sglebius *
1571290001Sglebius * While it is possible to write ISO calendar functions based on the
1572290001Sglebius * Gregorian calendar functions, the following implementation takes a
1573290001Sglebius * different approach, based directly on years and weeks.
1574290001Sglebius *
1575290001Sglebius * Analysis of the tabulated data shows that it is not possible to
1576290001Sglebius * interpolate from years to weeks over a full 400 year range; cyclic
1577290001Sglebius * shifts over 400 years do not provide a solution here. But it *is*
1578290001Sglebius * possible to interpolate over every single century of the 400-year
1579290001Sglebius * cycle. (The centennial leap year rule seems to be the culprit here.)
1580290001Sglebius *
1581290001Sglebius * It can be shown that a conversion from years to weeks can be done
1582290001Sglebius * using a linear transformation of the form
1583290001Sglebius *
1584290001Sglebius *   w = floor( y * a + b )
1585290001Sglebius *
1586290001Sglebius * where the slope a must hold to
1587290001Sglebius *
1588290001Sglebius *  52.1780821918 <= a < 52.1791044776
1589290001Sglebius *
1590290001Sglebius * and b must be chosen according to the selected slope and the number
1591290001Sglebius * of the century in a 400-year period.
1592290001Sglebius *
1593290001Sglebius * The inverse calculation can also be done in this way. Careful scaling
1594290001Sglebius * provides an unlimited set of integer coefficients a,k,b that enable
1595290001Sglebius * us to write the calulation in the form
1596290001Sglebius *
1597290001Sglebius *   w = (y * a	 + b ) / k
1598290001Sglebius *   y = (w * a' + b') / k'
1599290001Sglebius *
1600290001Sglebius * In this implementation the values of k and k' are chosen to be
1601290001Sglebius * smallest possible powers of two, so the division can be implemented
1602290001Sglebius * as shifts if the optimiser chooses to do so.
1603290001Sglebius *
1604290001Sglebius * ==================================================================
1605290001Sglebius */
1606290001Sglebius
1607290001Sglebius/*
1608290001Sglebius * Given a number of elapsed (ISO-)years since the begin of the
1609290001Sglebius * christian era, return the number of elapsed weeks corresponding to
1610290001Sglebius * the number of years.
1611290001Sglebius */
1612290001Sglebiusint32_t
1613290001Sglebiusisocal_weeks_in_years(
1614290001Sglebius	int32_t years
1615290001Sglebius	)
1616290001Sglebius{
1617290001Sglebius	/*
1618290001Sglebius	 * use: w = (y * 53431 + b[c]) / 1024 as interpolation
1619290001Sglebius	 */
1620290001Sglebius	static const uint16_t bctab[4] = { 157, 449, 597, 889 };
1621290001Sglebius
1622290001Sglebius	int32_t  cs, cw;
1623290001Sglebius	uint32_t cc, ci, yu, sflag;
1624290001Sglebius
1625290001Sglebius	sflag = int32_sflag(years);
1626290001Sglebius	yu    = int32_to_uint32_2cpl(years);
1627290001Sglebius
1628290001Sglebius	/* split off centuries, using floor division */
1629290001Sglebius	cc  = sflag ^ ((sflag ^ yu) / 100u);
1630290001Sglebius	yu -= cc * 100u;
1631290001Sglebius
1632290001Sglebius	/* calculate century cycles shift and cycle index:
1633290001Sglebius	 * Assuming a century is 5217 weeks, we have to add a cycle
1634290001Sglebius	 * shift that is 3 for every 4 centuries, because 3 of the four
1635290001Sglebius	 * centuries have 5218 weeks. So '(cc*3 + 1) / 4' is the actual
1636290001Sglebius	 * correction, and the second century is the defective one.
1637290001Sglebius	 *
1638290001Sglebius	 * Needs floor division by 4, which is done with masking and
1639290001Sglebius	 * shifting.
1640290001Sglebius	 */
1641290001Sglebius	ci = cc * 3u + 1;
1642290001Sglebius	cs = uint32_2cpl_to_int32(sflag ^ ((sflag ^ ci) / 4u));
1643290001Sglebius	ci = ci % 4u;
1644290001Sglebius
1645290001Sglebius	/* Get weeks in century. Can use plain division here as all ops
1646290001Sglebius	 * are >= 0,  and let the compiler sort out the possible
1647290001Sglebius	 * optimisations.
1648290001Sglebius	 */
1649290001Sglebius	cw = (yu * 53431u + bctab[ci]) / 1024u;
1650290001Sglebius
1651290001Sglebius	return uint32_2cpl_to_int32(cc) * 5217 + cs + cw;
1652290001Sglebius}
1653290001Sglebius
1654290001Sglebius/*
1655290001Sglebius * Given a number of elapsed weeks since the begin of the christian
1656290001Sglebius * era, split this number into the number of elapsed years in res.hi
1657290001Sglebius * and the excessive number of weeks in res.lo. (That is, res.lo is
1658290001Sglebius * the number of elapsed weeks in the remaining partial year.)
1659290001Sglebius */
1660290001Sglebiusntpcal_split
1661290001Sglebiusisocal_split_eraweeks(
1662290001Sglebius	int32_t weeks
1663290001Sglebius	)
1664290001Sglebius{
1665290001Sglebius	/*
1666290001Sglebius	 * use: y = (w * 157 + b[c]) / 8192 as interpolation
1667290001Sglebius	 */
1668290001Sglebius
1669290001Sglebius	static const uint16_t bctab[4] = { 85, 130, 17, 62 };
1670290001Sglebius
1671290001Sglebius	ntpcal_split res;
1672290001Sglebius	int32_t  cc, ci;
1673290001Sglebius	uint32_t sw, cy, Q, sflag;
1674290001Sglebius
1675290001Sglebius	/* Use two fast cycle-split divisions here. This is again
1676290001Sglebius	 * susceptible to internal overflow, so we check the range. This
1677290001Sglebius	 * still permits more than +/-20 million years, so this is
1678290001Sglebius	 * likely a pure academical problem.
1679290001Sglebius	 *
1680290001Sglebius	 * We want to execute '(weeks * 4 + 2) /% 20871' under floor
1681290001Sglebius	 * division rules in the first step.
1682290001Sglebius	 */
1683290001Sglebius	sflag = int32_sflag(weeks);
1684290001Sglebius	sw  = uint32_saturate(int32_to_uint32_2cpl(weeks), sflag);
1685290001Sglebius	sw  = 4u * sw + 2;
1686290001Sglebius	Q   = sflag ^ ((sflag ^ sw) / GREGORIAN_CYCLE_WEEKS);
1687290001Sglebius	sw -= Q * GREGORIAN_CYCLE_WEEKS;
1688290001Sglebius	ci  = Q % 4u;
1689290001Sglebius	cc  = uint32_2cpl_to_int32(Q);
1690290001Sglebius
1691290001Sglebius	/* Split off years; sw >= 0 here! The scaled weeks in the years
1692290001Sglebius	 * are scaled up by 157 afterwards.
1693290001Sglebius	 */
1694290001Sglebius	sw  = (sw / 4u) * 157u + bctab[ci];
1695290001Sglebius	cy  = sw / 8192u;	/* ws >> 13 , let the compiler sort it out */
1696290001Sglebius	sw  = sw % 8192u;	/* ws & 8191, let the compiler sort it out */
1697290001Sglebius
1698290001Sglebius	/* assemble elapsed years and downscale the elapsed weeks in
1699290001Sglebius	 * the year.
1700290001Sglebius	 */
1701290001Sglebius	res.hi = 100*cc + cy;
1702290001Sglebius	res.lo = sw / 157u;
1703290001Sglebius
1704290001Sglebius	return res;
1705290001Sglebius}
1706290001Sglebius
1707290001Sglebius/*
1708290001Sglebius * Given a second in the NTP time scale and a pivot, expand the NTP
1709290001Sglebius * time stamp around the pivot and convert into an ISO calendar time
1710290001Sglebius * stamp.
1711290001Sglebius */
1712290001Sglebiusint
1713290001Sglebiusisocal_ntp64_to_date(
1714290001Sglebius	struct isodate *id,
1715290001Sglebius	const vint64   *ntp
1716290001Sglebius	)
1717290001Sglebius{
1718290001Sglebius	ntpcal_split ds;
1719290001Sglebius	int32_t      ts[3];
1720290001Sglebius	uint32_t     uw, ud, sflag;
1721290001Sglebius
1722290001Sglebius	/*
1723290001Sglebius	 * Split NTP time into days and seconds, shift days into CE
1724290001Sglebius	 * domain and process the parts.
1725290001Sglebius	 */
1726290001Sglebius	ds = ntpcal_daysplit(ntp);
1727290001Sglebius
1728290001Sglebius	/* split time part */
1729290001Sglebius	ds.hi += priv_timesplit(ts, ds.lo);
1730290001Sglebius	id->hour   = (uint8_t)ts[0];
1731290001Sglebius	id->minute = (uint8_t)ts[1];
1732290001Sglebius	id->second = (uint8_t)ts[2];
1733290001Sglebius
1734290001Sglebius	/* split days into days and weeks, using floor division in unsigned */
1735290001Sglebius	ds.hi += DAY_NTP_STARTS - 1; /* shift from NTP to RDN */
1736290001Sglebius	sflag = int32_sflag(ds.hi);
1737290001Sglebius	ud  = int32_to_uint32_2cpl(ds.hi);
1738290001Sglebius	uw  = sflag ^ ((sflag ^ ud) / DAYSPERWEEK);
1739290001Sglebius	ud -= uw * DAYSPERWEEK;
1740290001Sglebius	ds.hi = uint32_2cpl_to_int32(uw);
1741290001Sglebius	ds.lo = ud;
1742290001Sglebius
1743290001Sglebius	id->weekday = (uint8_t)ds.lo + 1;	/* weekday result    */
1744290001Sglebius
1745290001Sglebius	/* get year and week in year */
1746290001Sglebius	ds = isocal_split_eraweeks(ds.hi);	/* elapsed years&week*/
1747290001Sglebius	id->year = (uint16_t)ds.hi + 1;		/* shift to current  */
1748290001Sglebius	id->week = (uint8_t )ds.lo + 1;
1749290001Sglebius
1750290001Sglebius	return (ds.hi >= 0 && ds.hi < 0x0000FFFF);
1751290001Sglebius}
1752290001Sglebius
1753290001Sglebiusint
1754290001Sglebiusisocal_ntp_to_date(
1755290001Sglebius	struct isodate *id,
1756290001Sglebius	uint32_t	ntp,
1757290001Sglebius	const time_t   *piv
1758290001Sglebius	)
1759290001Sglebius{
1760290001Sglebius	vint64	ntp64;
1761290001Sglebius
1762290001Sglebius	/*
1763290001Sglebius	 * Unfold ntp time around current time into NTP domain, then
1764290001Sglebius	 * convert the full time stamp.
1765290001Sglebius	 */
1766290001Sglebius	ntp64 = ntpcal_ntp_to_ntp(ntp, piv);
1767290001Sglebius	return isocal_ntp64_to_date(id, &ntp64);
1768290001Sglebius}
1769290001Sglebius
1770290001Sglebius/*
1771290001Sglebius * Convert a ISO date spec into a second in the NTP time scale,
1772290001Sglebius * properly truncated to 32 bit.
1773290001Sglebius */
1774290001Sglebiusvint64
1775290001Sglebiusisocal_date_to_ntp64(
1776290001Sglebius	const struct isodate *id
1777290001Sglebius	)
1778290001Sglebius{
1779290001Sglebius	int32_t weeks, days, secs;
1780290001Sglebius
1781290001Sglebius	weeks = isocal_weeks_in_years((int32_t)id->year - 1)
1782290001Sglebius	      + (int32_t)id->week - 1;
1783290001Sglebius	days = weeks * 7 + (int32_t)id->weekday;
1784290001Sglebius	/* days is RDN of ISO date now */
1785290001Sglebius	secs = ntpcal_etime_to_seconds(id->hour, id->minute, id->second);
1786290001Sglebius
1787290001Sglebius	return ntpcal_dayjoin(days - DAY_NTP_STARTS, secs);
1788290001Sglebius}
1789290001Sglebius
1790290001Sglebiusuint32_t
1791290001Sglebiusisocal_date_to_ntp(
1792290001Sglebius	const struct isodate *id
1793290001Sglebius	)
1794290001Sglebius{
1795290001Sglebius	/*
1796290001Sglebius	 * Get lower half of 64-bit NTP timestamp from date/time.
1797290001Sglebius	 */
1798290001Sglebius	return isocal_date_to_ntp64(id).d_s.lo;
1799290001Sglebius}
1800290001Sglebius
1801290001Sglebius/* -*-EOF-*- */
1802