/*
* 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 suffer from an internal
* overflow, and these conditions are checked: errno is set to EDOM and
* the results are clamped/saturated in this case. All other functions
* do not suffer from internal overflow and simply return the result
* truncated to 32 bits.
*
* This is a sacrifice made for execution speed. Since a 32-bit day
* counter covers +/- 5,879,610 years and the clamp limits the effective
* range to +/-2.9 million years, this should not pose a problem here.
*
*/
#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"
/* 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
/*
*---------------------------------------------------------------------
* 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 uint32_t
int32_to_uint32_2cpl(
const int32_t v)
{
uint32_t vu;
# if TARGET_HAS_2CPL
/* Just copy through the 32 bits from the signed value if we're
* on a two's complement target.
*/
vu = (uint32_t)v;
# else
/* Convert from signed int to unsigned int two's complement. Do
* not make any assumptions about the representation of signed
* integers, but make sure signed integer overflow cannot happen
* here. A compiler on a two's complement target *might* find
* out that this is just a complicated cast (as above), but your
* mileage might vary.
*/
if (v < 0)
vu = ~(uint32_t)(-(v + 1));
else
vu = (uint32_t)v;
# endif
return vu;
}
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;
}
/* Some of the calculations need to multiply the input by 4 before doing
* a division. This can cause overflow and strange results. Therefore we
* clamp / saturate the input operand. And since we do the calculations
* in unsigned int with an extra sign flag/mask, we only loose one bit
* of the input value range.
*/
static inline uint32_t
uint32_saturate(
uint32_t vu,
uint32_t mu)
{
static const uint32_t limit = UINT32_MAX/4u;
if ((mu ^ vu) > limit) {
vu = mu ^ limit;
errno = EDOM;
}
return vu;
}
/*
*---------------------------------------------------------------------
* 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
*---------------------------------------------------------------------
*/
/* month table for a year starting with March,1st */
static const uint16_t shift_month_table[13] = {
0, 31, 61, 92, 122, 153, 184, 214, 245, 275, 306, 337, 366
};
/* month tables for years starting with January,1st; regular & leap */
static const uint16_t real_month_table[2][13] = {
/* -*- table for regular years -*- */
{ 0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334, 365 },
/* -*- table for leap years -*- */
{ 0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335, 366 }
};
/*
* 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
*
* ====================================================================
*/
/*
*---------------------------------------------------------------------
* 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 32-bit version a 64bit version is NOT provided here.
*---------------------------------------------------------------------
*/
int32_t
ntpcal_periodic_extend(
int32_t pivot,
int32_t value,
int32_t cycle
)
{
uint32_t diff;
char cpl = 0; /* modulo complement flag */
char neg = 0; /* sign change flag */
/* make the cycle positive and adjust the flags */
if (cycle < 0) {
cycle = - cycle;
neg ^= 1;
cpl ^= 1;
}
/* guard against div by zero or one */
if (cycle > 1) {
/*
* Get absolute difference as unsigned quantity and
* the complement flag. This is done by always
* subtracting the smaller value from the bigger
* one.
*/
if (value >= pivot) {
diff = int32_to_uint32_2cpl(value)
- int32_to_uint32_2cpl(pivot);
} else {
diff = int32_to_uint32_2cpl(pivot)
- int32_to_uint32_2cpl(value);
cpl ^= 1;
}
diff %= (uint32_t)cycle;
if (diff) {
if (cpl)
diff = (uint32_t)cycle - diff;
if (neg)
diff = ~diff + 1;
pivot += uint32_2cpl_to_int32(diff);
}
}
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 befor 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;
# if defined(HAVE_INT64)
/* Manual floor division by SECSPERDAY. This uses the one's
* complement trick, too, but without an extra flag value: The
* flag would be 64bit, and that's a bit of overkill on a 32bit
* target that has to use a register pair for a 64bit number.
*/
if (ts->q_s < 0)
Q = ~(uint32_t)(~ts->Q_s / SECSPERDAY);
else
Q = (uint32_t)(ts->Q_s / SECSPERDAY);
# else
uint32_t ah, al, sflag, A;
/* get operand into ah/al (either ts or ts' one's complement,
* for later floor division)
*/
sflag = int32_sflag(ts->d_s.hi);
ah = sflag ^ ts->D_s.hi;
al = sflag ^ ts->D_s.lo;
/* Since 86400 == 128*675 we can drop the least 7 bits and
* divide by 675 instead of 86400. Then the maximum remainder
* after each devision step is 674, and we need 10 bits for
* that. So in the next step we can shift in 22 bits from the
* numerator.
*
* Therefore we load the accu with the top 13 bits (51..63) in
* the first shot. We don't have to remember the quotient -- it
* would be shifted out anyway.
*/
A = ah >> 19;
if (A >= 675)
A = (A % 675u);
/* Now assemble the remainder with bits 29..50 from the
* numerator and divide. This creates the upper ten bits of the
* quotient. (Well, the top 22 bits of a 44bit result. But that
* will be truncated to 32 bits anyway.)
*/
A = (A << 19) | (ah & 0x0007FFFFu);
A = (A << 3) | (al >> 29);
Q = A / 675u;
A = A % 675u;
/* Now assemble the remainder with bits 7..28 from the numerator
* and do a final division step.
*/
A = (A << 22) | ((al >> 7) & 0x003FFFFFu);
Q = (Q << 22) | (A / 675u);
/* The last 7 bits get simply dropped, as they have no affect on
* the quotient when dividing by 86400.
*/
/* apply sign correction and calculate the true floor
* remainder.
*/
Q ^= sflag;
# endif
res.hi = uint32_2cpl_to_int32(Q);
res.lo = ts->D_s.lo - Q * SECSPERDAY;
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, sflag;
sflag = int32_sflag(ts);
us = int32_to_uint32_2cpl(ts);
um = (sflag ^ us) / SECSPERMIN;
uh = um / MINSPERHR;
ud = uh / HRSPERDAY;
um ^= sflag;
uh ^= sflag;
ud ^= sflag;
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 cyclesplit 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 make sure the
* input operands are in the safe range; this still gives us
* approx +/-2.9 million years.
*/
ntpcal_split res;
int32_t n100, n001; /* calendar year cycles */
uint32_t uday, Q, sflag;
/* split off centuries first */
sflag = int32_sflag(days);
uday = uint32_saturate(int32_to_uint32_2cpl(days), sflag);
uday = (4u * uday) | 3u;
Q = sflag ^ ((sflag ^ uday) / GREGORIAN_CYCLE_DAYS);
uday = uday - Q * GREGORIAN_CYCLE_DAYS;
n100 = uint32_2cpl_to_int32(Q);
/* 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 = uday % GREGORIAN_NORMAL_LEAP_CYCLE_DAYS;
/* Assemble the year and day in year */
res.hi = n100 * 100 + n001;
res.lo = uday / 4u;
/* Eventually set the leap year flag. Note: 0 <= n001 <= 99 and
* Q is still the two's complement representation of the
* centuries: The modulo 4 ops can be done with masking here.
* We also shift the year and the century by one, so the tests
* can be done against zero instead of 3.
*/
if (isleapyear)
*isleapyear = !((n001+1) & 3)
&& ((n001 != 99) || !((Q+1) & 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 isleapyear
)
{
ntpcal_split res;
const uint16_t *lt; /* month length table */
/* check leap year flag and select proper table */
lt = real_month_table[(isleapyear != 0)];
if (0 <= eyd && eyd < lt[12]) {
/* get zero-based month by approximation & correction step */
res.hi = eyd >> 5; /* approx month; might be 1 too low */
if (lt[res.hi + 1] <= eyd) /* fixup approximative month value */
res.hi += 1;
res.lo = eyd - lt[res.hi];
} else {
res.lo = res.hi = -1;
}
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. Since rd is signed, the remainder can
* be in the range [-6..+6], but the assignment to an unsigned
* variable maps the negative values to positive values >=7.
* This makes the sign correction look strange, but adding 7
* causes the needed wrap-around into the desired value range of
* zero to six, both inclusive.
*/
jd->weekday = rd % DAYSPERWEEK;
if (jd->weekday >= DAYSPERWEEK) /* weekday is unsigned! */
jd->weekday += DAYSPERWEEK;
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 = rd % DAYSPERWEEK;
if (utm->tm_wday < 0)
utm->tm_wday += DAYSPERWEEK;
/* 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
*
* ====================================================================
*/
/*
*---------------------------------------------------------------------
* 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
uint32_t p1, p2;
int isneg;
/*
* res = days *86400 + secs, using manual 16/32 bit
* multiplications and shifts.
*/
isneg = (days < 0);
if (isneg)
days = -days;
/* assemble days * 675 */
res.D_s.lo = (days & 0xFFFF) * 675u;
res.D_s.hi = 0;
p1 = (days >> 16) * 675u;
p2 = p1 >> 16;
p1 = p1 << 16;
M_ADD(res.D_s.hi, res.D_s.lo, p2, p1);
/* mul by 128, using shift */
res.D_s.hi = (res.D_s.hi << 7) | (res.D_s.lo >> 25);
res.D_s.lo = (res.D_s.lo << 7);
/* fix sign */
if (isneg)
M_NEG(res.D_s.hi, res.D_s.lo);
/* properly add seconds */
p2 = 0;
if (secs < 0) {
p1 = (uint32_t)-secs;
M_NEG(p2, p1);
} else {
p1 = (uint32_t)secs;
}
M_ADD(res.D_s.hi, res.D_s.lo, p2, p1);
# 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 sflag, sum, uyear;
sflag = int32_sflag(years);
uyear = int32_to_uint32_2cpl(years);
uyear ^= sflag;
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(sflag ^ 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 and correct if needed. (It likely is...) */
res.lo = m + 10;
res.hi = (res.lo >= 12);
if (res.hi)
res.lo -= 12;
/* if still out of range, normalise by floor division ... */
if (res.lo < 0 || res.lo >= 12) {
uint32_t mu, Q, sflag;
sflag = int32_sflag(res.lo);
mu = int32_to_uint32_2cpl(res.lo);
Q = sflag ^ ((sflag ^ mu) / 12u);
res.hi += uint32_2cpl_to_int32(Q);
res.lo = mu - Q * 12u;
}
/* get cummulated days in year with unshift */
res.lo = shift_month_table[res.lo] - 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) {
years += 1;
mdays += real_month_table[is_leapyear(years)][mons];
} 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 64-bit 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
* 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, sflag;
sflag = int32_sflag(years);
yu = int32_to_uint32_2cpl(years);
/* split off centuries, using floor division */
cc = sflag ^ ((sflag ^ 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(sflag ^ ((sflag ^ ci) / 4u));
ci = ci % 4u;
/* 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, sflag;
/* Use two fast cycle-split divisions here. This is again
* susceptible to internal overflow, so we check the range. This
* still permits more than +/-20 million years, so this is
* likely a pure academical problem.
*
* We want to execute '(weeks * 4 + 2) /% 20871' under floor
* division rules in the first step.
*/
sflag = int32_sflag(weeks);
sw = uint32_saturate(int32_to_uint32_2cpl(weeks), sflag);
sw = 4u * sw + 2;
Q = sflag ^ ((sflag ^ sw) / GREGORIAN_CYCLE_WEEKS);
sw -= Q * GREGORIAN_CYCLE_WEEKS;
ci = Q % 4u;
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; /* ws >> 13 , let the compiler sort it out */
sw = sw % 8192u; /* ws & 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, sflag;
/*
* 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 */
sflag = int32_sflag(ds.hi);
ud = int32_to_uint32_2cpl(ds.hi);
uw = sflag ^ ((sflag ^ 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 64-bit NTP timestamp from date/time.
*/
return isocal_date_to_ntp64(id).d_s.lo;
}
/* -*-EOF-*- */
