diff --git a/SConstruct b/SConstruct index 48ce498..97ad28e 100644 --- a/SConstruct +++ b/SConstruct @@ -81,13 +81,14 @@ env.Alias(target='install', source=env.InstallVersionedLib(dir=os.path.join(pref ## Unit Tests ... if test: - print('Compiling unit tests ') + print('Compiling unit test ...') tenv = env.Clone() tenv['CXXFLAGS'] = ' '.join([ x for x in env['CXXFLAGS'].split() if 'inline' not in x]) cmp_error_fn = 'test/unit_tests/test_compilation_error.json' cerror_dct = {} if os.path.isfile(cmp_error_fn): os.remove(cmp_error_fn) test_sources = glob.glob(r"test/unit_tests/*.cpp") + test_sources += glob.glob(r"test/precision/*.cpp") tenv.Append(RPATH=root_dir) for tsource in test_sources: ttarget = os.path.join(os.path.dirname(tsource), os.path.basename(tsource).replace('_', '-').replace('.cpp', '.out')) @@ -99,6 +100,7 @@ if test: 'flags': '{:}'.format(' '.join(['-o', tenv['CXXFLAGS']])), 'exit': 1} else: + # print('adding target {:}'.format(ttarget)) tenv.Program(target=ttarget, source=tsource, CPPPATH='src/', LIBS=vlib, LIBPATH='.') with open(cmp_error_fn, 'w') as fjson: print(json.dumps(cerror_dct, indent = 4), file=fjson) diff --git a/src/datetime_tops.hpp b/src/datetime_tops.hpp index 6ef70ab..9873df5 100644 --- a/src/datetime_tops.hpp +++ b/src/datetime_tops.hpp @@ -387,11 +387,20 @@ template #else template > #endif -double fractional_days(S nsec) noexcept { +double to_fractional_days(S nsec) noexcept { const double sec = static_cast(nsec.__member_ref__()); return sec / static_cast(S::max_in_day); } +#if __cplusplus >= 202002L +template +#else +template > +#endif +double to_fractional_seconds(S nsec) noexcept { + return nsec. S::template cast_to() * S::sec_inv_factor(); +} + /** Explicit cast of any second type to another second type. * * Cast an instance of any second type (aka any instance for which diff --git a/src/dtfund.hpp b/src/dtfund.hpp index caaebfb..9574374 100644 --- a/src/dtfund.hpp +++ b/src/dtfund.hpp @@ -1182,6 +1182,9 @@ class seconds { return static_cast(1e0); } + /** The scale factor to transform from seconds to seconds. */ + static constexpr double sec_inv_factor() noexcept { return 1e0; } + /** Constructor; default seconds is 0, but any integral will do */ explicit constexpr seconds(underlying_type i = 0) noexcept : m_sec(i){}; @@ -1249,11 +1252,6 @@ class seconds { } */ - /** Cast to double (i.e. fractional seconds). */ - constexpr double to_fractional_seconds() const noexcept { - return this->cast_to(); - } - private: /** Cast to any arithmetic type. */ #if __cplusplus >= 202002L @@ -1322,10 +1320,19 @@ class milliseconds { /** MilliSeconds are a subdivision of seconds. */ static constexpr bool is_of_sec_type = true; - /** The scale factor to transform from seconds to milliseconds. */ + /** The scale factor to transform from seconds to milliseconds. i.e. + * milliseconds = seconds * sec_factor() + */ template static constexpr T sec_factor() noexcept { return static_cast(1000); } + + /** The scale factor to transform from milliseconds to seconds, i.e. + * seconds = milliseconds * sec_inv_factor() + */ + static constexpr double sec_inv_factor() noexcept { + return 1e-3; + } /** Max milliseconds in one day. */ static constexpr underlying_type max_in_day = 86400L * 1000L; @@ -1367,11 +1374,6 @@ class milliseconds { return static_cast(m_sec) / static_cast(max_in_day); } - /** Cast to fractional dso::seconds */ - constexpr double to_fractional_seconds() const noexcept { - return static_cast(m_sec) / sec_factor(); - } - /** Cast to fractional hours */ constexpr double to_fractional_hours() const noexcept { constexpr const underlying_type secinh = @@ -1379,7 +1381,6 @@ class milliseconds { return static_cast(m_sec) / static_cast(secinh); } -private: /** Cast to any arithmetic type. */ #if __cplusplus >= 202002L template @@ -1392,6 +1393,7 @@ class milliseconds { return static_cast(m_sec); } +private: /** Milliseconds as underlying type. */ underlying_type m_sec; }; /* class milliseconds */ @@ -1455,10 +1457,19 @@ class microseconds { static constexpr underlying_type max_in_day{86'400L * 1'000'000L}; static_assert(max_in_day < std::numeric_limits::max()); - /** The scale factor to transform from seconds to microseconds. */ + /** The scale factor to transform from seconds to microseconds. i.e. + * microseconds = seconds * sec_factor() + */ template static constexpr T sec_factor() noexcept { return static_cast(1'000'000); } + + /** The scale factor to transform from microseconds to seconds, i.e. + * seconds = microseconds * sec_inv_factor() + */ + static constexpr double sec_inv_factor() noexcept { + return 1e-6; + } /** Constructor; default microseconds is 0; any integral number will do */ explicit constexpr microseconds(underlying_type i = 0L) noexcept : m_sec(i){}; @@ -1491,12 +1502,6 @@ class microseconds { return m_sec; } - /** Cast to fractional seconds */ - constexpr double to_fractional_seconds() const noexcept { - return static_cast(m_sec) / sec_factor(); - } - -private: /** Cast to any arithmetic type. */ #if __cplusplus >= 202002L template @@ -1509,6 +1514,7 @@ class microseconds { return static_cast(m_sec); } +private: /** Microseconds as long ints. */ underlying_type m_sec; }; /* microseconds */ @@ -1573,10 +1579,19 @@ class nanoseconds { static constexpr underlying_type max_in_day{86'400L * 1'000'000'000L}; static_assert(max_in_day < std::numeric_limits::max()); - /** The scale factor to transform from seconds to nanoseconds. */ + /** The scale factor to transform from seconds to nanoseconds. i.e. + * nanoseconds = seconds * sec_factor() + */ template static constexpr T sec_factor() noexcept { return static_cast(1'000'000'000); } + + /** The scale factor to transform from nanoseconds to seconds, i.e. + * seconds = nanoseconds * sec_inv_factor() + */ + static constexpr double sec_inv_factor() noexcept { + return 1e-9; + } /** Constructor; default nanoseconds is 0. **/ explicit constexpr nanoseconds(underlying_type i = 0L) noexcept : m_sec(i){}; @@ -1616,12 +1631,6 @@ class nanoseconds { return static_cast(m_sec) / static_cast(secinh); } - /** Cast to fractional seconds */ - constexpr double to_fractional_seconds() const noexcept { - return static_cast(m_sec) / sec_factor(); - } - -private: /** Cast to any arithmetic type. */ #if __cplusplus >= 202002L template @@ -1634,6 +1643,7 @@ class nanoseconds { return static_cast(m_sec); } +private: /** Nanoseconds as long ints. */ underlying_type m_sec; }; /* class nanoseconds */ diff --git a/src/tpdate.hpp b/src/tpdate.hpp index 50cba7f..a672296 100644 --- a/src/tpdate.hpp +++ b/src/tpdate.hpp @@ -78,8 +78,7 @@ class TwoPartDate { #endif TwoPartDate(const datetime &d) noexcept : _mjd((double)d.imjd().as_underlying_type()), - //_fday(d.sec().fractional_days()) { - _fday(fractional_days(d.sec())) { + _fday(to_fractional_days(d.sec())) { this->normalize(); } @@ -95,6 +94,15 @@ class TwoPartDate { /** Get the fractional part of the MJD*/ double fday() const noexcept { return _fday; } +#if __cplusplus >= 202002L + template +#else + template > +#endif + double sec_of_day() const noexcept { + return fday() * static_cast(T::max_in_day); + } + /** Add seconds to instance. * @warning Does not take into account leap seconds. */ diff --git a/test/precision/tp2dt.cpp b/test/precision/tp2dt.cpp new file mode 100644 index 0000000..3cc0157 --- /dev/null +++ b/test/precision/tp2dt.cpp @@ -0,0 +1,112 @@ +#include "dtcalendar.hpp" +#include "tpdate.hpp" +#include +#include + +using namespace dso; + +constexpr const long num_tests = 1'000'000; +using nsec = dso::nanoseconds; +typedef nsec::underlying_type SecIntType; +constexpr const SecIntType TOSEC = nsec::sec_factor(); +constexpr const double SECINDAY = nsec::max_in_day; + +int main() { + /* Generators for random numbers ... */ + std::random_device rd; + std::mt19937 gen(rd()); + std::uniform_int_distribution<> ydstr(1972, 2050); /* range for years */ + std::uniform_int_distribution<> mdstr(1, 12); /* range for months */ + std::uniform_int_distribution<> ddstr(1, 31); /* range for day of month */ + std::uniform_int_distribution nsdstr( + 0, nsec::max_in_day); /* range for day of month */ + + int testnr = 0, ok; + datetime d1, d2, d3, max_diff_date; + double max_diff_nsec = 0e0; + double max_diff_sec = 0e0; + while (testnr < num_tests) { + /* do we have a valid date ? */ + try { + d1 = datetime{year(ydstr(gen)), month(mdstr(gen)), + day_of_month(ddstr(gen)), nsec(nsdstr(gen))}; + d2 = datetime{year(ydstr(gen)), month(mdstr(gen)), + day_of_month(ddstr(gen)), nsec(nsdstr(gen))}; + /* d3 on same day as d2 */ + d3 = datetime{d2.imjd(), nsec(nsdstr(gen))}; + ok = 1; + } catch (std::exception &) { + ok = 0; + } + if (ok) { + /* construct a TwoPartDate from a datetime */ + TwoPartDate tpd1(d1); + assert(tpd1.imjd() == d1.imjd().as_underlying_type()); + /* fraction of day to nanoseconds */ + const double ns1 = tpd1.fday() * nsec::max_in_day; + /* difference from original nanoseconds */ + const double fd1 = ns1 - (double)d1.sec().as_underlying_type(); + TwoPartDate tpd2(d2); + assert(tpd2.imjd() == d2.imjd().as_underlying_type()); + const double ns2 = tpd2.fday() * nsec::max_in_day; + const double fd2 = ns2 - (double)d2.sec().as_underlying_type(); + TwoPartDate tpd3(d3); + assert(tpd3.imjd() == d3.imjd().as_underlying_type()); + const double ns3 = tpd3.fday() * nsec::max_in_day; + const double fd3 = ns3 - (double)d3.sec().as_underlying_type(); + + /* maximum difference */ + if (std::abs(fd1) > std::abs(max_diff_nsec)) { + max_diff_nsec = fd1; + max_diff_date = d1; + } + if (std::abs(fd2) > std::abs(max_diff_nsec)) { + max_diff_nsec = fd2; + max_diff_date = d2; + } + if (std::abs(fd3) > std::abs(max_diff_nsec)) { + max_diff_nsec = fd3; + max_diff_date = d3; + } + + /* do the same, but with seconds instead of nsec */ + const double fs1 = + tpd1.sec_of_day() - to_fractional_seconds(d1.sec()); + const double fs2 = + tpd2.sec_of_day() - to_fractional_seconds(d2.sec()); + const double fs3 = + tpd3.sec_of_day() - to_fractional_seconds(d3.sec()); + + /* maximum difference */ + if (std::abs(fs1) > std::abs(max_diff_sec)) { + max_diff_sec = fs1; + } + if (std::abs(fs2) > std::abs(max_diff_sec)) { + max_diff_sec = fs2; + } + if (std::abs(fs3) > std::abs(max_diff_sec)) { + max_diff_sec = fs3; + } + + ++testnr; + if (testnr % 10) + printf("%d/%ld\r", testnr, num_tests); + } + } + + printf("Maximum difference between a datetime and a TwoPartDate " + "instances is:\n%+.9e[nanosec] at datetime: %ld[MJD] %ld[nanosec]\n", + max_diff_nsec, max_diff_date.imjd().as_underlying_type(), + max_diff_date.sec().as_underlying_type()); + printf("Maximum difference between a datetime and a TwoPartDate " + "instances is:\n%+.9e[sec]\n", + max_diff_sec); + printf("When transforming a datetime to a TwoPartDate instance, " + "the worst case scenario is an accuracy of about %.1e[nanosec]\n", + max_diff_nsec); + printf("When transforming a datetime to a TwoPartDate instance, " + "the worst case scenario is an accuracy of about %.1e[sec]\n", + max_diff_sec); + + return 0; +} diff --git a/test/unit_tests/tpdates1.cpp b/test/unit_tests/tpdates1.cpp index 6ff3347..49d7d9e 100644 --- a/test/unit_tests/tpdates1.cpp +++ b/test/unit_tests/tpdates1.cpp @@ -10,6 +10,8 @@ using nsec = dso::nanoseconds; typedef nsec::underlying_type SecIntType; constexpr const SecIntType TOSEC = nsec::sec_factor(); constexpr const double SECINDAY = nsec::max_in_day; +constexpr const double PRECISION_SEC = 5e-12; +constexpr const double PRECISION_NSEC = 5e-3; int main() { /* Generators for random numbers ... */ @@ -39,34 +41,29 @@ int main() { if (ok) { /* construct a TwoPartDate from a datetime */ TwoPartDate tpd1(d1); + /* the difference between the two dates, in nsec should be below + * PRECISION_NSEC + */ assert(tpd1.imjd() == d1.imjd().as_underlying_type()); - assert(tpd1.fday() == d1.sec().to_fractional_seconds()/SECINDAY); - /* constructor from MJD and fraction of day */ - tpd1 = TwoPartDate(d1.imjd().as_underlying_type(), - d1.sec().to_fractional_seconds() / SECINDAY); - assert(tpd1.imjd() == d1.imjd().as_underlying_type()); - assert(tpd1.fday() == d1.sec().to_fractional_seconds() / SECINDAY); + assert(std::abs(tpd1.sec_of_day() - + d1.sec().as_underlying_type()) < + PRECISION_NSEC); + //const double dd1 = std::abs(tpd1.sec_of_day() - d1.sec().to_fractional_seconds()); TwoPartDate tpd2(d2); + assert(tpd2.imjd() == d2.imjd().as_underlying_type()); + assert(std::abs(tpd2.sec_of_day() - + d2.sec().as_underlying_type()) < + PRECISION_NSEC); + //const double dd2 = std::abs(tpd2.sec_of_day() - d2.sec().to_fractional_seconds()); + TwoPartDate tpd3(d3); + assert(tpd3.imjd() == d3.imjd().as_underlying_type()); + assert(std::abs(tpd3.sec_of_day() - + d3.sec().as_underlying_type()) < + PRECISION_NSEC); + //const double dd3 = std::abs(tpd3.sec_of_day() - d3.sec().to_fractional_seconds()); - /* how its done using datetime_intervals */ - const auto interval = d2 - d1; - auto d = d1 + interval; - assert(d == d2); - /* same thing using TwoPartDate's */ - const auto i1 = tpd2 - tpd1; - auto tpd = tpd1 + i1; - assert(tpd == tpd2); - - /* using datetime_interval ... */ - const auto d32 = d3 - d2; - d = d2 + d32; - assert(d == d3); - /* same thing using TwoPartDate's */ - const auto tpd32 = tpd3 - tpd2; - tpd = tpd2 + tpd32; - assert(tpd == tpd3); ++testnr; if (testnr % 10)