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