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 potentially suffer from
44 * an internal overflow, and these are coded in a way that avoids
45 * it. All other functions do not suffer from internal overflow and
46 * simply return the result truncated to 32 bits.
47 */
48
49 #include <config.h>
50 #include <sys/types.h>
51
52 #include "ntp_types.h"
53 #include "ntp_calendar.h"
54 #include "ntp_stdlib.h"
55 #include "ntp_fp.h"
56 #include "ntp_unixtime.h"
57
58 #include "ntpd.h"
59 #include "lib_strbuf.h"
60
61 /* For now, let's take the conservative approach: if the target property
62 * macros are not defined, check a few well-known compiler/architecture
63 * settings. Default is to assume that the representation of signed
64 * integers is unknown and shift-arithmetic-right is not available.
65 */
66 #ifndef TARGET_HAS_2CPL
67 # if defined(__GNUC__)
68 # if defined(__i386__) || defined(__x86_64__) || defined(__arm__)
69 # define TARGET_HAS_2CPL 1
70 # else
71 # define TARGET_HAS_2CPL 0
72 # endif
73 # elif defined(_MSC_VER)
74 # if defined(_M_IX86) || defined(_M_X64) || defined(_M_ARM)
75 # define TARGET_HAS_2CPL 1
76 # else
77 # define TARGET_HAS_2CPL 0
78 # endif
79 # else
80 # define TARGET_HAS_2CPL 0
81 # endif
82 #endif
83
84 #ifndef TARGET_HAS_SAR
85 # define TARGET_HAS_SAR 0
86 #endif
87
88 #if !defined(HAVE_64BITREGS) && defined(UINT64_MAX) && (SIZE_MAX >= UINT64_MAX)
89 # define HAVE_64BITREGS
90 #endif
91
92 /*
93 *---------------------------------------------------------------------
94 * replacing the 'time()' function
95 *---------------------------------------------------------------------
96 */
97
98 static systime_func_ptr systime_func = &time;
99 static inline time_t now(void);
100
101
102 systime_func_ptr
ntpcal_set_timefunc(systime_func_ptr nfunc)103 ntpcal_set_timefunc(
104 systime_func_ptr nfunc
105 )
106 {
107 systime_func_ptr res;
108
109 res = systime_func;
110 if (NULL == nfunc)
111 nfunc = &time;
112 systime_func = nfunc;
113
114 return res;
115 }
116
117
118 static inline time_t
now(void)119 now(void)
120 {
121 return (*systime_func)(NULL);
122 }
123
124 /*
125 *---------------------------------------------------------------------
126 * Get sign extension mask and unsigned 2cpl rep for a signed integer
127 *---------------------------------------------------------------------
128 */
129
130 static inline uint32_t
int32_sflag(const int32_t v)131 int32_sflag(
132 const int32_t v)
133 {
134 # if TARGET_HAS_2CPL && TARGET_HAS_SAR && SIZEOF_INT >= 4
135
136 /* Let's assume that shift is the fastest way to get the sign
137 * extension of of a signed integer. This might not always be
138 * true, though -- On 8bit CPUs or machines without barrel
139 * shifter this will kill the performance. So we make sure
140 * we do this only if 'int' has at least 4 bytes.
141 */
142 return (uint32_t)(v >> 31);
143
144 # else
145
146 /* This should be a rather generic approach for getting a sign
147 * extension mask...
148 */
149 return UINT32_C(0) - (uint32_t)(v < 0);
150
151 # endif
152 }
153
154 static inline int32_t
uint32_2cpl_to_int32(const uint32_t vu)155 uint32_2cpl_to_int32(
156 const uint32_t vu)
157 {
158 int32_t v;
159
160 # if TARGET_HAS_2CPL
161
162 /* Just copy through the 32 bits from the unsigned value if
163 * we're on a two's complement target.
164 */
165 v = (int32_t)vu;
166
167 # else
168
169 /* Convert to signed integer, making sure signed integer
170 * overflow cannot happen. Again, the optimiser might or might
171 * not find out that this is just a copy of 32 bits on a target
172 * with two's complement representation for signed integers.
173 */
174 if (vu > INT32_MAX)
175 v = -(int32_t)(~vu) - 1;
176 else
177 v = (int32_t)vu;
178
179 # endif
180
181 return v;
182 }
183
184 /*
185 *---------------------------------------------------------------------
186 * Convert between 'time_t' and 'vint64'
187 *---------------------------------------------------------------------
188 */
189 vint64
time_to_vint64(const time_t * ptt)190 time_to_vint64(
191 const time_t * ptt
192 )
193 {
194 vint64 res;
195 time_t tt;
196
197 tt = *ptt;
198
199 # if SIZEOF_TIME_T <= 4
200
201 res.D_s.hi = 0;
202 if (tt < 0) {
203 res.D_s.lo = (uint32_t)-tt;
204 M_NEG(res.D_s.hi, res.D_s.lo);
205 } else {
206 res.D_s.lo = (uint32_t)tt;
207 }
208
209 # elif defined(HAVE_INT64)
210
211 res.q_s = tt;
212
213 # else
214 /*
215 * shifting negative signed quantities is compiler-dependent, so
216 * we better avoid it and do it all manually. And shifting more
217 * than the width of a quantity is undefined. Also a don't do!
218 */
219 if (tt < 0) {
220 tt = -tt;
221 res.D_s.lo = (uint32_t)tt;
222 res.D_s.hi = (uint32_t)(tt >> 32);
223 M_NEG(res.D_s.hi, res.D_s.lo);
224 } else {
225 res.D_s.lo = (uint32_t)tt;
226 res.D_s.hi = (uint32_t)(tt >> 32);
227 }
228
229 # endif
230
231 return res;
232 }
233
234
235 time_t
vint64_to_time(const vint64 * tv)236 vint64_to_time(
237 const vint64 *tv
238 )
239 {
240 time_t res;
241
242 # if SIZEOF_TIME_T <= 4
243
244 res = (time_t)tv->D_s.lo;
245
246 # elif defined(HAVE_INT64)
247
248 res = (time_t)tv->q_s;
249
250 # else
251
252 res = ((time_t)tv->d_s.hi << 32) | tv->D_s.lo;
253
254 # endif
255
256 return res;
257 }
258
259 /*
260 *---------------------------------------------------------------------
261 * Get the build date & time
262 *---------------------------------------------------------------------
263 */
264 int
ntpcal_get_build_date(struct calendar * jd)265 ntpcal_get_build_date(
266 struct calendar * jd
267 )
268 {
269 /* The C standard tells us the format of '__DATE__':
270 *
271 * __DATE__ The date of translation of the preprocessing
272 * translation unit: a character string literal of the form "Mmm
273 * dd yyyy", where the names of the months are the same as those
274 * generated by the asctime function, and the first character of
275 * dd is a space character if the value is less than 10. If the
276 * date of translation is not available, an
277 * implementation-defined valid date shall be supplied.
278 *
279 * __TIME__ The time of translation of the preprocessing
280 * translation unit: a character string literal of the form
281 * "hh:mm:ss" as in the time generated by the asctime
282 * function. If the time of translation is not available, an
283 * implementation-defined valid time shall be supplied.
284 *
285 * Note that MSVC declares DATE and TIME to be in the local time
286 * zone, while neither the C standard nor the GCC docs make any
287 * statement about this. As a result, we may be +/-12hrs off
288 * UTC. But for practical purposes, this should not be a
289 * problem.
290 *
291 */
292 # ifdef MKREPRO_DATE
293 static const char build[] = MKREPRO_TIME "/" MKREPRO_DATE;
294 # else
295 static const char build[] = __TIME__ "/" __DATE__;
296 # endif
297 static const char mlist[] = "JanFebMarAprMayJunJulAugSepOctNovDec";
298
299 char monstr[4];
300 const char * cp;
301 unsigned short hour, minute, second, day, year;
302 /* Note: The above quantities are used for sscanf 'hu' format,
303 * so using 'uint16_t' is contra-indicated!
304 */
305
306 # ifdef DEBUG
307 static int ignore = 0;
308 # endif
309
310 ZERO(*jd);
311 jd->year = 1970;
312 jd->month = 1;
313 jd->monthday = 1;
314
315 # ifdef DEBUG
316 /* check environment if build date should be ignored */
317 if (0 == ignore) {
318 const char * envstr;
319 envstr = getenv("NTPD_IGNORE_BUILD_DATE");
320 ignore = 1 + (envstr && (!*envstr || !strcasecmp(envstr, "yes")));
321 }
322 if (ignore > 1)
323 return FALSE;
324 # endif
325
326 if (6 == sscanf(build, "%hu:%hu:%hu/%3s %hu %hu",
327 &hour, &minute, &second, monstr, &day, &year)) {
328 cp = strstr(mlist, monstr);
329 if (NULL != cp) {
330 jd->year = year;
331 jd->month = (uint8_t)((cp - mlist) / 3 + 1);
332 jd->monthday = (uint8_t)day;
333 jd->hour = (uint8_t)hour;
334 jd->minute = (uint8_t)minute;
335 jd->second = (uint8_t)second;
336
337 return TRUE;
338 }
339 }
340
341 return FALSE;
342 }
343
344
345 /*
346 *---------------------------------------------------------------------
347 * basic calendar stuff
348 *---------------------------------------------------------------------
349 */
350
351 /*
352 * Some notes on the terminology:
353 *
354 * We use the proleptic Gregorian calendar, which is the Gregorian
355 * calendar extended in both directions ad infinitum. This totally
356 * disregards the fact that this calendar was invented in 1582, and
357 * was adopted at various dates over the world; sometimes even after
358 * the start of the NTP epoch.
359 *
360 * Normally date parts are given as current cycles, while time parts
361 * are given as elapsed cycles:
362 *
363 * 1970-01-01/03:04:05 means 'IN the 1970st. year, IN the first month,
364 * ON the first day, with 3hrs, 4minutes and 5 seconds elapsed.
365 *
366 * The basic calculations for this calendar implementation deal with
367 * ELAPSED date units, which is the number of full years, full months
368 * and full days before a date: 1970-01-01 would be (1969, 0, 0) in
369 * that notation.
370 *
371 * To ease the numeric computations, month and day values outside the
372 * normal range are acceptable: 2001-03-00 will be treated as the day
373 * before 2001-03-01, 2000-13-32 will give the same result as
374 * 2001-02-01 and so on.
375 *
376 * 'rd' or 'RD' is used as an abbreviation for the latin 'rata die'
377 * (day number). This is the number of days elapsed since 0000-12-31
378 * in the proleptic Gregorian calendar. The begin of the Christian Era
379 * (0001-01-01) is RD(1).
380 */
381
382 /*
383 * ====================================================================
384 *
385 * General algorithmic stuff
386 *
387 * ====================================================================
388 */
389
390 /*
391 *---------------------------------------------------------------------
392 * fast modulo 7 operations (floor/mathematical convention)
393 *---------------------------------------------------------------------
394 */
395 int
u32mod7(uint32_t x)396 u32mod7(
397 uint32_t x
398 )
399 {
400 /* This is a combination of tricks from "Hacker's Delight" with
401 * some modifications, like a multiplication that rounds up to
402 * drop the final adjustment stage.
403 *
404 * Do a partial reduction by digit sum to keep the value in the
405 * range permitted for the mul/shift stage. There are several
406 * possible and absolutely equivalent shift/mask combinations;
407 * this one is ARM-friendly because of a mask that fits into 16
408 * bit.
409 */
410 x = (x >> 15) + (x & UINT32_C(0x7FFF));
411 /* Take reminder as (mod 8) by mul/shift. Since the multiplier
412 * was calculated using ceil() instead of floor(), it skips the
413 * value '7' properly.
414 * M <- ceil(ldexp(8/7, 29))
415 */
416 return (int)((x * UINT32_C(0x24924925)) >> 29);
417 }
418
419 int
i32mod7(int32_t x)420 i32mod7(
421 int32_t x
422 )
423 {
424 /* We add (2**32 - 2**32 % 7), which is (2**32 - 4), to negative
425 * numbers to map them into the postive range. Only the term '-4'
426 * survives, obviously.
427 */
428 uint32_t ux = (uint32_t)x;
429 return u32mod7((x < 0) ? (ux - 4u) : ux);
430 }
431
432 uint32_t
i32fmod(int32_t x,uint32_t d)433 i32fmod(
434 int32_t x,
435 uint32_t d
436 )
437 {
438 uint32_t ux = (uint32_t)x;
439 uint32_t sf = UINT32_C(0) - (x < 0);
440 ux = (sf ^ ux ) % d;
441 return (d & sf) + (sf ^ ux);
442 }
443
444 /*
445 *---------------------------------------------------------------------
446 * Do a periodic extension of 'value' around 'pivot' with a period of
447 * 'cycle'.
448 *
449 * The result 'res' is a number that holds to the following properties:
450 *
451 * 1) res MOD cycle == value MOD cycle
452 * 2) pivot <= res < pivot + cycle
453 * (replace </<= with >/>= for negative cycles)
454 *
455 * where 'MOD' denotes the modulo operator for FLOOR DIVISION, which
456 * is not the same as the '%' operator in C: C requires division to be
457 * a truncated division, where remainder and dividend have the same
458 * sign if the remainder is not zero, whereas floor division requires
459 * divider and modulus to have the same sign for a non-zero modulus.
460 *
461 * This function has some useful applications:
462 *
463 * + let Y be a calendar year and V a truncated 2-digit year: then
464 * periodic_extend(Y-50, V, 100)
465 * is the closest expansion of the truncated year with respect to
466 * the full year, that is a 4-digit year with a difference of less
467 * than 50 years to the year Y. ("century unfolding")
468 *
469 * + let T be a UN*X time stamp and V be seconds-of-day: then
470 * perodic_extend(T-43200, V, 86400)
471 * is a time stamp that has the same seconds-of-day as the input
472 * value, with an absolute difference to T of <= 12hrs. ("day
473 * unfolding")
474 *
475 * + Wherever you have a truncated periodic value and a non-truncated
476 * base value and you want to match them somehow...
477 *
478 * Basically, the function delivers 'pivot + (value - pivot) % cycle',
479 * but the implementation takes some pains to avoid internal signed
480 * integer overflows in the '(value - pivot) % cycle' part and adheres
481 * to the floor division convention.
482 *
483 * If 64bit scalars where available on all intended platforms, writing a
484 * version that uses 64 bit ops would be easy; writing a general
485 * division routine for 64bit ops on a platform that can only do
486 * 32/16bit divisions and is still performant is a bit more
487 * difficult. Since most usecases can be coded in a way that does only
488 * require the 32bit version a 64bit version is NOT provided here.
489 *---------------------------------------------------------------------
490 */
491 int32_t
ntpcal_periodic_extend(int32_t pivot,int32_t value,int32_t cycle)492 ntpcal_periodic_extend(
493 int32_t pivot,
494 int32_t value,
495 int32_t cycle
496 )
497 {
498 /* Implement a 4-quadrant modulus calculation by 2 2-quadrant
499 * branches, one for positive and one for negative dividers.
500 * Everything else can be handled by bit level logic and
501 * conditional one's complement arithmetic. By convention, we
502 * assume
503 *
504 * x % b == 0 if |b| < 2
505 *
506 * that is, we don't actually divide for cycles of -1,0,1 and
507 * return the pivot value in that case.
508 */
509 uint32_t uv = (uint32_t)value;
510 uint32_t up = (uint32_t)pivot;
511 uint32_t uc, sf;
512
513 if (cycle > 1)
514 {
515 uc = (uint32_t)cycle;
516 sf = UINT32_C(0) - (value < pivot);
517
518 uv = sf ^ (uv - up);
519 uv %= uc;
520 pivot += (uc & sf) + (sf ^ uv);
521 }
522 else if (cycle < -1)
523 {
524 uc = ~(uint32_t)cycle + 1;
525 sf = UINT32_C(0) - (value > pivot);
526
527 uv = sf ^ (up - uv);
528 uv %= uc;
529 pivot -= (uc & sf) + (sf ^ uv);
530 }
531 return pivot;
532 }
533
534 /*---------------------------------------------------------------------
535 * Note to the casual reader
536 *
537 * In the next two functions you will find (or would have found...)
538 * the expression
539 *
540 * res.Q_s -= 0x80000000;
541 *
542 * There was some ruckus about a possible programming error due to
543 * integer overflow and sign propagation.
544 *
545 * This assumption is based on a lack of understanding of the C
546 * standard. (Though this is admittedly not one of the most 'natural'
547 * aspects of the 'C' language and easily to get wrong.)
548 *
549 * see
550 * http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1570.pdf
551 * "ISO/IEC 9899:201x Committee Draft — April 12, 2011"
552 * 6.4.4.1 Integer constants, clause 5
553 *
554 * why there is no sign extension/overflow problem here.
555 *
556 * But to ease the minds of the doubtful, I added back the 'u' qualifiers
557 * that somehow got lost over the last years.
558 */
559
560
561 /*
562 *---------------------------------------------------------------------
563 * Convert a timestamp in NTP scale to a 64bit seconds value in the UN*X
564 * scale with proper epoch unfolding around a given pivot or the current
565 * system time. This function happily accepts negative pivot values as
566 * timestamps before 1970-01-01, so be aware of possible trouble on
567 * platforms with 32bit 'time_t'!
568 *
569 * This is also a periodic extension, but since the cycle is 2^32 and
570 * the shift is 2^31, we can do some *very* fast math without explicit
571 * divisions.
572 *---------------------------------------------------------------------
573 */
574 vint64
ntpcal_ntp_to_time(uint32_t ntp,const time_t * pivot)575 ntpcal_ntp_to_time(
576 uint32_t ntp,
577 const time_t * pivot
578 )
579 {
580 vint64 res;
581
582 # if defined(HAVE_INT64)
583
584 res.q_s = (pivot != NULL)
585 ? *pivot
586 : now();
587 res.Q_s -= 0x80000000u; /* unshift of half range */
588 ntp -= (uint32_t)JAN_1970; /* warp into UN*X domain */
589 ntp -= res.D_s.lo; /* cycle difference */
590 res.Q_s += (uint64_t)ntp; /* get expanded time */
591
592 # else /* no 64bit scalars */
593
594 time_t tmp;
595
596 tmp = (pivot != NULL)
597 ? *pivot
598 : now();
599 res = time_to_vint64(&tmp);
600 M_SUB(res.D_s.hi, res.D_s.lo, 0, 0x80000000u);
601 ntp -= (uint32_t)JAN_1970; /* warp into UN*X domain */
602 ntp -= res.D_s.lo; /* cycle difference */
603 M_ADD(res.D_s.hi, res.D_s.lo, 0, ntp);
604
605 # endif /* no 64bit scalars */
606
607 return res;
608 }
609
610 /*
611 *---------------------------------------------------------------------
612 * Convert a timestamp in NTP scale to a 64bit seconds value in the NTP
613 * scale with proper epoch unfolding around a given pivot or the current
614 * system time.
615 *
616 * Note: The pivot must be given in the UN*X time domain!
617 *
618 * This is also a periodic extension, but since the cycle is 2^32 and
619 * the shift is 2^31, we can do some *very* fast math without explicit
620 * divisions.
621 *---------------------------------------------------------------------
622 */
623 vint64
ntpcal_ntp_to_ntp(uint32_t ntp,const time_t * pivot)624 ntpcal_ntp_to_ntp(
625 uint32_t ntp,
626 const time_t *pivot
627 )
628 {
629 vint64 res;
630
631 # if defined(HAVE_INT64)
632
633 res.q_s = (pivot)
634 ? *pivot
635 : now();
636 res.Q_s -= 0x80000000u; /* unshift of half range */
637 res.Q_s += (uint32_t)JAN_1970; /* warp into NTP domain */
638 ntp -= res.D_s.lo; /* cycle difference */
639 res.Q_s += (uint64_t)ntp; /* get expanded time */
640
641 # else /* no 64bit scalars */
642
643 time_t tmp;
644
645 tmp = (pivot)
646 ? *pivot
647 : now();
648 res = time_to_vint64(&tmp);
649 M_SUB(res.D_s.hi, res.D_s.lo, 0, 0x80000000u);
650 M_ADD(res.D_s.hi, res.D_s.lo, 0, (uint32_t)JAN_1970);/*into NTP */
651 ntp -= res.D_s.lo; /* cycle difference */
652 M_ADD(res.D_s.hi, res.D_s.lo, 0, ntp);
653
654 # endif /* no 64bit scalars */
655
656 return res;
657 }
658
659
660 /*
661 * ====================================================================
662 *
663 * Splitting values to composite entities
664 *
665 * ====================================================================
666 */
667
668 /*
669 *---------------------------------------------------------------------
670 * Split a 64bit seconds value into elapsed days in 'res.hi' and
671 * elapsed seconds since midnight in 'res.lo' using explicit floor
672 * division. This function happily accepts negative time values as
673 * timestamps before the respective epoch start.
674 *---------------------------------------------------------------------
675 */
676 ntpcal_split
ntpcal_daysplit(const vint64 * ts)677 ntpcal_daysplit(
678 const vint64 *ts
679 )
680 {
681 ntpcal_split res;
682 uint32_t Q, R;
683
684 # if defined(HAVE_64BITREGS)
685
686 /* Assume we have 64bit registers an can do a divison by
687 * constant reasonably fast using the one's complement trick..
688 */
689 uint64_t sf64 = (uint64_t)-(ts->q_s < 0);
690 Q = (uint32_t)(sf64 ^ ((sf64 ^ ts->Q_s) / SECSPERDAY));
691 R = (uint32_t)(ts->Q_s - Q * SECSPERDAY);
692
693 # elif defined(UINT64_MAX) && !defined(__arm__)
694
695 /* We rely on the compiler to do efficient 64bit divisions as
696 * good as possible. Which might or might not be true. At least
697 * for ARM CPUs, the sum-by-digit code in the next section is
698 * faster for many compilers. (This might change over time, but
699 * the 64bit-by-32bit division will never outperform the exact
700 * division by a substantial factor....)
701 */
702 if (ts->q_s < 0)
703 Q = ~(uint32_t)(~ts->Q_s / SECSPERDAY);
704 else
705 Q = (uint32_t)( ts->Q_s / SECSPERDAY);
706 R = ts->D_s.lo - Q * SECSPERDAY;
707
708 # else
709
710 /* We don't have 64bit regs. That hurts a bit.
711 *
712 * Here we use a mean trick to get away with just one explicit
713 * modulo operation and pure 32bit ops.
714 *
715 * Remember: 86400 <--> 128 * 675
716 *
717 * So we discard the lowest 7 bit and do an exact division by
718 * 675, modulo 2**32.
719 *
720 * First we shift out the lower 7 bits.
721 *
722 * Then we use a digit-wise pseudo-reduction, where a 'digit' is
723 * actually a 16-bit group. This is followed by a full reduction
724 * with a 'true' division step. This yields the modulus of the
725 * full 64bit value. The sign bit gets some extra treatment.
726 *
727 * Then we decrement the lower limb by that modulus, so it is
728 * exactly divisible by 675. [*]
729 *
730 * Then we multiply with the modular inverse of 675 (mod 2**32)
731 * and voila, we have the result.
732 *
733 * Special Thanks to Henry S. Warren and his "Hacker's delight"
734 * for giving that idea.
735 *
736 * (Note[*]: that's not the full truth. We would have to
737 * subtract the modulus from the full 64 bit number to get a
738 * number that is divisible by 675. But since we use the
739 * multiplicative inverse (mod 2**32) there's no reason to carry
740 * the subtraction into the upper bits!)
741 */
742 uint32_t al = ts->D_s.lo;
743 uint32_t ah = ts->D_s.hi;
744
745 /* shift out the lower 7 bits, smash sign bit */
746 al = (al >> 7) | (ah << 25);
747 ah = (ah >> 7) & 0x00FFFFFFu;
748
749 R = (ts->d_s.hi < 0) ? 239 : 0;/* sign bit value */
750 R += (al & 0xFFFF);
751 R += (al >> 16 ) * 61u; /* 2**16 % 675 */
752 R += (ah & 0xFFFF) * 346u; /* 2**32 % 675 */
753 R += (ah >> 16 ) * 181u; /* 2**48 % 675 */
754 R %= 675u; /* final reduction */
755 Q = (al - R) * 0x2D21C10Bu; /* modinv(675, 2**32) */
756 R = (R << 7) | (ts->d_s.lo & 0x07F);
757
758 # endif
759
760 res.hi = uint32_2cpl_to_int32(Q);
761 res.lo = R;
762
763 return res;
764 }
765
766 /*
767 *---------------------------------------------------------------------
768 * Split a 64bit seconds value into elapsed weeks in 'res.hi' and
769 * elapsed seconds since week start in 'res.lo' using explicit floor
770 * division. This function happily accepts negative time values as
771 * timestamps before the respective epoch start.
772 *---------------------------------------------------------------------
773 */
774 ntpcal_split
ntpcal_weeksplit(const vint64 * ts)775 ntpcal_weeksplit(
776 const vint64 *ts
777 )
778 {
779 ntpcal_split res;
780 uint32_t Q, R;
781
782 /* This is a very close relative to the day split function; for
783 * details, see there!
784 */
785
786 # if defined(HAVE_64BITREGS)
787
788 uint64_t sf64 = (uint64_t)-(ts->q_s < 0);
789 Q = (uint32_t)(sf64 ^ ((sf64 ^ ts->Q_s) / SECSPERWEEK));
790 R = (uint32_t)(ts->Q_s - Q * SECSPERWEEK);
791
792 # elif defined(UINT64_MAX) && !defined(__arm__)
793
794 if (ts->q_s < 0)
795 Q = ~(uint32_t)(~ts->Q_s / SECSPERWEEK);
796 else
797 Q = (uint32_t)( ts->Q_s / SECSPERWEEK);
798 R = ts->D_s.lo - Q * SECSPERWEEK;
799
800 # else
801
802 /* Remember: 7*86400 <--> 604800 <--> 128 * 4725 */
803 uint32_t al = ts->D_s.lo;
804 uint32_t ah = ts->D_s.hi;
805
806 al = (al >> 7) | (ah << 25);
807 ah = (ah >> 7) & 0x00FFFFFF;
808
809 R = (ts->d_s.hi < 0) ? 2264 : 0;/* sign bit value */
810 R += (al & 0xFFFF);
811 R += (al >> 16 ) * 4111u; /* 2**16 % 4725 */
812 R += (ah & 0xFFFF) * 3721u; /* 2**32 % 4725 */
813 R += (ah >> 16 ) * 2206u; /* 2**48 % 4725 */
814 R %= 4725u; /* final reduction */
815 Q = (al - R) * 0x98BBADDDu; /* modinv(4725, 2**32) */
816 R = (R << 7) | (ts->d_s.lo & 0x07F);
817
818 # endif
819
820 res.hi = uint32_2cpl_to_int32(Q);
821 res.lo = R;
822
823 return res;
824 }
825
826 /*
827 *---------------------------------------------------------------------
828 * Split a 32bit seconds value into h/m/s and excessive days. This
829 * function happily accepts negative time values as timestamps before
830 * midnight.
831 *---------------------------------------------------------------------
832 */
833 static int32_t
priv_timesplit(int32_t split[3],int32_t ts)834 priv_timesplit(
835 int32_t split[3],
836 int32_t ts
837 )
838 {
839 /* Do 3 chained floor divisions by positive constants, using the
840 * one's complement trick and factoring out the intermediate XOR
841 * ops to reduce the number of operations.
842 */
843 uint32_t us, um, uh, ud, sf32;
844
845 sf32 = int32_sflag(ts);
846
847 us = (uint32_t)ts;
848 um = (sf32 ^ us) / SECSPERMIN;
849 uh = um / MINSPERHR;
850 ud = uh / HRSPERDAY;
851
852 um ^= sf32;
853 uh ^= sf32;
854 ud ^= sf32;
855
856 split[0] = (int32_t)(uh - ud * HRSPERDAY );
857 split[1] = (int32_t)(um - uh * MINSPERHR );
858 split[2] = (int32_t)(us - um * SECSPERMIN);
859
860 return uint32_2cpl_to_int32(ud);
861 }
862
863 /*
864 *---------------------------------------------------------------------
865 * Given the number of elapsed days in the calendar era, split this
866 * number into the number of elapsed years in 'res.hi' and the number
867 * of elapsed days of that year in 'res.lo'.
868 *
869 * if 'isleapyear' is not NULL, it will receive an integer that is 0 for
870 * regular years and a non-zero value for leap years.
871 *---------------------------------------------------------------------
872 */
873 ntpcal_split
ntpcal_split_eradays(int32_t days,int * isleapyear)874 ntpcal_split_eradays(
875 int32_t days,
876 int *isleapyear
877 )
878 {
879 /* Use the fast cycle split algorithm here, to calculate the
880 * centuries and years in a century with one division each. This
881 * reduces the number of division operations to two, but is
882 * susceptible to internal range overflow. We take some extra
883 * steps to avoid the gap.
884 */
885 ntpcal_split res;
886 int32_t n100, n001; /* calendar year cycles */
887 uint32_t uday, Q;
888
889 /* split off centuries first
890 *
891 * We want to execute '(days * 4 + 3) /% 146097' under floor
892 * division rules in the first step. Well, actually we want to
893 * calculate 'floor((days + 0.75) / 36524.25)', but we want to
894 * do it in scaled integer calculation.
895 */
896 # if defined(HAVE_64BITREGS)
897
898 /* not too complicated with an intermediate 64bit value */
899 uint64_t ud64, sf64;
900 ud64 = ((uint64_t)days << 2) | 3u;
901 sf64 = (uint64_t)-(days < 0);
902 Q = (uint32_t)(sf64 ^ ((sf64 ^ ud64) / GREGORIAN_CYCLE_DAYS));
903 uday = (uint32_t)(ud64 - Q * GREGORIAN_CYCLE_DAYS);
904 n100 = uint32_2cpl_to_int32(Q);
905
906 # else
907
908 /* '4*days+3' suffers from range overflow when going to the
909 * limits. We solve this by doing an exact division (mod 2^32)
910 * after caclulating the remainder first.
911 *
912 * We start with a partial reduction by digit sums, extracting
913 * the upper bits from the original value before they get lost
914 * by scaling, and do one full division step to get the true
915 * remainder. Then a final multiplication with the
916 * multiplicative inverse of 146097 (mod 2^32) gives us the full
917 * quotient.
918 *
919 * (-2^33) % 146097 --> 130717 : the sign bit value
920 * ( 2^20) % 146097 --> 25897 : the upper digit value
921 * modinv(146097, 2^32) --> 660721233 : the inverse
922 */
923 uint32_t ux = ((uint32_t)days << 2) | 3;
924 uday = (days < 0) ? 130717u : 0u; /* sign dgt */
925 uday += ((days >> 18) & 0x01FFFu) * 25897u; /* hi dgt (src!) */
926 uday += (ux & 0xFFFFFu); /* lo dgt */
927 uday %= GREGORIAN_CYCLE_DAYS; /* full reduction */
928 Q = (ux - uday) * 660721233u; /* exact div */
929 n100 = uint32_2cpl_to_int32(Q);
930
931 # endif
932
933 /* Split off years in century -- days >= 0 here, and we're far
934 * away from integer overflow trouble now. */
935 uday |= 3;
936 n001 = uday / GREGORIAN_NORMAL_LEAP_CYCLE_DAYS;
937 uday -= n001 * GREGORIAN_NORMAL_LEAP_CYCLE_DAYS;
938
939 /* Assemble the year and day in year */
940 res.hi = n100 * 100 + n001;
941 res.lo = uday / 4u;
942
943 /* Possibly set the leap year flag */
944 if (isleapyear) {
945 uint32_t tc = (uint32_t)n100 + 1;
946 uint32_t ty = (uint32_t)n001 + 1;
947 *isleapyear = !(ty & 3)
948 && ((ty != 100) || !(tc & 3));
949 }
950 return res;
951 }
952
953 /*
954 *---------------------------------------------------------------------
955 * Given a number of elapsed days in a year and a leap year indicator,
956 * split the number of elapsed days into the number of elapsed months in
957 * 'res.hi' and the number of elapsed days of that month in 'res.lo'.
958 *
959 * This function will fail and return {-1,-1} if the number of elapsed
960 * days is not in the valid range!
961 *---------------------------------------------------------------------
962 */
963 ntpcal_split
ntpcal_split_yeardays(int32_t eyd,int isleap)964 ntpcal_split_yeardays(
965 int32_t eyd,
966 int isleap
967 )
968 {
969 /* Use the unshifted-year, February-with-30-days approach here.
970 * Fractional interpolations are used in both directions, with
971 * the smallest power-of-two divider to avoid any true division.
972 */
973 ntpcal_split res = {-1, -1};
974
975 /* convert 'isleap' to number of defective days */
976 isleap = 1 + !isleap;
977 /* adjust for February of 30 nominal days */
978 if (eyd >= 61 - isleap)
979 eyd += isleap;
980 /* if in range, convert to months and days in month */
981 if (eyd >= 0 && eyd < 367) {
982 res.hi = (eyd * 67 + 32) >> 11;
983 res.lo = eyd - ((489 * res.hi + 8) >> 4);
984 }
985
986 return res;
987 }
988
989 /*
990 *---------------------------------------------------------------------
991 * Convert a RD into the date part of a 'struct calendar'.
992 *---------------------------------------------------------------------
993 */
994 int
ntpcal_rd_to_date(struct calendar * jd,int32_t rd)995 ntpcal_rd_to_date(
996 struct calendar *jd,
997 int32_t rd
998 )
999 {
1000 ntpcal_split split;
1001 int leapy;
1002 u_int ymask;
1003
1004 /* Get day-of-week first. It's simply the RD (mod 7)... */
1005 jd->weekday = i32mod7(rd);
1006
1007 split = ntpcal_split_eradays(rd - 1, &leapy);
1008 /* Get year and day-of-year, with overflow check. If any of the
1009 * upper 16 bits is set after shifting to unity-based years, we
1010 * will have an overflow when converting to an unsigned 16bit
1011 * year. Shifting to the right is OK here, since it does not
1012 * matter if the shift is logic or arithmetic.
1013 */
1014 split.hi += 1;
1015 ymask = 0u - ((split.hi >> 16) == 0);
1016 jd->year = (uint16_t)(split.hi & ymask);
1017 jd->yearday = (uint16_t)split.lo + 1;
1018
1019 /* convert to month and mday */
1020 split = ntpcal_split_yeardays(split.lo, leapy);
1021 jd->month = (uint8_t)split.hi + 1;
1022 jd->monthday = (uint8_t)split.lo + 1;
1023
1024 return ymask ? leapy : -1;
1025 }
1026
1027 /*
1028 *---------------------------------------------------------------------
1029 * Convert a RD into the date part of a 'struct tm'.
1030 *---------------------------------------------------------------------
1031 */
1032 int
ntpcal_rd_to_tm(struct tm * utm,int32_t rd)1033 ntpcal_rd_to_tm(
1034 struct tm *utm,
1035 int32_t rd
1036 )
1037 {
1038 ntpcal_split split;
1039 int leapy;
1040
1041 /* get day-of-week first */
1042 utm->tm_wday = i32mod7(rd);
1043
1044 /* get year and day-of-year */
1045 split = ntpcal_split_eradays(rd - 1, &leapy);
1046 utm->tm_year = split.hi - 1899;
1047 utm->tm_yday = split.lo; /* 0-based */
1048
1049 /* convert to month and mday */
1050 split = ntpcal_split_yeardays(split.lo, leapy);
1051 utm->tm_mon = split.hi; /* 0-based */
1052 utm->tm_mday = split.lo + 1; /* 1-based */
1053
1054 return leapy;
1055 }
1056
1057 /*
1058 *---------------------------------------------------------------------
1059 * Take a value of seconds since midnight and split it into hhmmss in a
1060 * 'struct calendar'.
1061 *---------------------------------------------------------------------
1062 */
1063 int32_t
ntpcal_daysec_to_date(struct calendar * jd,int32_t sec)1064 ntpcal_daysec_to_date(
1065 struct calendar *jd,
1066 int32_t sec
1067 )
1068 {
1069 int32_t days;
1070 int ts[3];
1071
1072 days = priv_timesplit(ts, sec);
1073 jd->hour = (uint8_t)ts[0];
1074 jd->minute = (uint8_t)ts[1];
1075 jd->second = (uint8_t)ts[2];
1076
1077 return days;
1078 }
1079
1080 /*
1081 *---------------------------------------------------------------------
1082 * Take a value of seconds since midnight and split it into hhmmss in a
1083 * 'struct tm'.
1084 *---------------------------------------------------------------------
1085 */
1086 int32_t
ntpcal_daysec_to_tm(struct tm * utm,int32_t sec)1087 ntpcal_daysec_to_tm(
1088 struct tm *utm,
1089 int32_t sec
1090 )
1091 {
1092 int32_t days;
1093 int32_t ts[3];
1094
1095 days = priv_timesplit(ts, sec);
1096 utm->tm_hour = ts[0];
1097 utm->tm_min = ts[1];
1098 utm->tm_sec = ts[2];
1099
1100 return days;
1101 }
1102
1103 /*
1104 *---------------------------------------------------------------------
1105 * take a split representation for day/second-of-day and day offset
1106 * and convert it to a 'struct calendar'. The seconds will be normalised
1107 * into the range of a day, and the day will be adjusted accordingly.
1108 *
1109 * returns >0 if the result is in a leap year, 0 if in a regular
1110 * year and <0 if the result did not fit into the calendar struct.
1111 *---------------------------------------------------------------------
1112 */
1113 int
ntpcal_daysplit_to_date(struct calendar * jd,const ntpcal_split * ds,int32_t dof)1114 ntpcal_daysplit_to_date(
1115 struct calendar *jd,
1116 const ntpcal_split *ds,
1117 int32_t dof
1118 )
1119 {
1120 dof += ntpcal_daysec_to_date(jd, ds->lo);
1121 return ntpcal_rd_to_date(jd, ds->hi + dof);
1122 }
1123
1124 /*
1125 *---------------------------------------------------------------------
1126 * take a split representation for day/second-of-day and day offset
1127 * and convert it to a 'struct tm'. The seconds will be normalised
1128 * into the range of a day, and the day will be adjusted accordingly.
1129 *
1130 * returns 1 if the result is in a leap year and zero if in a regular
1131 * year.
1132 *---------------------------------------------------------------------
1133 */
1134 int
ntpcal_daysplit_to_tm(struct tm * utm,const ntpcal_split * ds,int32_t dof)1135 ntpcal_daysplit_to_tm(
1136 struct tm *utm,
1137 const ntpcal_split *ds ,
1138 int32_t dof
1139 )
1140 {
1141 dof += ntpcal_daysec_to_tm(utm, ds->lo);
1142
1143 return ntpcal_rd_to_tm(utm, ds->hi + dof);
1144 }
1145
1146 /*
1147 *---------------------------------------------------------------------
1148 * Take a UN*X time and convert to a calendar structure.
1149 *---------------------------------------------------------------------
1150 */
1151 int
ntpcal_time_to_date(struct calendar * jd,const vint64 * ts)1152 ntpcal_time_to_date(
1153 struct calendar *jd,
1154 const vint64 *ts
1155 )
1156 {
1157 ntpcal_split ds;
1158
1159 ds = ntpcal_daysplit(ts);
1160 ds.hi += ntpcal_daysec_to_date(jd, ds.lo);
1161 ds.hi += DAY_UNIX_STARTS;
1162
1163 return ntpcal_rd_to_date(jd, ds.hi);
1164 }
1165
1166
1167 /*
1168 * ====================================================================
1169 *
1170 * merging composite entities
1171 *
1172 * ====================================================================
1173 */
1174
1175 #if !defined(HAVE_INT64)
1176 /* multiplication helper. Seconds in days and weeks are multiples of 128,
1177 * and without that factor fit well into 16 bit. So a multiplication
1178 * of 32bit by 16bit and some shifting can be used on pure 32bit machines
1179 * with compilers that do not support 64bit integers.
1180 *
1181 * Calculate ( hi * mul * 128 ) + lo
1182 */
1183 static vint64
_dwjoin(uint16_t mul,int32_t hi,int32_t lo)1184 _dwjoin(
1185 uint16_t mul,
1186 int32_t hi,
1187 int32_t lo
1188 )
1189 {
1190 vint64 res;
1191 uint32_t p1, p2, sf;
1192
1193 /* get sign flag and absolute value of 'hi' in p1 */
1194 sf = (uint32_t)-(hi < 0);
1195 p1 = ((uint32_t)hi + sf) ^ sf;
1196
1197 /* assemble major units: res <- |hi| * mul */
1198 res.D_s.lo = (p1 & 0xFFFF) * mul;
1199 res.D_s.hi = 0;
1200 p1 = (p1 >> 16) * mul;
1201 p2 = p1 >> 16;
1202 p1 = p1 << 16;
1203 M_ADD(res.D_s.hi, res.D_s.lo, p2, p1);
1204
1205 /* mul by 128, using shift: res <-- res << 7 */
1206 res.D_s.hi = (res.D_s.hi << 7) | (res.D_s.lo >> 25);
1207 res.D_s.lo = (res.D_s.lo << 7);
1208
1209 /* fix up sign: res <-- (res + [sf|sf]) ^ [sf|sf] */
1210 M_ADD(res.D_s.hi, res.D_s.lo, sf, sf);
1211 res.D_s.lo ^= sf;
1212 res.D_s.hi ^= sf;
1213
1214 /* properly add seconds: res <-- res + [sx(lo)|lo] */
1215 p2 = (uint32_t)-(lo < 0);
1216 p1 = (uint32_t)lo;
1217 M_ADD(res.D_s.hi, res.D_s.lo, p2, p1);
1218 return res;
1219 }
1220 #endif
1221
1222 /*
1223 *---------------------------------------------------------------------
1224 * Merge a number of days and a number of seconds into seconds,
1225 * expressed in 64 bits to avoid overflow.
1226 *---------------------------------------------------------------------
1227 */
1228 vint64
ntpcal_dayjoin(int32_t days,int32_t secs)1229 ntpcal_dayjoin(
1230 int32_t days,
1231 int32_t secs
1232 )
1233 {
1234 vint64 res;
1235
1236 # if defined(HAVE_INT64)
1237
1238 res.q_s = days;
1239 res.q_s *= SECSPERDAY;
1240 res.q_s += secs;
1241
1242 # else
1243
1244 res = _dwjoin(675, days, secs);
1245
1246 # endif
1247
1248 return res;
1249 }
1250
1251 /*
1252 *---------------------------------------------------------------------
1253 * Merge a number of weeks and a number of seconds into seconds,
1254 * expressed in 64 bits to avoid overflow.
1255 *---------------------------------------------------------------------
1256 */
1257 vint64
ntpcal_weekjoin(int32_t week,int32_t secs)1258 ntpcal_weekjoin(
1259 int32_t week,
1260 int32_t secs
1261 )
1262 {
1263 vint64 res;
1264
1265 # if defined(HAVE_INT64)
1266
1267 res.q_s = week;
1268 res.q_s *= SECSPERWEEK;
1269 res.q_s += secs;
1270
1271 # else
1272
1273 res = _dwjoin(4725, week, secs);
1274
1275 # endif
1276
1277 return res;
1278 }
1279
1280 /*
1281 *---------------------------------------------------------------------
1282 * get leap years since epoch in elapsed years
1283 *---------------------------------------------------------------------
1284 */
1285 int32_t
ntpcal_leapyears_in_years(int32_t years)1286 ntpcal_leapyears_in_years(
1287 int32_t years
1288 )
1289 {
1290 /* We use the in-out-in algorithm here, using the one's
1291 * complement division trick for negative numbers. The chained
1292 * division sequence by 4/25/4 gives the compiler the chance to
1293 * get away with only one true division and doing shifts otherwise.
1294 */
1295
1296 uint32_t sf32, sum, uyear;
1297
1298 sf32 = int32_sflag(years);
1299 uyear = (uint32_t)years;
1300 uyear ^= sf32;
1301
1302 sum = (uyear /= 4u); /* 4yr rule --> IN */
1303 sum -= (uyear /= 25u); /* 100yr rule --> OUT */
1304 sum += (uyear /= 4u); /* 400yr rule --> IN */
1305
1306 /* Thanks to the alternation of IN/OUT/IN we can do the sum
1307 * directly and have a single one's complement operation
1308 * here. (Only if the years are negative, of course.) Otherwise
1309 * the one's complement would have to be done when
1310 * adding/subtracting the terms.
1311 */
1312 return uint32_2cpl_to_int32(sf32 ^ sum);
1313 }
1314
1315 /*
1316 *---------------------------------------------------------------------
1317 * Convert elapsed years in Era into elapsed days in Era.
1318 *---------------------------------------------------------------------
1319 */
1320 int32_t
ntpcal_days_in_years(int32_t years)1321 ntpcal_days_in_years(
1322 int32_t years
1323 )
1324 {
1325 return years * DAYSPERYEAR + ntpcal_leapyears_in_years(years);
1326 }
1327
1328 /*
1329 *---------------------------------------------------------------------
1330 * Convert a number of elapsed month in a year into elapsed days in year.
1331 *
1332 * The month will be normalized, and 'res.hi' will contain the
1333 * excessive years that must be considered when converting the years,
1334 * while 'res.lo' will contain the number of elapsed days since start
1335 * of the year.
1336 *
1337 * This code uses the shifted-month-approach to convert month to days,
1338 * because then there is no need to have explicit leap year
1339 * information. The slight disadvantage is that for most month values
1340 * the result is a negative value, and the year excess is one; the
1341 * conversion is then simply based on the start of the following year.
1342 *---------------------------------------------------------------------
1343 */
1344 ntpcal_split
ntpcal_days_in_months(int32_t m)1345 ntpcal_days_in_months(
1346 int32_t m
1347 )
1348 {
1349 ntpcal_split res;
1350
1351 /* Add ten months with proper year adjustment. */
1352 if (m < 2) {
1353 res.lo = m + 10;
1354 res.hi = 0;
1355 } else {
1356 res.lo = m - 2;
1357 res.hi = 1;
1358 }
1359
1360 /* Possibly normalise by floor division. This does not hapen for
1361 * input in normal range. */
1362 if (res.lo < 0 || res.lo >= 12) {
1363 uint32_t mu, Q, sf32;
1364 sf32 = int32_sflag(res.lo);
1365 mu = (uint32_t)res.lo;
1366 Q = sf32 ^ ((sf32 ^ mu) / 12u);
1367
1368 res.hi += uint32_2cpl_to_int32(Q);
1369 res.lo = mu - Q * 12u;
1370 }
1371
1372 /* Get cummulated days in year with unshift. Use the fractional
1373 * interpolation with smallest possible power of two in the
1374 * divider.
1375 */
1376 res.lo = ((res.lo * 979 + 16) >> 5) - 306;
1377
1378 return res;
1379 }
1380
1381 /*
1382 *---------------------------------------------------------------------
1383 * Convert ELAPSED years/months/days of gregorian calendar to elapsed
1384 * days in Gregorian epoch.
1385 *
1386 * If you want to convert years and days-of-year, just give a month of
1387 * zero.
1388 *---------------------------------------------------------------------
1389 */
1390 int32_t
ntpcal_edate_to_eradays(int32_t years,int32_t mons,int32_t mdays)1391 ntpcal_edate_to_eradays(
1392 int32_t years,
1393 int32_t mons,
1394 int32_t mdays
1395 )
1396 {
1397 ntpcal_split tmp;
1398 int32_t res;
1399
1400 if (mons) {
1401 tmp = ntpcal_days_in_months(mons);
1402 res = ntpcal_days_in_years(years + tmp.hi) + tmp.lo;
1403 } else
1404 res = ntpcal_days_in_years(years);
1405 res += mdays;
1406
1407 return res;
1408 }
1409
1410 /*
1411 *---------------------------------------------------------------------
1412 * Convert ELAPSED years/months/days of gregorian calendar to elapsed
1413 * days in year.
1414 *
1415 * Note: This will give the true difference to the start of the given
1416 * year, even if months & days are off-scale.
1417 *---------------------------------------------------------------------
1418 */
1419 int32_t
ntpcal_edate_to_yeardays(int32_t years,int32_t mons,int32_t mdays)1420 ntpcal_edate_to_yeardays(
1421 int32_t years,
1422 int32_t mons,
1423 int32_t mdays
1424 )
1425 {
1426 ntpcal_split tmp;
1427
1428 if (0 <= mons && mons < 12) {
1429 if (mons >= 2)
1430 mdays -= 2 - is_leapyear(years+1);
1431 mdays += (489 * mons + 8) >> 4;
1432 } else {
1433 tmp = ntpcal_days_in_months(mons);
1434 mdays += tmp.lo
1435 + ntpcal_days_in_years(years + tmp.hi)
1436 - ntpcal_days_in_years(years);
1437 }
1438
1439 return mdays;
1440 }
1441
1442 /*
1443 *---------------------------------------------------------------------
1444 * Convert elapsed days and the hour/minute/second information into
1445 * total seconds.
1446 *
1447 * If 'isvalid' is not NULL, do a range check on the time specification
1448 * and tell if the time input is in the normal range, permitting for a
1449 * single leapsecond.
1450 *---------------------------------------------------------------------
1451 */
1452 int32_t
ntpcal_etime_to_seconds(int32_t hours,int32_t minutes,int32_t seconds)1453 ntpcal_etime_to_seconds(
1454 int32_t hours,
1455 int32_t minutes,
1456 int32_t seconds
1457 )
1458 {
1459 int32_t res;
1460
1461 res = (hours * MINSPERHR + minutes) * SECSPERMIN + seconds;
1462
1463 return res;
1464 }
1465
1466 /*
1467 *---------------------------------------------------------------------
1468 * Convert the date part of a 'struct tm' (that is, year, month,
1469 * day-of-month) into the RD of that day.
1470 *---------------------------------------------------------------------
1471 */
1472 int32_t
ntpcal_tm_to_rd(const struct tm * utm)1473 ntpcal_tm_to_rd(
1474 const struct tm *utm
1475 )
1476 {
1477 return ntpcal_edate_to_eradays(utm->tm_year + 1899,
1478 utm->tm_mon,
1479 utm->tm_mday - 1) + 1;
1480 }
1481
1482 /*
1483 *---------------------------------------------------------------------
1484 * Convert the date part of a 'struct calendar' (that is, year, month,
1485 * day-of-month) into the RD of that day.
1486 *---------------------------------------------------------------------
1487 */
1488 int32_t
ntpcal_date_to_rd(const struct calendar * jd)1489 ntpcal_date_to_rd(
1490 const struct calendar *jd
1491 )
1492 {
1493 return ntpcal_edate_to_eradays((int32_t)jd->year - 1,
1494 (int32_t)jd->month - 1,
1495 (int32_t)jd->monthday - 1) + 1;
1496 }
1497
1498 /*
1499 *---------------------------------------------------------------------
1500 * convert a year number to rata die of year start
1501 *---------------------------------------------------------------------
1502 */
1503 int32_t
ntpcal_year_to_ystart(int32_t year)1504 ntpcal_year_to_ystart(
1505 int32_t year
1506 )
1507 {
1508 return ntpcal_days_in_years(year - 1) + 1;
1509 }
1510
1511 /*
1512 *---------------------------------------------------------------------
1513 * For a given RD, get the RD of the associated year start,
1514 * that is, the RD of the last January,1st on or before that day.
1515 *---------------------------------------------------------------------
1516 */
1517 int32_t
ntpcal_rd_to_ystart(int32_t rd)1518 ntpcal_rd_to_ystart(
1519 int32_t rd
1520 )
1521 {
1522 /*
1523 * Rather simple exercise: split the day number into elapsed
1524 * years and elapsed days, then remove the elapsed days from the
1525 * input value. Nice'n sweet...
1526 */
1527 return rd - ntpcal_split_eradays(rd - 1, NULL).lo;
1528 }
1529
1530 /*
1531 *---------------------------------------------------------------------
1532 * For a given RD, get the RD of the associated month start.
1533 *---------------------------------------------------------------------
1534 */
1535 int32_t
ntpcal_rd_to_mstart(int32_t rd)1536 ntpcal_rd_to_mstart(
1537 int32_t rd
1538 )
1539 {
1540 ntpcal_split split;
1541 int leaps;
1542
1543 split = ntpcal_split_eradays(rd - 1, &leaps);
1544 split = ntpcal_split_yeardays(split.lo, leaps);
1545
1546 return rd - split.lo;
1547 }
1548
1549 /*
1550 *---------------------------------------------------------------------
1551 * take a 'struct calendar' and get the seconds-of-day from it.
1552 *---------------------------------------------------------------------
1553 */
1554 int32_t
ntpcal_date_to_daysec(const struct calendar * jd)1555 ntpcal_date_to_daysec(
1556 const struct calendar *jd
1557 )
1558 {
1559 return ntpcal_etime_to_seconds(jd->hour, jd->minute,
1560 jd->second);
1561 }
1562
1563 /*
1564 *---------------------------------------------------------------------
1565 * take a 'struct tm' and get the seconds-of-day from it.
1566 *---------------------------------------------------------------------
1567 */
1568 int32_t
ntpcal_tm_to_daysec(const struct tm * utm)1569 ntpcal_tm_to_daysec(
1570 const struct tm *utm
1571 )
1572 {
1573 return ntpcal_etime_to_seconds(utm->tm_hour, utm->tm_min,
1574 utm->tm_sec);
1575 }
1576
1577 /*
1578 *---------------------------------------------------------------------
1579 * take a 'struct calendar' and convert it to a 'time_t'
1580 *---------------------------------------------------------------------
1581 */
1582 time_t
ntpcal_date_to_time(const struct calendar * jd)1583 ntpcal_date_to_time(
1584 const struct calendar *jd
1585 )
1586 {
1587 vint64 join;
1588 int32_t days, secs;
1589
1590 days = ntpcal_date_to_rd(jd) - DAY_UNIX_STARTS;
1591 secs = ntpcal_date_to_daysec(jd);
1592 join = ntpcal_dayjoin(days, secs);
1593
1594 return vint64_to_time(&join);
1595 }
1596
1597
1598 /*
1599 * ====================================================================
1600 *
1601 * extended and unchecked variants of caljulian/caltontp
1602 *
1603 * ====================================================================
1604 */
1605 int
ntpcal_ntp64_to_date(struct calendar * jd,const vint64 * ntp)1606 ntpcal_ntp64_to_date(
1607 struct calendar *jd,
1608 const vint64 *ntp
1609 )
1610 {
1611 ntpcal_split ds;
1612
1613 ds = ntpcal_daysplit(ntp);
1614 ds.hi += ntpcal_daysec_to_date(jd, ds.lo);
1615
1616 return ntpcal_rd_to_date(jd, ds.hi + DAY_NTP_STARTS);
1617 }
1618
1619 int
ntpcal_ntp_to_date(struct calendar * jd,uint32_t ntp,const time_t * piv)1620 ntpcal_ntp_to_date(
1621 struct calendar *jd,
1622 uint32_t ntp,
1623 const time_t *piv
1624 )
1625 {
1626 vint64 ntp64;
1627
1628 /*
1629 * Unfold ntp time around current time into NTP domain. Split
1630 * into days and seconds, shift days into CE domain and
1631 * process the parts.
1632 */
1633 ntp64 = ntpcal_ntp_to_ntp(ntp, piv);
1634 return ntpcal_ntp64_to_date(jd, &ntp64);
1635 }
1636
1637
1638 vint64
ntpcal_date_to_ntp64(const struct calendar * jd)1639 ntpcal_date_to_ntp64(
1640 const struct calendar *jd
1641 )
1642 {
1643 /*
1644 * Convert date to NTP. Ignore yearday, use d/m/y only.
1645 */
1646 return ntpcal_dayjoin(ntpcal_date_to_rd(jd) - DAY_NTP_STARTS,
1647 ntpcal_date_to_daysec(jd));
1648 }
1649
1650
1651 uint32_t
ntpcal_date_to_ntp(const struct calendar * jd)1652 ntpcal_date_to_ntp(
1653 const struct calendar *jd
1654 )
1655 {
1656 /*
1657 * Get lower half of 64bit NTP timestamp from date/time.
1658 */
1659 return ntpcal_date_to_ntp64(jd).d_s.lo;
1660 }
1661
1662
1663
1664 /*
1665 * ====================================================================
1666 *
1667 * day-of-week calculations
1668 *
1669 * ====================================================================
1670 */
1671 /*
1672 * Given a RataDie and a day-of-week, calculate a RDN that is reater-than,
1673 * greater-or equal, closest, less-or-equal or less-than the given RDN
1674 * and denotes the given day-of-week
1675 */
1676 int32_t
ntpcal_weekday_gt(int32_t rdn,int32_t dow)1677 ntpcal_weekday_gt(
1678 int32_t rdn,
1679 int32_t dow
1680 )
1681 {
1682 return ntpcal_periodic_extend(rdn+1, dow, 7);
1683 }
1684
1685 int32_t
ntpcal_weekday_ge(int32_t rdn,int32_t dow)1686 ntpcal_weekday_ge(
1687 int32_t rdn,
1688 int32_t dow
1689 )
1690 {
1691 return ntpcal_periodic_extend(rdn, dow, 7);
1692 }
1693
1694 int32_t
ntpcal_weekday_close(int32_t rdn,int32_t dow)1695 ntpcal_weekday_close(
1696 int32_t rdn,
1697 int32_t dow
1698 )
1699 {
1700 return ntpcal_periodic_extend(rdn-3, dow, 7);
1701 }
1702
1703 int32_t
ntpcal_weekday_le(int32_t rdn,int32_t dow)1704 ntpcal_weekday_le(
1705 int32_t rdn,
1706 int32_t dow
1707 )
1708 {
1709 return ntpcal_periodic_extend(rdn, dow, -7);
1710 }
1711
1712 int32_t
ntpcal_weekday_lt(int32_t rdn,int32_t dow)1713 ntpcal_weekday_lt(
1714 int32_t rdn,
1715 int32_t dow
1716 )
1717 {
1718 return ntpcal_periodic_extend(rdn-1, dow, -7);
1719 }
1720
1721 /*
1722 * ====================================================================
1723 *
1724 * ISO week-calendar conversions
1725 *
1726 * The ISO8601 calendar defines a calendar of years, weeks and weekdays.
1727 * It is related to the Gregorian calendar, and a ISO year starts at the
1728 * Monday closest to Jan,1st of the corresponding Gregorian year. A ISO
1729 * calendar year has always 52 or 53 weeks, and like the Grogrian
1730 * calendar the ISO8601 calendar repeats itself every 400 years, or
1731 * 146097 days, or 20871 weeks.
1732 *
1733 * While it is possible to write ISO calendar functions based on the
1734 * Gregorian calendar functions, the following implementation takes a
1735 * different approach, based directly on years and weeks.
1736 *
1737 * Analysis of the tabulated data shows that it is not possible to
1738 * interpolate from years to weeks over a full 400 year range; cyclic
1739 * shifts over 400 years do not provide a solution here. But it *is*
1740 * possible to interpolate over every single century of the 400-year
1741 * cycle. (The centennial leap year rule seems to be the culprit here.)
1742 *
1743 * It can be shown that a conversion from years to weeks can be done
1744 * using a linear transformation of the form
1745 *
1746 * w = floor( y * a + b )
1747 *
1748 * where the slope a must hold to
1749 *
1750 * 52.1780821918 <= a < 52.1791044776
1751 *
1752 * and b must be chosen according to the selected slope and the number
1753 * of the century in a 400-year period.
1754 *
1755 * The inverse calculation can also be done in this way. Careful scaling
1756 * provides an unlimited set of integer coefficients a,k,b that enable
1757 * us to write the calulation in the form
1758 *
1759 * w = (y * a + b ) / k
1760 * y = (w * a' + b') / k'
1761 *
1762 * In this implementation the values of k and k' are chosen to be the
1763 * smallest possible powers of two, so the division can be implemented
1764 * as shifts if the optimiser chooses to do so.
1765 *
1766 * ====================================================================
1767 */
1768
1769 /*
1770 * Given a number of elapsed (ISO-)years since the begin of the
1771 * christian era, return the number of elapsed weeks corresponding to
1772 * the number of years.
1773 */
1774 int32_t
isocal_weeks_in_years(int32_t years)1775 isocal_weeks_in_years(
1776 int32_t years
1777 )
1778 {
1779 /*
1780 * use: w = (y * 53431 + b[c]) / 1024 as interpolation
1781 */
1782 static const uint16_t bctab[4] = { 157, 449, 597, 889 };
1783
1784 int32_t cs, cw;
1785 uint32_t cc, ci, yu, sf32;
1786
1787 sf32 = int32_sflag(years);
1788 yu = (uint32_t)years;
1789
1790 /* split off centuries, using floor division */
1791 cc = sf32 ^ ((sf32 ^ yu) / 100u);
1792 yu -= cc * 100u;
1793
1794 /* calculate century cycles shift and cycle index:
1795 * Assuming a century is 5217 weeks, we have to add a cycle
1796 * shift that is 3 for every 4 centuries, because 3 of the four
1797 * centuries have 5218 weeks. So '(cc*3 + 1) / 4' is the actual
1798 * correction, and the second century is the defective one.
1799 *
1800 * Needs floor division by 4, which is done with masking and
1801 * shifting.
1802 */
1803 ci = cc * 3u + 1;
1804 cs = uint32_2cpl_to_int32(sf32 ^ ((sf32 ^ ci) >> 2));
1805 ci = ci & 3u;
1806
1807 /* Get weeks in century. Can use plain division here as all ops
1808 * are >= 0, and let the compiler sort out the possible
1809 * optimisations.
1810 */
1811 cw = (yu * 53431u + bctab[ci]) / 1024u;
1812
1813 return uint32_2cpl_to_int32(cc) * 5217 + cs + cw;
1814 }
1815
1816 /*
1817 * Given a number of elapsed weeks since the begin of the christian
1818 * era, split this number into the number of elapsed years in res.hi
1819 * and the excessive number of weeks in res.lo. (That is, res.lo is
1820 * the number of elapsed weeks in the remaining partial year.)
1821 */
1822 ntpcal_split
isocal_split_eraweeks(int32_t weeks)1823 isocal_split_eraweeks(
1824 int32_t weeks
1825 )
1826 {
1827 /*
1828 * use: y = (w * 157 + b[c]) / 8192 as interpolation
1829 */
1830
1831 static const uint16_t bctab[4] = { 85, 130, 17, 62 };
1832
1833 ntpcal_split res;
1834 int32_t cc, ci;
1835 uint32_t sw, cy, Q;
1836
1837 /* Use two fast cycle-split divisions again. Herew e want to
1838 * execute '(weeks * 4 + 2) /% 20871' under floor division rules
1839 * in the first step.
1840 *
1841 * This is of course (again) susceptible to internal overflow if
1842 * coded directly in 32bit. And again we use 64bit division on
1843 * a 64bit target and exact division after calculating the
1844 * remainder first on a 32bit target. With the smaller divider,
1845 * that's even a bit neater.
1846 */
1847 # if defined(HAVE_64BITREGS)
1848
1849 /* Full floor division with 64bit values. */
1850 uint64_t sf64, sw64;
1851 sf64 = (uint64_t)-(weeks < 0);
1852 sw64 = ((uint64_t)weeks << 2) | 2u;
1853 Q = (uint32_t)(sf64 ^ ((sf64 ^ sw64) / GREGORIAN_CYCLE_WEEKS));
1854 sw = (uint32_t)(sw64 - Q * GREGORIAN_CYCLE_WEEKS);
1855
1856 # else
1857
1858 /* Exact division after calculating the remainder via partial
1859 * reduction by digit sum.
1860 * (-2^33) % 20871 --> 5491 : the sign bit value
1861 * ( 2^20) % 20871 --> 5026 : the upper digit value
1862 * modinv(20871, 2^32) --> 330081335 : the inverse
1863 */
1864 uint32_t ux = ((uint32_t)weeks << 2) | 2;
1865 sw = (weeks < 0) ? 5491u : 0u; /* sign dgt */
1866 sw += ((weeks >> 18) & 0x01FFFu) * 5026u; /* hi dgt (src!) */
1867 sw += (ux & 0xFFFFFu); /* lo dgt */
1868 sw %= GREGORIAN_CYCLE_WEEKS; /* full reduction */
1869 Q = (ux - sw) * 330081335u; /* exact div */
1870
1871 # endif
1872
1873 ci = Q & 3u;
1874 cc = uint32_2cpl_to_int32(Q);
1875
1876 /* Split off years; sw >= 0 here! The scaled weeks in the years
1877 * are scaled up by 157 afterwards.
1878 */
1879 sw = (sw / 4u) * 157u + bctab[ci];
1880 cy = sw / 8192u; /* sw >> 13 , let the compiler sort it out */
1881 sw = sw % 8192u; /* sw & 8191, let the compiler sort it out */
1882
1883 /* assemble elapsed years and downscale the elapsed weeks in
1884 * the year.
1885 */
1886 res.hi = 100*cc + cy;
1887 res.lo = sw / 157u;
1888
1889 return res;
1890 }
1891
1892 /*
1893 * Given a second in the NTP time scale and a pivot, expand the NTP
1894 * time stamp around the pivot and convert into an ISO calendar time
1895 * stamp.
1896 */
1897 int
isocal_ntp64_to_date(struct isodate * id,const vint64 * ntp)1898 isocal_ntp64_to_date(
1899 struct isodate *id,
1900 const vint64 *ntp
1901 )
1902 {
1903 ntpcal_split ds;
1904 int32_t ts[3];
1905 uint32_t uw, ud, sf32;
1906
1907 /*
1908 * Split NTP time into days and seconds, shift days into CE
1909 * domain and process the parts.
1910 */
1911 ds = ntpcal_daysplit(ntp);
1912
1913 /* split time part */
1914 ds.hi += priv_timesplit(ts, ds.lo);
1915 id->hour = (uint8_t)ts[0];
1916 id->minute = (uint8_t)ts[1];
1917 id->second = (uint8_t)ts[2];
1918
1919 /* split days into days and weeks, using floor division in unsigned */
1920 ds.hi += DAY_NTP_STARTS - 1; /* shift from NTP to RDN */
1921 sf32 = int32_sflag(ds.hi);
1922 ud = (uint32_t)ds.hi;
1923 uw = sf32 ^ ((sf32 ^ ud) / DAYSPERWEEK);
1924 ud -= uw * DAYSPERWEEK;
1925
1926 ds.hi = uint32_2cpl_to_int32(uw);
1927 ds.lo = ud;
1928
1929 id->weekday = (uint8_t)ds.lo + 1; /* weekday result */
1930
1931 /* get year and week in year */
1932 ds = isocal_split_eraweeks(ds.hi); /* elapsed years&week*/
1933 id->year = (uint16_t)ds.hi + 1; /* shift to current */
1934 id->week = (uint8_t )ds.lo + 1;
1935
1936 return (ds.hi >= 0 && ds.hi < 0x0000FFFF);
1937 }
1938
1939 int
isocal_ntp_to_date(struct isodate * id,uint32_t ntp,const time_t * piv)1940 isocal_ntp_to_date(
1941 struct isodate *id,
1942 uint32_t ntp,
1943 const time_t *piv
1944 )
1945 {
1946 vint64 ntp64;
1947
1948 /*
1949 * Unfold ntp time around current time into NTP domain, then
1950 * convert the full time stamp.
1951 */
1952 ntp64 = ntpcal_ntp_to_ntp(ntp, piv);
1953 return isocal_ntp64_to_date(id, &ntp64);
1954 }
1955
1956 /*
1957 * Convert a ISO date spec into a second in the NTP time scale,
1958 * properly truncated to 32 bit.
1959 */
1960 vint64
isocal_date_to_ntp64(const struct isodate * id)1961 isocal_date_to_ntp64(
1962 const struct isodate *id
1963 )
1964 {
1965 int32_t weeks, days, secs;
1966
1967 weeks = isocal_weeks_in_years((int32_t)id->year - 1)
1968 + (int32_t)id->week - 1;
1969 days = weeks * 7 + (int32_t)id->weekday;
1970 /* days is RDN of ISO date now */
1971 secs = ntpcal_etime_to_seconds(id->hour, id->minute, id->second);
1972
1973 return ntpcal_dayjoin(days - DAY_NTP_STARTS, secs);
1974 }
1975
1976 uint32_t
isocal_date_to_ntp(const struct isodate * id)1977 isocal_date_to_ntp(
1978 const struct isodate *id
1979 )
1980 {
1981 /*
1982 * Get lower half of 64bit NTP timestamp from date/time.
1983 */
1984 return isocal_date_to_ntp64(id).d_s.lo;
1985 }
1986
1987 /*
1988 * ====================================================================
1989 * 'basedate' support functions
1990 * ====================================================================
1991 */
1992
1993 static int32_t s_baseday = NTP_TO_UNIX_DAYS;
1994 static int32_t s_gpsweek = 0;
1995
1996 int32_t
basedate_eval_buildstamp(void)1997 basedate_eval_buildstamp(void)
1998 {
1999 struct calendar jd;
2000 int32_t ed;
2001
2002 if (!ntpcal_get_build_date(&jd))
2003 return NTP_TO_UNIX_DAYS;
2004
2005 /* The time zone of the build stamp is unspecified; we remove
2006 * one day to provide a certain slack. And in case somebody
2007 * fiddled with the system clock, we make sure we do not go
2008 * before the UNIX epoch (1970-01-01). It's probably not possible
2009 * to do this to the clock on most systems, but there are other
2010 * ways to tweak the build stamp.
2011 */
2012 jd.monthday -= 1;
2013 ed = ntpcal_date_to_rd(&jd) - DAY_NTP_STARTS;
2014 return (ed < NTP_TO_UNIX_DAYS) ? NTP_TO_UNIX_DAYS : ed;
2015 }
2016
2017 int32_t
basedate_eval_string(const char * str)2018 basedate_eval_string(
2019 const char * str
2020 )
2021 {
2022 u_short y,m,d;
2023 u_long ned;
2024 int rc, nc;
2025 size_t sl;
2026
2027 sl = strlen(str);
2028 rc = sscanf(str, "%4hu-%2hu-%2hu%n", &y, &m, &d, &nc);
2029 if (rc == 3 && (size_t)nc == sl) {
2030 if (m >= 1 && m <= 12 && d >= 1 && d <= 31)
2031 return ntpcal_edate_to_eradays(y-1, m-1, d)
2032 - DAY_NTP_STARTS;
2033 goto buildstamp;
2034 }
2035
2036 rc = sscanf(str, "%lu%n", &ned, &nc);
2037 if (rc == 1 && (size_t)nc == sl) {
2038 if (ned <= INT32_MAX)
2039 return (int32_t)ned;
2040 goto buildstamp;
2041 }
2042
2043 buildstamp:
2044 msyslog(LOG_WARNING,
2045 "basedate string \"%s\" invalid, build date substituted!",
2046 str);
2047 return basedate_eval_buildstamp();
2048 }
2049
2050 uint32_t
basedate_get_day(void)2051 basedate_get_day(void)
2052 {
2053 return s_baseday;
2054 }
2055
2056 int32_t
basedate_set_day(int32_t day)2057 basedate_set_day(
2058 int32_t day
2059 )
2060 {
2061 struct calendar jd;
2062 int32_t retv;
2063
2064 /* set NTP base date for NTP era unfolding */
2065 if (day < NTP_TO_UNIX_DAYS) {
2066 msyslog(LOG_WARNING,
2067 "baseday_set_day: invalid day (%lu), UNIX epoch substituted",
2068 (unsigned long)day);
2069 day = NTP_TO_UNIX_DAYS;
2070 }
2071 retv = s_baseday;
2072 s_baseday = day;
2073 ntpcal_rd_to_date(&jd, day + DAY_NTP_STARTS);
2074 msyslog(LOG_INFO, "basedate set to %04hu-%02hu-%02hu",
2075 jd.year, (u_short)jd.month, (u_short)jd.monthday);
2076
2077 /* set GPS base week for GPS week unfolding */
2078 day = ntpcal_weekday_ge(day + DAY_NTP_STARTS, CAL_SUNDAY)
2079 - DAY_NTP_STARTS;
2080 if (day < NTP_TO_GPS_DAYS)
2081 day = NTP_TO_GPS_DAYS;
2082 s_gpsweek = (day - NTP_TO_GPS_DAYS) / DAYSPERWEEK;
2083 ntpcal_rd_to_date(&jd, day + DAY_NTP_STARTS);
2084 msyslog(LOG_INFO, "gps base set to %04hu-%02hu-%02hu (week %d)",
2085 jd.year, (u_short)jd.month, (u_short)jd.monthday, s_gpsweek);
2086
2087 return retv;
2088 }
2089
2090 time_t
basedate_get_eracenter(void)2091 basedate_get_eracenter(void)
2092 {
2093 time_t retv;
2094 retv = (time_t)(s_baseday - NTP_TO_UNIX_DAYS);
2095 retv *= SECSPERDAY;
2096 retv += (UINT32_C(1) << 31);
2097 return retv;
2098 }
2099
2100 time_t
basedate_get_erabase(void)2101 basedate_get_erabase(void)
2102 {
2103 time_t retv;
2104 retv = (time_t)(s_baseday - NTP_TO_UNIX_DAYS);
2105 retv *= SECSPERDAY;
2106 return retv;
2107 }
2108
2109 uint32_t
basedate_get_gpsweek(void)2110 basedate_get_gpsweek(void)
2111 {
2112 return s_gpsweek;
2113 }
2114
2115 uint32_t
basedate_expand_gpsweek(unsigned short weekno)2116 basedate_expand_gpsweek(
2117 unsigned short weekno
2118 )
2119 {
2120 /* We do a fast modulus expansion here. Since all quantities are
2121 * unsigned and we cannot go before the start of the GPS epoch
2122 * anyway, and since the truncated GPS week number is 10 bit, the
2123 * expansion becomes a simple sub/and/add sequence.
2124 */
2125 #if GPSWEEKS != 1024
2126 # error GPSWEEKS defined wrong -- should be 1024!
2127 #endif
2128
2129 uint32_t diff;
2130 diff = ((uint32_t)weekno - s_gpsweek) & (GPSWEEKS - 1);
2131 return s_gpsweek + diff;
2132 }
2133
2134 /*
2135 * ====================================================================
2136 * misc. helpers
2137 * ====================================================================
2138 */
2139
2140 /* --------------------------------------------------------------------
2141 * reconstruct the centrury from a truncated date and a day-of-week
2142 *
2143 * Given a date with truncated year (2-digit, 0..99) and a day-of-week
2144 * from 1(Mon) to 7(Sun), recover the full year between 1900AD and 2300AD.
2145 */
2146 int32_t
ntpcal_expand_century(uint32_t y,uint32_t m,uint32_t d,uint32_t wd)2147 ntpcal_expand_century(
2148 uint32_t y,
2149 uint32_t m,
2150 uint32_t d,
2151 uint32_t wd)
2152 {
2153 /* This algorithm is short but tricky... It's related to
2154 * Zeller's congruence, partially done backwards.
2155 *
2156 * A few facts to remember:
2157 * 1) The Gregorian calendar has a cycle of 400 years.
2158 * 2) The weekday of the 1st day of a century shifts by 5 days
2159 * during a great cycle.
2160 * 3) For calendar math, a century starts with the 1st year,
2161 * which is year 1, !not! zero.
2162 *
2163 * So we start with taking the weekday difference (mod 7)
2164 * between the truncated date (which is taken as an absolute
2165 * date in the 1st century in the proleptic calendar) and the
2166 * weekday given.
2167 *
2168 * When dividing this residual by 5, we obtain the number of
2169 * centuries to add to the base. But since the residual is (mod
2170 * 7), we have to make this an exact division by multiplication
2171 * with the modular inverse of 5 (mod 7), which is 3:
2172 * 3*5 === 1 (mod 7).
2173 *
2174 * If this yields a result of 4/5/6, the given date/day-of-week
2175 * combination is impossible, and we return zero as resulting
2176 * year to indicate failure.
2177 *
2178 * Then we remap the century to the range starting with year
2179 * 1900.
2180 */
2181
2182 uint32_t c;
2183
2184 /* check basic constraints */
2185 if ((y >= 100u) || (--m >= 12u) || (--d >= 31u))
2186 return 0;
2187
2188 if ((m += 10u) >= 12u) /* shift base to prev. March,1st */
2189 m -= 12u;
2190 else if (--y >= 100u)
2191 y += 100u;
2192 d += y + (y >> 2) + 2u; /* year share */
2193 d += (m * 83u + 16u) >> 5; /* month share */
2194
2195 /* get (wd - d), shifted to positive value, and multiply with
2196 * 3(mod 7). (Exact division, see to comment)
2197 * Note: 1) d <= 184 at this point.
2198 * 2) 252 % 7 == 0, but 'wd' is off by one since we did
2199 * '--d' above, so we add just 251 here!
2200 */
2201 c = u32mod7(3 * (251u + wd - d));
2202 if (c > 3u)
2203 return 0;
2204
2205 if ((m > 9u) && (++y >= 100u)) {/* undo base shift */
2206 y -= 100u;
2207 c = (c + 1) & 3u;
2208 }
2209 y += (c * 100u); /* combine into 1st cycle */
2210 y += (y < 300u) ? 2000 : 1600; /* map to destination era */
2211 return (int)y;
2212 }
2213
2214 char *
ntpcal_iso8601std(char * buf,size_t len,TcCivilDate * cdp)2215 ntpcal_iso8601std(
2216 char * buf,
2217 size_t len,
2218 TcCivilDate * cdp
2219 )
2220 {
2221 if (!buf) {
2222 LIB_GETBUF(buf);
2223 len = LIB_BUFLENGTH;
2224 }
2225 if (len) {
2226 len = snprintf(buf, len, "%04u-%02u-%02uT%02u:%02u:%02u",
2227 cdp->year, cdp->month, cdp->monthday,
2228 cdp->hour, cdp->minute, cdp->second);
2229 if (len < 0)
2230 *buf = '\0';
2231 }
2232 return buf;
2233 }
2234
2235 /* -*-EOF-*- */
2236