xref: /netbsd/src/external/bsd/ntp/dist/libntp/ntp_calendar.c

/*        $NetBSD: ntp_calendar.c,v 1.12 2024/08/18 20:47:13 christos Exp $     */

/*
 * ntp_calendar.c - calendar and helper functions
 *
 * Written by Juergen Perlinger (perlinger@ntp.org) for the NTP project.
 * The contents of 'html/copyright.html' apply.
 *
 * --------------------------------------------------------------------
 * Some notes on the implementation:
 *
 * Calendar algorithms thrive on the division operation, which is one of
 * the slowest numerical operations in any CPU. What saves us here from
 * abysmal performance is the fact that all divisions are divisions by
 * constant numbers, and most compilers can do this by a multiplication
 * operation.  But this might not work when using the div/ldiv/lldiv
 * function family, because many compilers are not able to do inline
 * expansion of the code with following optimisation for the
 * constant-divider case.
 *
 * Also div/ldiv/lldiv are defined in terms of int/long/longlong, which
 * are inherently target dependent. Nothing that could not be cured with
 * autoconf, but still a mess...
 *
 * Furthermore, we need floor division in many places. C either leaves
 * the division behaviour undefined (< C99) or demands truncation to
 * zero (>= C99), so additional steps are required to make sure the
 * algorithms work. The {l,ll}div function family is requested to
 * truncate towards zero, which is also the wrong direction for our
 * purpose.
 *
 * For all this, all divisions by constant are coded manually, even when
 * there is a joined div/mod operation: The optimiser should sort that
 * out, if possible. Most of the calculations are done with unsigned
 * types, explicitely using two's complement arithmetics where
 * necessary. This minimises the dependecies to compiler and target,
 * while still giving reasonable to good performance.
 *
 * The implementation uses a few tricks that exploit properties of the
 * two's complement: Floor division on negative dividents can be
 * executed by using the one's complement of the divident. One's
 * complement can be easily created using XOR and a mask.
 *
 * Finally, check for overflow conditions is minimal. There are only two
 * calculation steps in the whole calendar that potentially suffer from
 * an internal overflow, and these are coded in a way that avoids
 * it. All other functions do not suffer from internal overflow and
 * simply return the result truncated to 32 bits.
 */

#include <config.h>
#include <sys/types.h>

#include "ntp_types.h"
#include "ntp_calendar.h"
#include "ntp_stdlib.h"
#include "ntp_fp.h"
#include "ntp_unixtime.h"

#include "ntpd.h"

/* For now, let's take the conservative approach: if the target property
 * macros are not defined, check a few well-known compiler/architecture
 * settings. Default is to assume that the representation of signed
 * integers is unknown and shift-arithmetic-right is not available.
 */
#ifndef TARGET_HAS_2CPL
# if defined(__GNUC__)
#  if defined(__i386__) || defined(__x86_64__) || defined(__arm__)
#   define TARGET_HAS_2CPL 1
#  else
#   define TARGET_HAS_2CPL 0
#  endif
# elif defined(_MSC_VER)
#  if defined(_M_IX86) || defined(_M_X64) || defined(_M_ARM)
#   define TARGET_HAS_2CPL 1
#  else
#   define TARGET_HAS_2CPL 0
#  endif
# else
#  define TARGET_HAS_2CPL 0
# endif
#endif

#ifndef TARGET_HAS_SAR
# define TARGET_HAS_SAR 0
#endif

#if !defined(HAVE_64BITREGS) && defined(UINT64_MAX) && (SIZE_MAX >= UINT64_MAX)
# define HAVE_64BITREGS
#endif

/*
 *---------------------------------------------------------------------
 * replacing the 'time()' function
 *---------------------------------------------------------------------
 */

static systime_func_ptr systime_func = &time;
static inline time_t now(void);


systime_func_ptr
ntpcal_set_timefunc(
          systime_func_ptr nfunc
          )
{
          systime_func_ptr res;

          res = systime_func;
          if (NULL == nfunc)
                    nfunc = &time;
          systime_func = nfunc;

          return res;
}


static inline time_t
now(void)
{
          return (*systime_func)(NULL);
}

/*
 *---------------------------------------------------------------------
 * Get sign extension mask and unsigned 2cpl rep for a signed integer
 *---------------------------------------------------------------------
 */

static inline uint32_t
int32_sflag(
          const int32_t v)
{
#   if TARGET_HAS_2CPL && TARGET_HAS_SAR && SIZEOF_INT >= 4

          /* Let's assume that shift is the fastest way to get the sign
           * extension of of a signed integer. This might not always be
           * true, though -- On 8bit CPUs or machines without barrel
           * shifter this will kill the performance. So we make sure
           * we do this only if 'int' has at least 4 bytes.
           */
          return (uint32_t)(v >> 31);

#   else

          /* This should be a rather generic approach for getting a sign
           * extension mask...
           */
          return UINT32_C(0) - (uint32_t)(v < 0);

#   endif
}

static inline int32_t
uint32_2cpl_to_int32(
          const uint32_t vu)
{
          int32_t v;

#   if TARGET_HAS_2CPL

          /* Just copy through the 32 bits from the unsigned value if
           * we're on a two's complement target.
           */
          v = (int32_t)vu;

#   else

          /* Convert to signed integer, making sure signed integer
           * overflow cannot happen. Again, the optimiser might or might
           * not find out that this is just a copy of 32 bits on a target
           * with two's complement representation for signed integers.
           */
          if (vu > INT32_MAX)
                    v = -(int32_t)(~vu) - 1;
          else
                    v = (int32_t)vu;

#   endif

          return v;
}

/*
 *---------------------------------------------------------------------
 * Convert between 'time_t' and 'vint64'
 *---------------------------------------------------------------------
 */
vint64
time_to_vint64(
          const time_t * ptt
          )
{
          vint64 res;
          time_t tt;

          tt = *ptt;

#   if SIZEOF_TIME_T <= 4

          res.D_s.hi = 0;
          if (tt < 0) {
                    res.D_s.lo = (uint32_t)-tt;
                    M_NEG(res.D_s.hi, res.D_s.lo);
          } else {
                    res.D_s.lo = (uint32_t)tt;
          }

#   elif defined(HAVE_INT64)

          res.q_s = tt;

#   else
          /*
           * shifting negative signed quantities is compiler-dependent, so
           * we better avoid it and do it all manually. And shifting more
           * than the width of a quantity is undefined. Also a don't do!
           */
          if (tt < 0) {
                    tt = -tt;
                    res.D_s.lo = (uint32_t)tt;
                    res.D_s.hi = (uint32_t)(tt >> 32);
                    M_NEG(res.D_s.hi, res.D_s.lo);
          } else {
                    res.D_s.lo = (uint32_t)tt;
                    res.D_s.hi = (uint32_t)(tt >> 32);
          }

#   endif

          return res;
}


time_t
vint64_to_time(
          const vint64 *tv
          )
{
          time_t res;

#   if SIZEOF_TIME_T <= 4

          res = (time_t)tv->D_s.lo;

#   elif defined(HAVE_INT64)

          res = (time_t)tv->q_s;

#   else

          res = ((time_t)tv->d_s.hi << 32) | tv->D_s.lo;

#   endif

          return res;
}

/*
 *---------------------------------------------------------------------
 * Get the build date & time
 *---------------------------------------------------------------------
 */
int
ntpcal_get_build_date(
          struct calendar * jd
          )
{
          /* The C standard tells us the format of '__DATE__':
           *
           * __DATE__ The date of translation of the preprocessing
           * translation unit: a character string literal of the form "Mmm
           * dd yyyy", where the names of the months are the same as those
           * generated by the asctime function, and the first character of
           * dd is a space character if the value is less than 10. If the
           * date of translation is not available, an
           * implementation-defined valid date shall be supplied.
           *
           * __TIME__ The time of translation of the preprocessing
           * translation unit: a character string literal of the form
           * "hh:mm:ss" as in the time generated by the asctime
           * function. If the time of translation is not available, an
           * implementation-defined valid time shall be supplied.
           *
           * Note that MSVC declares DATE and TIME to be in the local time
           * zone, while neither the C standard nor the GCC docs make any
           * statement about this. As a result, we may be +/-12hrs off
           * UTC.    But for practical purposes, this should not be a
           * problem.
           *
           */
#   ifdef MKREPRO_DATE
          static const char build[] = MKREPRO_TIME "/" MKREPRO_DATE;
#   else
          static const char build[] = __TIME__ "/" __DATE__;
#   endif
          static const char mlist[] = "JanFebMarAprMayJunJulAugSepOctNovDec";

          char                  monstr[4];
          const char *          cp;
          unsigned short        hour, minute, second, day, year;
          /* Note: The above quantities are used for sscanf 'hu' format,
           * so using 'uint16_t' is contra-indicated!
           */

#   ifdef DEBUG
          static int            ignore  = 0;
#   endif

          ZERO(*jd);
          jd->year     = 1970;
          jd->month    = 1;
          jd->monthday = 1;

#   ifdef DEBUG
          /* check environment if build date should be ignored */
          if (0 == ignore) {
              const char * envstr;
              envstr = getenv("NTPD_IGNORE_BUILD_DATE");
              ignore = 1 + (envstr && (!*envstr || !strcasecmp(envstr, "yes")));
          }
          if (ignore > 1)
              return FALSE;
#   endif

          if (6 == sscanf(build, "%hu:%hu:%hu/%3s %hu %hu",
                              &hour, &minute, &second, monstr, &day, &year)) {
                    cp = strstr(mlist, monstr);
                    if (NULL != cp) {
                              jd->year     = year;
                              jd->month    = (uint8_t)((cp - mlist) / 3 + 1);
                              jd->monthday = (uint8_t)day;
                              jd->hour     = (uint8_t)hour;
                              jd->minute   = (uint8_t)minute;
                              jd->second   = (uint8_t)second;

                              return TRUE;
                    }
          }

          return FALSE;
}


/*
 *---------------------------------------------------------------------
 * basic calendar stuff
 *---------------------------------------------------------------------
 */

/*
 * Some notes on the terminology:
 *
 * We use the proleptic Gregorian calendar, which is the Gregorian
 * calendar extended in both directions ad infinitum. This totally
 * disregards the fact that this calendar was invented in 1582, and
 * was adopted at various dates over the world; sometimes even after
 * the start of the NTP epoch.
 *
 * Normally date parts are given as current cycles, while time parts
 * are given as elapsed cycles:
 *
 * 1970-01-01/03:04:05 means 'IN the 1970st. year, IN the first month,
 * ON the first day, with 3hrs, 4minutes and 5 seconds elapsed.
 *
 * The basic calculations for this calendar implementation deal with
 * ELAPSED date units, which is the number of full years, full months
 * and full days before a date: 1970-01-01 would be (1969, 0, 0) in
 * that notation.
 *
 * To ease the numeric computations, month and day values outside the
 * normal range are acceptable: 2001-03-00 will be treated as the day
 * before 2001-03-01, 2000-13-32 will give the same result as
 * 2001-02-01 and so on.
 *
 * 'rd' or 'RD' is used as an abbreviation for the latin 'rata die'
 * (day number).  This is the number of days elapsed since 0000-12-31
 * in the proleptic Gregorian calendar. The begin of the Christian Era
 * (0001-01-01) is RD(1).
 */

/*
 * ====================================================================
 *
 * General algorithmic stuff
 *
 * ====================================================================
 */

/*
 *---------------------------------------------------------------------
 * fast modulo 7 operations (floor/mathematical convention)
 *---------------------------------------------------------------------
 */
int
u32mod7(
          uint32_t x
          )
{
          /* This is a combination of tricks from "Hacker's Delight" with
           * some modifications, like a multiplication that rounds up to
           * drop the final adjustment stage.
           *
           * Do a partial reduction by digit sum to keep the value in the
           * range permitted for the mul/shift stage. There are several
           * possible and absolutely equivalent shift/mask combinations;
           * this one is ARM-friendly because of a mask that fits into 16
           * bit.
           */
          x = (x >> 15) + (x & UINT32_C(0x7FFF));
          /* Take reminder as (mod 8) by mul/shift. Since the multiplier
           * was calculated using ceil() instead of floor(), it skips the
           * value '7' properly.
           *    M <- ceil(ldexp(8/7, 29))
           */
          return (int)((x * UINT32_C(0x24924925)) >> 29);
}

int
i32mod7(
          int32_t x
          )
{
          /* We add (2**32 - 2**32 % 7), which is (2**32 - 4), to negative
           * numbers to map them into the postive range. Only the term '-4'
           * survives, obviously.
           */
          uint32_t ux = (uint32_t)x;
          return u32mod7((x < 0) ? (ux - 4u) : ux);
}

uint32_t
i32fmod(
          int32_t    x,
          uint32_t d
          )
{
          uint32_t ux = (uint32_t)x;
          uint32_t sf = UINT32_C(0) - (x < 0);
          ux = (sf ^ ux ) % d;
          return (d & sf) + (sf ^ ux);
}

/*
 *---------------------------------------------------------------------
 * Do a periodic extension of 'value' around 'pivot' with a period of
 * 'cycle'.
 *
 * The result 'res' is a number that holds to the following properties:
 *
 *   1)    res MOD cycle == value MOD cycle
 *   2)    pivot <= res < pivot + cycle
 *         (replace </<= with >/>= for negative cycles)
 *
 * where 'MOD' denotes the modulo operator for FLOOR DIVISION, which
 * is not the same as the '%' operator in C: C requires division to be
 * a truncated division, where remainder and dividend have the same
 * sign if the remainder is not zero, whereas floor division requires
 * divider and modulus to have the same sign for a non-zero modulus.
 *
 * This function has some useful applications:
 *
 * + let Y be a calendar year and V a truncated 2-digit year: then
 *        periodic_extend(Y-50, V, 100)
 *   is the closest expansion of the truncated year with respect to
 *   the full year, that is a 4-digit year with a difference of less
 *   than 50 years to the year Y. ("century unfolding")
 *
 * + let T be a UN*X time stamp and V be seconds-of-day: then
 *        perodic_extend(T-43200, V, 86400)
 *   is a time stamp that has the same seconds-of-day as the input
 *   value, with an absolute difference to T of <= 12hrs.  ("day
 *   unfolding")
 *
 * + Wherever you have a truncated periodic value and a non-truncated
 *   base value and you want to match them somehow...
 *
 * Basically, the function delivers 'pivot + (value - pivot) % cycle',
 * but the implementation takes some pains to avoid internal signed
 * integer overflows in the '(value - pivot) % cycle' part and adheres
 * to the floor division convention.
 *
 * If 64bit scalars where available on all intended platforms, writing a
 * version that uses 64 bit ops would be easy; writing a general
 * division routine for 64bit ops on a platform that can only do
 * 32/16bit divisions and is still performant is a bit more
 * difficult. Since most usecases can be coded in a way that does only
 * require the 32bit version a 64bit version is NOT provided here.
 *---------------------------------------------------------------------
 */
int32_t
ntpcal_periodic_extend(
          int32_t pivot,
          int32_t value,
          int32_t cycle
          )
{
          /* Implement a 4-quadrant modulus calculation by 2 2-quadrant
           * branches, one for positive and one for negative dividers.
           * Everything else can be handled by bit level logic and
           * conditional one's complement arithmetic.  By convention, we
           * assume
           *
           * x % b == 0  if  |b| < 2
           *
           * that is, we don't actually divide for cycles of -1,0,1 and
           * return the pivot value in that case.
           */
          uint32_t  uv = (uint32_t)value;
          uint32_t  up = (uint32_t)pivot;
          uint32_t  uc, sf;

          if (cycle > 1)
          {
                    uc = (uint32_t)cycle;
                    sf = UINT32_C(0) - (value < pivot);

                    uv = sf ^ (uv - up);
                    uv %= uc;
                    pivot += (uc & sf) + (sf ^ uv);
          }
          else if (cycle < -1)
          {
                    uc = ~(uint32_t)cycle + 1;
                    sf = UINT32_C(0) - (value > pivot);

                    uv = sf ^ (up - uv);
                    uv %= uc;
                    pivot -= (uc & sf) + (sf ^ uv);
          }
          return pivot;
}

/*---------------------------------------------------------------------
 * Note to the casual reader
 *
 * In the next two functions you will find (or would have found...)
 * the expression
 *
 *   res.Q_s -= 0x80000000;
 *
 * There was some ruckus about a possible programming error due to
 * integer overflow and sign propagation.
 *
 * This assumption is based on a lack of understanding of the C
 * standard. (Though this is admittedly not one of the most 'natural'
 * aspects of the 'C' language and easily to get wrong.)
 *
 * see
 *        http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1570.pdf
 *        "ISO/IEC 9899:201x Committee Draft — April 12, 2011"
 *        6.4.4.1 Integer constants, clause 5
 *
 * why there is no sign extension/overflow problem here.
 *
 * But to ease the minds of the doubtful, I added back the 'u' qualifiers
 * that somehow got lost over the last years.
 */


/*
 *---------------------------------------------------------------------
 * Convert a timestamp in NTP scale to a 64bit seconds value in the UN*X
 * scale with proper epoch unfolding around a given pivot or the current
 * system time. This function happily accepts negative pivot values as
 * timestamps before 1970-01-01, so be aware of possible trouble on
 * platforms with 32bit 'time_t'!
 *
 * This is also a periodic extension, but since the cycle is 2^32 and
 * the shift is 2^31, we can do some *very* fast math without explicit
 * divisions.
 *---------------------------------------------------------------------
 */
vint64
ntpcal_ntp_to_time(
          uint32_t  ntp,
          const time_t *      pivot
          )
{
          vint64 res;

#   if defined(HAVE_INT64)

          res.q_s = (pivot != NULL)
                          ? *pivot
                          : now();
          res.Q_s -= 0x80000000u;                 /* unshift of half range */
          ntp       -= (uint32_t)JAN_1970;        /* warp into UN*X domain */
          ntp       -= res.D_s.lo;                /* cycle difference  */
          res.Q_s += (uint64_t)ntp;     /* get expanded time           */

#   else /* no 64bit scalars */

          time_t tmp;

          tmp = (pivot != NULL)
                      ? *pivot
                      : now();
          res = time_to_vint64(&tmp);
          M_SUB(res.D_s.hi, res.D_s.lo, 0, 0x80000000u);
          ntp -= (uint32_t)JAN_1970;    /* warp into UN*X domain */
          ntp -= res.D_s.lo;            /* cycle difference  */
          M_ADD(res.D_s.hi, res.D_s.lo, 0, ntp);

#   endif /* no 64bit scalars */

          return res;
}

/*
 *---------------------------------------------------------------------
 * Convert a timestamp in NTP scale to a 64bit seconds value in the NTP
 * scale with proper epoch unfolding around a given pivot or the current
 * system time.
 *
 * Note: The pivot must be given in the UN*X time domain!
 *
 * This is also a periodic extension, but since the cycle is 2^32 and
 * the shift is 2^31, we can do some *very* fast math without explicit
 * divisions.
 *---------------------------------------------------------------------
 */
vint64
ntpcal_ntp_to_ntp(
          uint32_t      ntp,
          const time_t *pivot
          )
{
          vint64 res;

#   if defined(HAVE_INT64)

          res.q_s = (pivot)
                          ? *pivot
                          : now();
          res.Q_s -= 0x80000000u;                 /* unshift of half range */
          res.Q_s += (uint32_t)JAN_1970;          /* warp into NTP domain        */
          ntp       -= res.D_s.lo;                /* cycle difference  */
          res.Q_s += (uint64_t)ntp;     /* get expanded time           */

#   else /* no 64bit scalars */

          time_t tmp;

          tmp = (pivot)
                      ? *pivot
                      : now();
          res = time_to_vint64(&tmp);
          M_SUB(res.D_s.hi, res.D_s.lo, 0, 0x80000000u);
          M_ADD(res.D_s.hi, res.D_s.lo, 0, (uint32_t)JAN_1970);/*into NTP */
          ntp -= res.D_s.lo;            /* cycle difference  */
          M_ADD(res.D_s.hi, res.D_s.lo, 0, ntp);

#   endif /* no 64bit scalars */

          return res;
}


/*
 * ====================================================================
 *
 * Splitting values to composite entities
 *
 * ====================================================================
 */

/*
 *---------------------------------------------------------------------
 * Split a 64bit seconds value into elapsed days in 'res.hi' and
 * elapsed seconds since midnight in 'res.lo' using explicit floor
 * division. This function happily accepts negative time values as
 * timestamps before the respective epoch start.
 *---------------------------------------------------------------------
 */
ntpcal_split
ntpcal_daysplit(
          const vint64 *ts
          )
{
          ntpcal_split res;
          uint32_t Q, R;

#   if defined(HAVE_64BITREGS)

          /* Assume we have 64bit registers an can do a divison by
           * constant reasonably fast using the one's complement trick..
           */
          uint64_t sf64 = (uint64_t)-(ts->q_s < 0);
          Q = (uint32_t)(sf64 ^ ((sf64 ^ ts->Q_s) / SECSPERDAY));
          R = (uint32_t)(ts->Q_s - Q * SECSPERDAY);

#   elif defined(UINT64_MAX) && !defined(__arm__)

          /* We rely on the compiler to do efficient 64bit divisions as
           * good as possible. Which might or might not be true. At least
           * for ARM CPUs, the sum-by-digit code in the next section is
           * faster for many compilers. (This might change over time, but
           * the 64bit-by-32bit division will never outperform the exact
           * division by a substantial factor....)
           */
          if (ts->q_s < 0)
                    Q = ~(uint32_t)(~ts->Q_s / SECSPERDAY);
          else
                    Q =  (uint32_t)( ts->Q_s / SECSPERDAY);
          R = ts->D_s.lo - Q * SECSPERDAY;

#   else

          /* We don't have 64bit regs. That hurts a bit.
           *
           * Here we use a mean trick to get away with just one explicit
           * modulo operation and pure 32bit ops.
           *
           * Remember: 86400 <--> 128 * 675
           *
           * So we discard the lowest 7 bit and do an exact division by
           * 675, modulo 2**32.
           *
           * First we shift out the lower 7 bits.
           *
           * Then we use a digit-wise pseudo-reduction, where a 'digit' is
           * actually a 16-bit group. This is followed by a full reduction
           * with a 'true' division step. This yields the modulus of the
           * full 64bit value. The sign bit gets some extra treatment.
           *
           * Then we decrement the lower limb by that modulus, so it is
           * exactly divisible by 675. [*]
           *
           * Then we multiply with the modular inverse of 675 (mod 2**32)
           * and voila, we have the result.
           *
           * Special Thanks to Henry S. Warren and his "Hacker's delight"
           * for giving that idea.
           *
           * (Note[*]: that's not the full truth. We would have to
           * subtract the modulus from the full 64 bit number to get a
           * number that is divisible by 675. But since we use the
           * multiplicative inverse (mod 2**32) there's no reason to carry
           * the subtraction into the upper bits!)
           */
          uint32_t al = ts->D_s.lo;
          uint32_t ah = ts->D_s.hi;

          /* shift out the lower 7 bits, smash sign bit */
          al = (al >> 7) | (ah << 25);
          ah = (ah >> 7) & 0x00FFFFFFu;

          R  = (ts->d_s.hi < 0) ? 239 : 0;/* sign bit value */
          R += (al & 0xFFFF);
          R += (al >> 16       ) * 61u; /* 2**16 % 675 */
          R += (ah & 0xFFFF) * 346u;    /* 2**32 % 675 */
          R += (ah >> 16       ) * 181u;          /* 2**48 % 675 */
          R %= 675u;                              /* final reduction */
          Q  = (al - R) * 0x2D21C10Bu;  /* modinv(675, 2**32) */
          R  = (R << 7) | (ts->d_s.lo & 0x07F);

#   endif

          res.hi = uint32_2cpl_to_int32(Q);
          res.lo = R;

          return res;
}

/*
 *---------------------------------------------------------------------
 * Split a 64bit seconds value into elapsed weeks in 'res.hi' and
 * elapsed seconds since week start in 'res.lo' using explicit floor
 * division. This function happily accepts negative time values as
 * timestamps before the respective epoch start.
 *---------------------------------------------------------------------
 */
ntpcal_split
ntpcal_weeksplit(
          const vint64 *ts
          )
{
          ntpcal_split res;
          uint32_t Q, R;

          /* This is a very close relative to the day split function; for
           * details, see there!
           */

#   if defined(HAVE_64BITREGS)

          uint64_t sf64 = (uint64_t)-(ts->q_s < 0);
          Q = (uint32_t)(sf64 ^ ((sf64 ^ ts->Q_s) / SECSPERWEEK));
          R = (uint32_t)(ts->Q_s - Q * SECSPERWEEK);

#   elif defined(UINT64_MAX) && !defined(__arm__)

          if (ts->q_s < 0)
                    Q = ~(uint32_t)(~ts->Q_s / SECSPERWEEK);
          else
                    Q =  (uint32_t)( ts->Q_s / SECSPERWEEK);
          R = ts->D_s.lo - Q * SECSPERWEEK;

#   else

          /* Remember: 7*86400 <--> 604800 <--> 128 * 4725 */
          uint32_t al = ts->D_s.lo;
          uint32_t ah = ts->D_s.hi;

          al = (al >> 7) | (ah << 25);
          ah = (ah >> 7) & 0x00FFFFFF;

          R  = (ts->d_s.hi < 0) ? 2264 : 0;/* sign bit value */
          R += (al & 0xFFFF);
          R += (al >> 16       ) * 4111u;         /* 2**16 % 4725 */
          R += (ah & 0xFFFF) * 3721u;   /* 2**32 % 4725 */
          R += (ah >> 16       ) * 2206u;         /* 2**48 % 4725 */
          R %= 4725u;                             /* final reduction */
          Q  = (al - R) * 0x98BBADDDu;  /* modinv(4725, 2**32) */
          R  = (R << 7) | (ts->d_s.lo & 0x07F);

#   endif

          res.hi = uint32_2cpl_to_int32(Q);
          res.lo = R;

          return res;
}

/*
 *---------------------------------------------------------------------
 * Split a 32bit seconds value into h/m/s and excessive days.  This
 * function happily accepts negative time values as timestamps before
 * midnight.
 *---------------------------------------------------------------------
 */
static int32_t
priv_timesplit(
          int32_t split[3],
          int32_t ts
          )
{
          /* Do 3 chained floor divisions by positive constants, using the
           * one's complement trick and factoring out the intermediate XOR
           * ops to reduce the number of operations.
           */
          uint32_t us, um, uh, ud, sf32;

          sf32 = int32_sflag(ts);

          us = (uint32_t)ts;
          um = (sf32 ^ us) / SECSPERMIN;
          uh = um / MINSPERHR;
          ud = uh / HRSPERDAY;

          um ^= sf32;
          uh ^= sf32;
          ud ^= sf32;

          split[0] = (int32_t)(uh - ud * HRSPERDAY );
          split[1] = (int32_t)(um - uh * MINSPERHR );
          split[2] = (int32_t)(us - um * SECSPERMIN);

          return uint32_2cpl_to_int32(ud);
}

/*
 *---------------------------------------------------------------------
 * Given the number of elapsed days in the calendar era, split this
 * number into the number of elapsed years in 'res.hi' and the number
 * of elapsed days of that year in 'res.lo'.
 *
 * if 'isleapyear' is not NULL, it will receive an integer that is 0 for
 * regular years and a non-zero value for leap years.
 *---------------------------------------------------------------------
 */
ntpcal_split
ntpcal_split_eradays(
          int32_t days,
          int  *isleapyear
          )
{
          /* Use the fast cycle split algorithm here, to calculate the
           * centuries and years in a century with one division each. This
           * reduces the number of division operations to two, but is
           * susceptible to internal range overflow. We take some extra
           * steps to avoid the gap.
           */
          ntpcal_split res;
          int32_t    n100, n001; /* calendar year cycles */
          uint32_t uday, Q;

          /* split off centuries first
           *
           * We want to execute '(days * 4 + 3) /% 146097' under floor
           * division rules in the first step. Well, actually we want to
           * calculate 'floor((days + 0.75) / 36524.25)', but we want to
           * do it in scaled integer calculation.
           */
#   if defined(HAVE_64BITREGS)

          /* not too complicated with an intermediate 64bit value */
          uint64_t  ud64, sf64;
          ud64 = ((uint64_t)days << 2) | 3u;
          sf64 = (uint64_t)-(days < 0);
          Q    = (uint32_t)(sf64 ^ ((sf64 ^ ud64) / GREGORIAN_CYCLE_DAYS));
          uday = (uint32_t)(ud64 - Q * GREGORIAN_CYCLE_DAYS);
          n100 = uint32_2cpl_to_int32(Q);

#   else

          /* '4*days+3' suffers from range overflow when going to the
           * limits. We solve this by doing an exact division (mod 2^32)
           * after caclulating the remainder first.
           *
           * We start with a partial reduction by digit sums, extracting
           * the upper bits from the original value before they get lost
           * by scaling, and do one full division step to get the true
           * remainder.  Then a final multiplication with the
           * multiplicative inverse of 146097 (mod 2^32) gives us the full
           * quotient.
           *
           * (-2^33) % 146097 --> 130717    : the sign bit value
           * ( 2^20) % 146097 --> 25897     : the upper digit value
           * modinv(146097, 2^32) --> 660721233 : the inverse
           */
          uint32_t ux = ((uint32_t)days << 2) | 3;
          uday  = (days < 0) ? 130717u : 0u;          /* sign dgt */
          uday += ((days >> 18) & 0x01FFFu) * 25897u; /* hi dgt (src!) */
          uday += (ux & 0xFFFFFu);                    /* lo dgt */
          uday %= GREGORIAN_CYCLE_DAYS;               /* full reduction */
          Q     = (ux  - uday) * 660721233u;          /* exact div */
          n100  = uint32_2cpl_to_int32(Q);

#   endif

          /* Split off years in century -- days >= 0 here, and we're far
           * away from integer overflow trouble now. */
          uday |= 3;
          n001  = uday / GREGORIAN_NORMAL_LEAP_CYCLE_DAYS;
          uday -= n001 * GREGORIAN_NORMAL_LEAP_CYCLE_DAYS;

          /* Assemble the year and day in year */
          res.hi = n100 * 100 + n001;
          res.lo = uday / 4u;

          /* Possibly set the leap year flag */
          if (isleapyear) {
                    uint32_t tc = (uint32_t)n100 + 1;
                    uint32_t ty = (uint32_t)n001 + 1;
                    *isleapyear = !(ty & 3)
                        && ((ty != 100) || !(tc & 3));
          }
          return res;
}

/*
 *---------------------------------------------------------------------
 * Given a number of elapsed days in a year and a leap year indicator,
 * split the number of elapsed days into the number of elapsed months in
 * 'res.hi' and the number of elapsed days of that month in 'res.lo'.
 *
 * This function will fail and return {-1,-1} if the number of elapsed
 * days is not in the valid range!
 *---------------------------------------------------------------------
 */
ntpcal_split
ntpcal_split_yeardays(
          int32_t eyd,
          int       isleap
          )
{
          /* Use the unshifted-year, February-with-30-days approach here.
           * Fractional interpolations are used in both directions, with
           * the smallest power-of-two divider to avoid any true division.
           */
          ntpcal_split        res = {-1, -1};

          /* convert 'isleap' to number of defective days */
          isleap = 1 + !isleap;
          /* adjust for February of 30 nominal days */
          if (eyd >= 61 - isleap)
                    eyd += isleap;
          /* if in range, convert to months and days in month */
          if (eyd >= 0 && eyd < 367) {
                    res.hi = (eyd * 67 + 32) >> 11;
                    res.lo = eyd - ((489 * res.hi + 8) >> 4);
          }

          return res;
}

/*
 *---------------------------------------------------------------------
 * Convert a RD into the date part of a 'struct calendar'.
 *---------------------------------------------------------------------
 */
int
ntpcal_rd_to_date(
          struct calendar *jd,
          int32_t              rd
          )
{
          ntpcal_split split;
          int            leapy;
          u_int          ymask;

          /* Get day-of-week first. It's simply the RD (mod 7)... */
          jd->weekday = i32mod7(rd);

          split = ntpcal_split_eradays(rd - 1, &leapy);
          /* Get year and day-of-year, with overflow check. If any of the
           * upper 16 bits is set after shifting to unity-based years, we
           * will have an overflow when converting to an unsigned 16bit
           * year. Shifting to the right is OK here, since it does not
           * matter if the shift is logic or arithmetic.
           */
          split.hi += 1;
          ymask = 0u - ((split.hi >> 16) == 0);
          jd->year = (uint16_t)(split.hi & ymask);
          jd->yearday = (uint16_t)split.lo + 1;

          /* convert to month and mday */
          split = ntpcal_split_yeardays(split.lo, leapy);
          jd->month    = (uint8_t)split.hi + 1;
          jd->monthday = (uint8_t)split.lo + 1;

          return ymask ? leapy : -1;
}

/*
 *---------------------------------------------------------------------
 * Convert a RD into the date part of a 'struct tm'.
 *---------------------------------------------------------------------
 */
int
ntpcal_rd_to_tm(
          struct tm  *utm,
          int32_t       rd
          )
{
          ntpcal_split split;
          int            leapy;

          /* get day-of-week first */
          utm->tm_wday = i32mod7(rd);

          /* get year and day-of-year */
          split = ntpcal_split_eradays(rd - 1, &leapy);
          utm->tm_year = split.hi - 1899;
          utm->tm_yday = split.lo;      /* 0-based */

          /* convert to month and mday */
          split = ntpcal_split_yeardays(split.lo, leapy);
          utm->tm_mon  = split.hi;      /* 0-based */
          utm->tm_mday = split.lo + 1;  /* 1-based */

          return leapy;
}

/*
 *---------------------------------------------------------------------
 * Take a value of seconds since midnight and split it into hhmmss in a
 * 'struct calendar'.
 *---------------------------------------------------------------------
 */
int32_t
ntpcal_daysec_to_date(
          struct calendar *jd,
          int32_t             sec
          )
{
          int32_t days;
          int   ts[3];

          days = priv_timesplit(ts, sec);
          jd->hour   = (uint8_t)ts[0];
          jd->minute = (uint8_t)ts[1];
          jd->second = (uint8_t)ts[2];

          return days;
}

/*
 *---------------------------------------------------------------------
 * Take a value of seconds since midnight and split it into hhmmss in a
 * 'struct tm'.
 *---------------------------------------------------------------------
 */
int32_t
ntpcal_daysec_to_tm(
          struct tm *utm,
          int32_t      sec
          )
{
          int32_t days;
          int32_t ts[3];

          days = priv_timesplit(ts, sec);
          utm->tm_hour = ts[0];
          utm->tm_min  = ts[1];
          utm->tm_sec  = ts[2];

          return days;
}

/*
 *---------------------------------------------------------------------
 * take a split representation for day/second-of-day and day offset
 * and convert it to a 'struct calendar'. The seconds will be normalised
 * into the range of a day, and the day will be adjusted accordingly.
 *
 * returns >0 if the result is in a leap year, 0 if in a regular
 * year and <0 if the result did not fit into the calendar struct.
 *---------------------------------------------------------------------
 */
int
ntpcal_daysplit_to_date(
          struct calendar        *jd,
          const ntpcal_split *ds,
          int32_t                 dof
          )
{
          dof += ntpcal_daysec_to_date(jd, ds->lo);
          return ntpcal_rd_to_date(jd, ds->hi + dof);
}

/*
 *---------------------------------------------------------------------
 * take a split representation for day/second-of-day and day offset
 * and convert it to a 'struct tm'. The seconds will be normalised
 * into the range of a day, and the day will be adjusted accordingly.
 *
 * returns 1 if the result is in a leap year and zero if in a regular
 * year.
 *---------------------------------------------------------------------
 */
int
ntpcal_daysplit_to_tm(
          struct tm    *utm,
          const ntpcal_split *ds ,
          int32_t                 dof
          )
{
          dof += ntpcal_daysec_to_tm(utm, ds->lo);

          return ntpcal_rd_to_tm(utm, ds->hi + dof);
}

/*
 *---------------------------------------------------------------------
 * Take a UN*X time and convert to a calendar structure.
 *---------------------------------------------------------------------
 */
int
ntpcal_time_to_date(
          struct calendar     *jd,
          const vint64        *ts
          )
{
          ntpcal_split ds;

          ds = ntpcal_daysplit(ts);
          ds.hi += ntpcal_daysec_to_date(jd, ds.lo);
          ds.hi += DAY_UNIX_STARTS;

          return ntpcal_rd_to_date(jd, ds.hi);
}


/*
 * ====================================================================
 *
 * merging composite entities
 *
 * ====================================================================
 */

#if !defined(HAVE_INT64)
/* multiplication helper. Seconds in days and weeks are multiples of 128,
 * and without that factor fit well into 16 bit. So a multiplication
 * of 32bit by 16bit and some shifting can be used on pure 32bit machines
 * with compilers that do not support 64bit integers.
 *
 * Calculate ( hi * mul * 128 ) + lo
 */
static vint64
_dwjoin(
          uint16_t  mul,
          int32_t             hi,
          int32_t             lo
          )
{
          vint64              res;
          uint32_t  p1, p2, sf;

          /* get sign flag and absolute value of 'hi' in p1 */
          sf = (uint32_t)-(hi < 0);
          p1 = ((uint32_t)hi + sf) ^ sf;

          /* assemble major units: res <- |hi| * mul */
          res.D_s.lo = (p1 & 0xFFFF) * mul;
          res.D_s.hi = 0;
          p1 = (p1 >> 16) * mul;
          p2 = p1 >> 16;
          p1 = p1 << 16;
          M_ADD(res.D_s.hi, res.D_s.lo, p2, p1);

          /* mul by 128, using shift: res <-- res << 7 */
          res.D_s.hi = (res.D_s.hi << 7) | (res.D_s.lo >> 25);
          res.D_s.lo = (res.D_s.lo << 7);

          /* fix up sign: res <-- (res + [sf|sf]) ^ [sf|sf] */
          M_ADD(res.D_s.hi, res.D_s.lo, sf, sf);
          res.D_s.lo ^= sf;
          res.D_s.hi ^= sf;

          /* properly add seconds: res <-- res + [sx(lo)|lo] */
          p2 = (uint32_t)-(lo < 0);
          p1 = (uint32_t)lo;
          M_ADD(res.D_s.hi, res.D_s.lo, p2, p1);
          return res;
}
#endif

/*
 *---------------------------------------------------------------------
 * Merge a number of days and a number of seconds into seconds,
 * expressed in 64 bits to avoid overflow.
 *---------------------------------------------------------------------
 */
vint64
ntpcal_dayjoin(
          int32_t days,
          int32_t secs
          )
{
          vint64 res;

#   if defined(HAVE_INT64)

          res.q_s    = days;
          res.q_s *= SECSPERDAY;
          res.q_s += secs;

#   else

          res = _dwjoin(675, days, secs);

#   endif

          return res;
}

/*
 *---------------------------------------------------------------------
 * Merge a number of weeks and a number of seconds into seconds,
 * expressed in 64 bits to avoid overflow.
 *---------------------------------------------------------------------
 */
vint64
ntpcal_weekjoin(
          int32_t week,
          int32_t secs
          )
{
          vint64 res;

#   if defined(HAVE_INT64)

          res.q_s    = week;
          res.q_s *= SECSPERWEEK;
          res.q_s += secs;

#   else

          res = _dwjoin(4725, week, secs);

#   endif

          return res;
}

/*
 *---------------------------------------------------------------------
 * get leap years since epoch in elapsed years
 *---------------------------------------------------------------------
 */
int32_t
ntpcal_leapyears_in_years(
          int32_t years
          )
{
          /* We use the in-out-in algorithm here, using the one's
           * complement division trick for negative numbers. The chained
           * division sequence by 4/25/4 gives the compiler the chance to
           * get away with only one true division and doing shifts otherwise.
           */

          uint32_t sf32, sum, uyear;

          sf32  = int32_sflag(years);
          uyear = (uint32_t)years;
          uyear ^= sf32;

          sum  = (uyear /=  4u);        /*   4yr rule --> IN  */
          sum -= (uyear /= 25u);        /* 100yr rule --> OUT */
          sum += (uyear /=  4u);        /* 400yr rule --> IN  */

          /* Thanks to the alternation of IN/OUT/IN we can do the sum
           * directly and have a single one's complement operation
           * here. (Only if the years are negative, of course.) Otherwise
           * the one's complement would have to be done when
           * adding/subtracting the terms.
           */
          return uint32_2cpl_to_int32(sf32 ^ sum);
}

/*
 *---------------------------------------------------------------------
 * Convert elapsed years in Era into elapsed days in Era.
 *---------------------------------------------------------------------
 */
int32_t
ntpcal_days_in_years(
          int32_t years
          )
{
          return years * DAYSPERYEAR + ntpcal_leapyears_in_years(years);
}

/*
 *---------------------------------------------------------------------
 * Convert a number of elapsed month in a year into elapsed days in year.
 *
 * The month will be normalized, and 'res.hi' will contain the
 * excessive years that must be considered when converting the years,
 * while 'res.lo' will contain the number of elapsed days since start
 * of the year.
 *
 * This code uses the shifted-month-approach to convert month to days,
 * because then there is no need to have explicit leap year
 * information.      The slight disadvantage is that for most month values
 * the result is a negative value, and the year excess is one; the
 * conversion is then simply based on the start of the following year.
 *---------------------------------------------------------------------
 */
ntpcal_split
ntpcal_days_in_months(
          int32_t m
          )
{
          ntpcal_split res;

          /* Add ten months with proper year adjustment. */
          if (m < 2) {
              res.lo  = m + 10;
              res.hi  = 0;
          } else {
              res.lo  = m - 2;
              res.hi  = 1;
          }

          /* Possibly normalise by floor division. This does not hapen for
           * input in normal range. */
          if (res.lo < 0 || res.lo >= 12) {
                    uint32_t mu, Q, sf32;
                    sf32 = int32_sflag(res.lo);
                    mu   = (uint32_t)res.lo;
                    Q    = sf32 ^ ((sf32 ^ mu) / 12u);

                    res.hi += uint32_2cpl_to_int32(Q);
                    res.lo    = mu - Q * 12u;
          }

          /* Get cummulated days in year with unshift. Use the fractional
           * interpolation with smallest possible power of two in the
           * divider.
           */
          res.lo = ((res.lo * 979 + 16) >> 5) - 306;

          return res;
}

/*
 *---------------------------------------------------------------------
 * Convert ELAPSED years/months/days of gregorian calendar to elapsed
 * days in Gregorian epoch.
 *
 * If you want to convert years and days-of-year, just give a month of
 * zero.
 *---------------------------------------------------------------------
 */
int32_t
ntpcal_edate_to_eradays(
          int32_t years,
          int32_t mons,
          int32_t mdays
          )
{
          ntpcal_split tmp;
          int32_t        res;

          if (mons) {
                    tmp = ntpcal_days_in_months(mons);
                    res = ntpcal_days_in_years(years + tmp.hi) + tmp.lo;
          } else
                    res = ntpcal_days_in_years(years);
          res += mdays;

          return res;
}

/*
 *---------------------------------------------------------------------
 * Convert ELAPSED years/months/days of gregorian calendar to elapsed
 * days in year.
 *
 * Note: This will give the true difference to the start of the given
 * year, even if months & days are off-scale.
 *---------------------------------------------------------------------
 */
int32_t
ntpcal_edate_to_yeardays(
          int32_t years,
          int32_t mons,
          int32_t mdays
          )
{
          ntpcal_split tmp;

          if (0 <= mons && mons < 12) {
                    if (mons >= 2)
                              mdays -= 2 - is_leapyear(years+1);
                    mdays += (489 * mons + 8) >> 4;
          } else {
                    tmp = ntpcal_days_in_months(mons);
                    mdays += tmp.lo
                           + ntpcal_days_in_years(years + tmp.hi)
                           - ntpcal_days_in_years(years);
          }

          return mdays;
}

/*
 *---------------------------------------------------------------------
 * Convert elapsed days and the hour/minute/second information into
 * total seconds.
 *
 * If 'isvalid' is not NULL, do a range check on the time specification
 * and tell if the time input is in the normal range, permitting for a
 * single leapsecond.
 *---------------------------------------------------------------------
 */
int32_t
ntpcal_etime_to_seconds(
          int32_t hours,
          int32_t minutes,
          int32_t seconds
          )
{
          int32_t res;

          res = (hours * MINSPERHR + minutes) * SECSPERMIN + seconds;

          return res;
}

/*
 *---------------------------------------------------------------------
 * Convert the date part of a 'struct tm' (that is, year, month,
 * day-of-month) into the RD of that day.
 *---------------------------------------------------------------------
 */
int32_t
ntpcal_tm_to_rd(
          const struct tm *utm
          )
{
          return ntpcal_edate_to_eradays(utm->tm_year + 1899,
                                               utm->tm_mon,
                                               utm->tm_mday - 1) + 1;
}

/*
 *---------------------------------------------------------------------
 * Convert the date part of a 'struct calendar' (that is, year, month,
 * day-of-month) into the RD of that day.
 *---------------------------------------------------------------------
 */
int32_t
ntpcal_date_to_rd(
          const struct calendar *jd
          )
{
          return ntpcal_edate_to_eradays((int32_t)jd->year - 1,
                                               (int32_t)jd->month - 1,
                                               (int32_t)jd->monthday - 1) + 1;
}

/*
 *---------------------------------------------------------------------
 * convert a year number to rata die of year start
 *---------------------------------------------------------------------
 */
int32_t
ntpcal_year_to_ystart(
          int32_t year
          )
{
          return ntpcal_days_in_years(year - 1) + 1;
}

/*
 *---------------------------------------------------------------------
 * For a given RD, get the RD of the associated year start,
 * that is, the RD of the last January,1st on or before that day.
 *---------------------------------------------------------------------
 */
int32_t
ntpcal_rd_to_ystart(
          int32_t rd
          )
{
          /*
           * Rather simple exercise: split the day number into elapsed
           * years and elapsed days, then remove the elapsed days from the
           * input value. Nice'n sweet...
           */
          return rd - ntpcal_split_eradays(rd - 1, NULL).lo;
}

/*
 *---------------------------------------------------------------------
 * For a given RD, get the RD of the associated month start.
 *---------------------------------------------------------------------
 */
int32_t
ntpcal_rd_to_mstart(
          int32_t rd
          )
{
          ntpcal_split split;
          int            leaps;

          split = ntpcal_split_eradays(rd - 1, &leaps);
          split = ntpcal_split_yeardays(split.lo, leaps);

          return rd - split.lo;
}

/*
 *---------------------------------------------------------------------
 * take a 'struct calendar' and get the seconds-of-day from it.
 *---------------------------------------------------------------------
 */
int32_t
ntpcal_date_to_daysec(
          const struct calendar *jd
          )
{
          return ntpcal_etime_to_seconds(jd->hour, jd->minute,
                                               jd->second);
}

/*
 *---------------------------------------------------------------------
 * take a 'struct tm' and get the seconds-of-day from it.
 *---------------------------------------------------------------------
 */
int32_t
ntpcal_tm_to_daysec(
          const struct tm *utm
          )
{
          return ntpcal_etime_to_seconds(utm->tm_hour, utm->tm_min,
                                               utm->tm_sec);
}

/*
 *---------------------------------------------------------------------
 * take a 'struct calendar' and convert it to a 'time_t'
 *---------------------------------------------------------------------
 */
time_t
ntpcal_date_to_time(
          const struct calendar *jd
          )
{
          vint64    join;
          int32_t days, secs;

          days = ntpcal_date_to_rd(jd) - DAY_UNIX_STARTS;
          secs = ntpcal_date_to_daysec(jd);
          join = ntpcal_dayjoin(days, secs);

          return vint64_to_time(&join);
}


/*
 * ====================================================================
 *
 * extended and unchecked variants of caljulian/caltontp
 *
 * ====================================================================
 */
int
ntpcal_ntp64_to_date(
          struct calendar *jd,
          const vint64        *ntp
          )
{
          ntpcal_split ds;

          ds = ntpcal_daysplit(ntp);
          ds.hi += ntpcal_daysec_to_date(jd, ds.lo);

          return ntpcal_rd_to_date(jd, ds.hi + DAY_NTP_STARTS);
}

int
ntpcal_ntp_to_date(
          struct calendar *jd,
          uint32_t   ntp,
          const time_t        *piv
          )
{
          vint64    ntp64;

          /*
           * Unfold ntp time around current time into NTP domain. Split
           * into days and seconds, shift days into CE domain and
           * process the parts.
           */
          ntp64 = ntpcal_ntp_to_ntp(ntp, piv);
          return ntpcal_ntp64_to_date(jd, &ntp64);
}


vint64
ntpcal_date_to_ntp64(
          const struct calendar *jd
          )
{
          /*
           * Convert date to NTP. Ignore yearday, use d/m/y only.
           */
          return ntpcal_dayjoin(ntpcal_date_to_rd(jd) - DAY_NTP_STARTS,
                                    ntpcal_date_to_daysec(jd));
}


uint32_t
ntpcal_date_to_ntp(
          const struct calendar *jd
          )
{
          /*
           * Get lower half of 64bit NTP timestamp from date/time.
           */
          return ntpcal_date_to_ntp64(jd).d_s.lo;
}



/*
 * ====================================================================
 *
 * day-of-week calculations
 *
 * ====================================================================
 */
/*
 * Given a RataDie and a day-of-week, calculate a RDN that is reater-than,
 * greater-or equal, closest, less-or-equal or less-than the given RDN
 * and denotes the given day-of-week
 */
int32_t
ntpcal_weekday_gt(
          int32_t rdn,
          int32_t dow
          )
{
          return ntpcal_periodic_extend(rdn+1, dow, 7);
}

int32_t
ntpcal_weekday_ge(
          int32_t rdn,
          int32_t dow
          )
{
          return ntpcal_periodic_extend(rdn, dow, 7);
}

int32_t
ntpcal_weekday_close(
          int32_t rdn,
          int32_t dow
          )
{
          return ntpcal_periodic_extend(rdn-3, dow, 7);
}

int32_t
ntpcal_weekday_le(
          int32_t rdn,
          int32_t dow
          )
{
          return ntpcal_periodic_extend(rdn, dow, -7);
}

int32_t
ntpcal_weekday_lt(
          int32_t rdn,
          int32_t dow
          )
{
          return ntpcal_periodic_extend(rdn-1, dow, -7);
}

/*
 * ====================================================================
 *
 * ISO week-calendar conversions
 *
 * The ISO8601 calendar defines a calendar of years, weeks and weekdays.
 * It is related to the Gregorian calendar, and a ISO year starts at the
 * Monday closest to Jan,1st of the corresponding Gregorian year.  A ISO
 * calendar year has always 52 or 53 weeks, and like the Grogrian
 * calendar the ISO8601 calendar repeats itself every 400 years, or
 * 146097 days, or 20871 weeks.
 *
 * While it is possible to write ISO calendar functions based on the
 * Gregorian calendar functions, the following implementation takes a
 * different approach, based directly on years and weeks.
 *
 * Analysis of the tabulated data shows that it is not possible to
 * interpolate from years to weeks over a full 400 year range; cyclic
 * shifts over 400 years do not provide a solution here. But it *is*
 * possible to interpolate over every single century of the 400-year
 * cycle. (The centennial leap year rule seems to be the culprit here.)
 *
 * It can be shown that a conversion from years to weeks can be done
 * using a linear transformation of the form
 *
 *   w = floor( y * a + b )
 *
 * where the slope a must hold to
 *
 *  52.1780821918 <= a < 52.1791044776
 *
 * and b must be chosen according to the selected slope and the number
 * of the century in a 400-year period.
 *
 * The inverse calculation can also be done in this way. Careful scaling
 * provides an unlimited set of integer coefficients a,k,b that enable
 * us to write the calulation in the form
 *
 *   w = (y * a      + b ) / k
 *   y = (w * a' + b') / k'
 *
 * In this implementation the values of k and k' are chosen to be the
 * smallest possible powers of two, so the division can be implemented
 * as shifts if the optimiser chooses to do so.
 *
 * ====================================================================
 */

/*
 * Given a number of elapsed (ISO-)years since the begin of the
 * christian era, return the number of elapsed weeks corresponding to
 * the number of years.
 */
int32_t
isocal_weeks_in_years(
          int32_t years
          )
{
          /*
           * use: w = (y * 53431 + b[c]) / 1024 as interpolation
           */
          static const uint16_t bctab[4] = { 157, 449, 597, 889 };

          int32_t    cs, cw;
          uint32_t cc, ci, yu, sf32;

          sf32 = int32_sflag(years);
          yu   = (uint32_t)years;

          /* split off centuries, using floor division */
          cc  = sf32 ^ ((sf32 ^ yu) / 100u);
          yu -= cc * 100u;

          /* calculate century cycles shift and cycle index:
           * Assuming a century is 5217 weeks, we have to add a cycle
           * shift that is 3 for every 4 centuries, because 3 of the four
           * centuries have 5218 weeks. So '(cc*3 + 1) / 4' is the actual
           * correction, and the second century is the defective one.
           *
           * Needs floor division by 4, which is done with masking and
           * shifting.
           */
          ci = cc * 3u + 1;
          cs = uint32_2cpl_to_int32(sf32 ^ ((sf32 ^ ci) >> 2));
          ci = ci & 3u;

          /* Get weeks in century. Can use plain division here as all ops
           * are >= 0,  and let the compiler sort out the possible
           * optimisations.
           */
          cw = (yu * 53431u + bctab[ci]) / 1024u;

          return uint32_2cpl_to_int32(cc) * 5217 + cs + cw;
}

/*
 * Given a number of elapsed weeks since the begin of the christian
 * era, split this number into the number of elapsed years in res.hi
 * and the excessive number of weeks in res.lo. (That is, res.lo is
 * the number of elapsed weeks in the remaining partial year.)
 */
ntpcal_split
isocal_split_eraweeks(
          int32_t weeks
          )
{
          /*
           * use: y = (w * 157 + b[c]) / 8192 as interpolation
           */

          static const uint16_t bctab[4] = { 85, 130, 17, 62 };

          ntpcal_split res;
          int32_t    cc, ci;
          uint32_t sw, cy, Q;

          /* Use two fast cycle-split divisions again. Herew e want to
           * execute '(weeks * 4 + 2) /% 20871' under floor division rules
           * in the first step.
           *
           * This is of course (again) susceptible to internal overflow if
           * coded directly in 32bit. And again we use 64bit division on
           * a 64bit target and exact division after calculating the
           * remainder first on a 32bit target. With the smaller divider,
           * that's even a bit neater.
           */
#   if defined(HAVE_64BITREGS)

          /* Full floor division with 64bit values. */
          uint64_t sf64, sw64;
          sf64 = (uint64_t)-(weeks < 0);
          sw64 = ((uint64_t)weeks << 2) | 2u;
          Q    = (uint32_t)(sf64 ^ ((sf64 ^ sw64) / GREGORIAN_CYCLE_WEEKS));
          sw   = (uint32_t)(sw64 - Q * GREGORIAN_CYCLE_WEEKS);

#   else

          /* Exact division after calculating the remainder via partial
           * reduction by digit sum.
           * (-2^33) % 20871     --> 5491              : the sign bit value
           * ( 2^20) % 20871     --> 5026              : the upper digit value
           * modinv(20871, 2^32) --> 330081335 : the inverse
           */
          uint32_t ux = ((uint32_t)weeks << 2) | 2;
          sw  = (weeks < 0) ? 5491u : 0u;                     /* sign dgt */
          sw += ((weeks >> 18) & 0x01FFFu) * 5026u; /* hi dgt (src!) */
          sw += (ux & 0xFFFFFu);                              /* lo dgt */
          sw %= GREGORIAN_CYCLE_WEEKS;              /* full reduction */
          Q   = (ux  - sw) * 330081335u;                      /* exact div */

#   endif

          ci  = Q & 3u;
          cc  = uint32_2cpl_to_int32(Q);

          /* Split off years; sw >= 0 here! The scaled weeks in the years
           * are scaled up by 157 afterwards.
           */
          sw  = (sw / 4u) * 157u + bctab[ci];
          cy  = sw / 8192u;   /* sw >> 13 , let the compiler sort it out */
          sw  = sw % 8192u;   /* sw & 8191, let the compiler sort it out */

          /* assemble elapsed years and downscale the elapsed weeks in
           * the year.
           */
          res.hi = 100*cc + cy;
          res.lo = sw / 157u;

          return res;
}

/*
 * Given a second in the NTP time scale and a pivot, expand the NTP
 * time stamp around the pivot and convert into an ISO calendar time
 * stamp.
 */
int
isocal_ntp64_to_date(
          struct isodate *id,
          const vint64   *ntp
          )
{
          ntpcal_split ds;
          int32_t        ts[3];
          uint32_t     uw, ud, sf32;

          /*
           * Split NTP time into days and seconds, shift days into CE
           * domain and process the parts.
           */
          ds = ntpcal_daysplit(ntp);

          /* split time part */
          ds.hi += priv_timesplit(ts, ds.lo);
          id->hour   = (uint8_t)ts[0];
          id->minute = (uint8_t)ts[1];
          id->second = (uint8_t)ts[2];

          /* split days into days and weeks, using floor division in unsigned */
          ds.hi += DAY_NTP_STARTS - 1; /* shift from NTP to RDN */
          sf32 = int32_sflag(ds.hi);
          ud   = (uint32_t)ds.hi;
          uw   = sf32 ^ ((sf32 ^ ud) / DAYSPERWEEK);
          ud  -= uw * DAYSPERWEEK;

          ds.hi = uint32_2cpl_to_int32(uw);
          ds.lo = ud;

          id->weekday = (uint8_t)ds.lo + 1;       /* weekday result    */

          /* get year and week in year */
          ds = isocal_split_eraweeks(ds.hi);      /* elapsed years&week*/
          id->year = (uint16_t)ds.hi + 1;                   /* shift to current  */
          id->week = (uint8_t )ds.lo + 1;

          return (ds.hi >= 0 && ds.hi < 0x0000FFFF);
}

int
isocal_ntp_to_date(
          struct isodate *id,
          uint32_t  ntp,
          const time_t   *piv
          )
{
          vint64    ntp64;

          /*
           * Unfold ntp time around current time into NTP domain, then
           * convert the full time stamp.
           */
          ntp64 = ntpcal_ntp_to_ntp(ntp, piv);
          return isocal_ntp64_to_date(id, &ntp64);
}

/*
 * Convert a ISO date spec into a second in the NTP time scale,
 * properly truncated to 32 bit.
 */
vint64
isocal_date_to_ntp64(
          const struct isodate *id
          )
{
          int32_t weeks, days, secs;

          weeks = isocal_weeks_in_years((int32_t)id->year - 1)
                + (int32_t)id->week - 1;
          days = weeks * 7 + (int32_t)id->weekday;
          /* days is RDN of ISO date now */
          secs = ntpcal_etime_to_seconds(id->hour, id->minute, id->second);

          return ntpcal_dayjoin(days - DAY_NTP_STARTS, secs);
}

uint32_t
isocal_date_to_ntp(
          const struct isodate *id
          )
{
          /*
           * Get lower half of 64bit NTP timestamp from date/time.
           */
          return isocal_date_to_ntp64(id).d_s.lo;
}

/*
 * ====================================================================
 * 'basedate' support functions
 * ====================================================================
 */

static int32_t s_baseday = NTP_TO_UNIX_DAYS;
static int32_t s_gpsweek = 0;

int32_t
basedate_eval_buildstamp(void)
{
          struct calendar jd;
          int32_t             ed;

          if (!ntpcal_get_build_date(&jd))
                    return NTP_TO_UNIX_DAYS;

          /* The time zone of the build stamp is unspecified; we remove
           * one day to provide a certain slack. And in case somebody
           * fiddled with the system clock, we make sure we do not go
           * before the UNIX epoch (1970-01-01). It's probably not possible
           * to do this to the clock on most systems, but there are other
           * ways to tweak the build stamp.
           */
          jd.monthday -= 1;
          ed = ntpcal_date_to_rd(&jd) - DAY_NTP_STARTS;
          return (ed < NTP_TO_UNIX_DAYS) ? NTP_TO_UNIX_DAYS : ed;
}

int32_t
basedate_eval_string(
          const char * str
          )
{
          u_short   y,m,d;
          u_long    ned;
          int       rc, nc;
          size_t    sl;

          sl = strlen(str);
          rc = sscanf(str, "%4hu-%2hu-%2hu%n", &y, &m, &d, &nc);
          if (rc == 3 && (size_t)nc == sl) {
                    if (m >= 1 && m <= 12 && d >= 1 && d <= 31)
                              return ntpcal_edate_to_eradays(y-1, m-1, d)
                                  - DAY_NTP_STARTS;
                    goto buildstamp;
          }

          rc = sscanf(str, "%lu%n", &ned, &nc);
          if (rc == 1 && (size_t)nc == sl) {
                    if (ned <= INT32_MAX)
                              return (int32_t)ned;
                    goto buildstamp;
          }

  buildstamp:
          msyslog(LOG_WARNING,
                    "basedate string \"%s\" invalid, build date substituted!",
                    str);
          return basedate_eval_buildstamp();
}

uint32_t
basedate_get_day(void)
{
          return s_baseday;
}

int32_t
basedate_set_day(
          int32_t day
          )
{
          struct calendar     jd;
          int32_t             retv;

          /* set NTP base date for NTP era unfolding */
          if (day < NTP_TO_UNIX_DAYS) {
                    msyslog(LOG_WARNING,
                              "baseday_set_day: invalid day (%lu), UNIX epoch substituted",
                              (unsigned long)day);
                    day = NTP_TO_UNIX_DAYS;
          }
          retv = s_baseday;
          s_baseday = day;
          ntpcal_rd_to_date(&jd, day + DAY_NTP_STARTS);
          msyslog(LOG_INFO, "basedate set to %04hu-%02hu-%02hu",
                    jd.year, (u_short)jd.month, (u_short)jd.monthday);

          /* set GPS base week for GPS week unfolding */
          day = ntpcal_weekday_ge(day + DAY_NTP_STARTS, CAL_SUNDAY)
              - DAY_NTP_STARTS;
          if (day < NTP_TO_GPS_DAYS)
              day = NTP_TO_GPS_DAYS;
          s_gpsweek = (day - NTP_TO_GPS_DAYS) / DAYSPERWEEK;
          ntpcal_rd_to_date(&jd, day + DAY_NTP_STARTS);
          msyslog(LOG_INFO, "gps base set to %04hu-%02hu-%02hu (week %d)",
                    jd.year, (u_short)jd.month, (u_short)jd.monthday, s_gpsweek);

          return retv;
}

time_t
basedate_get_eracenter(void)
{
          time_t retv;
          retv  = (time_t)(s_baseday - NTP_TO_UNIX_DAYS);
          retv *= SECSPERDAY;
          retv += (UINT32_C(1) << 31);
          return retv;
}

time_t
basedate_get_erabase(void)
{
          time_t retv;
          retv  = (time_t)(s_baseday - NTP_TO_UNIX_DAYS);
          retv *= SECSPERDAY;
          return retv;
}

uint32_t
basedate_get_gpsweek(void)
{
    return s_gpsweek;
}

uint32_t
basedate_expand_gpsweek(
    unsigned short weekno
    )
{
    /* We do a fast modulus expansion here. Since all quantities are
     * unsigned and we cannot go before the start of the GPS epoch
     * anyway, and since the truncated GPS week number is 10 bit, the
     * expansion becomes a simple sub/and/add sequence.
     */
    #if GPSWEEKS != 1024
    # error GPSWEEKS defined wrong -- should be 1024!
    #endif

    uint32_t diff;
    diff = ((uint32_t)weekno - s_gpsweek) & (GPSWEEKS - 1);
    return s_gpsweek + diff;
}

/*
 * ====================================================================
 * misc. helpers
 * ====================================================================
 */

/* --------------------------------------------------------------------
 * reconstruct the centrury from a truncated date and a day-of-week
 *
 * Given a date with truncated year (2-digit, 0..99) and a day-of-week
 * from 1(Mon) to 7(Sun), recover the full year between 1900AD and 2300AD.
 */
int32_t
ntpcal_expand_century(
          uint32_t y,
          uint32_t m,
          uint32_t d,
          uint32_t wd)
{
          /* This algorithm is short but tricky... It's related to
           * Zeller's congruence, partially done backwards.
           *
           * A few facts to remember:
           *  1) The Gregorian calendar has a cycle of 400 years.
           *  2) The weekday of the 1st day of a century shifts by 5 days
           *     during a great cycle.
           *  3) For calendar math, a century starts with the 1st year,
           *     which is year 1, !not! zero.
           *
           * So we start with taking the weekday difference (mod 7)
           * between the truncated date (which is taken as an absolute
           * date in the 1st century in the proleptic calendar) and the
           * weekday given.
           *
           * When dividing this residual by 5, we obtain the number of
           * centuries to add to the base. But since the residual is (mod
           * 7), we have to make this an exact division by multiplication
           * with the modular inverse of 5 (mod 7), which is 3:
           *    3*5 === 1 (mod 7).
           *
           * If this yields a result of 4/5/6, the given date/day-of-week
           * combination is impossible, and we return zero as resulting
           * year to indicate failure.
           *
           * Then we remap the century to the range starting with year
           * 1900.
           */

          uint32_t c;

          /* check basic constraints */
          if ((y >= 100u) || (--m >= 12u) || (--d >= 31u))
                    return 0;

          if ((m += 10u) >= 12u)                  /* shift base to prev. March,1st */
                    m -= 12u;
          else if (--y >= 100u)
                    y += 100u;
          d += y + (y >> 2) + 2u;                 /* year share */
          d += (m * 83u + 16u) >> 5;    /* month share */

          /* get (wd - d), shifted to positive value, and multiply with
           * 3(mod 7). (Exact division, see to comment)
           * Note: 1) d <= 184 at this point.
           *         2) 252 % 7 == 0, but 'wd' is off by one since we did
           *            '--d' above, so we add just 251 here!
           */
          c = u32mod7(3 * (251u + wd - d));
          if (c > 3u)
                    return 0;

          if ((m > 9u) && (++y >= 100u)) {/* undo base shift */
                    y -= 100u;
                    c = (c + 1) & 3u;
          }
          y += (c * 100u);              /* combine into 1st cycle */
          y += (y < 300u) ? 2000 : 1600;          /* map to destination era */
          return (int)y;
}

char *
ntpcal_iso8601std(
          char *              buf,
          size_t              len,
          TcCivilDate *       cdp
          )
{
          if (!buf) {
                    LIB_GETBUF(buf);
                    len = LIB_BUFLENGTH;
          }
          if (len) {
                    int slen = snprintf(buf, len, "%04u-%02u-%02uT%02u:%02u:%02u",
                                     cdp->year, cdp->month, cdp->monthday,
                                     cdp->hour, cdp->minute, cdp->second);
                    if (slen < 0)
                              *buf = '\0';
          }
          return buf;
}

/* -*-EOF-*- */