added precision

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)
 

src/datetime_tops.hpp

@@ -387,11 +387,20 @@ template <gconcepts::is_sec_dt S>
 #else
 template <typename S, typename = std::enable_if_t<S::is_of_sec_type>>
 #endif
-double fractional_days(S nsec) noexcept {
+double to_fractional_days(S nsec) noexcept {
   const double sec = static_cast<double>(nsec.__member_ref__());
   return sec / static_cast<double>(S::max_in_day);
 }
 
+#if __cplusplus >= 202002L
+template <gconcepts::is_sec_dt S>
+#else
+template <typename S, typename = std::enable_if_t<S::is_of_sec_type>>
+#endif
+double to_fractional_seconds(S nsec) noexcept {
+  return nsec. S::template cast_to<double>() * 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

src/dtfund.hpp

@@ -1182,6 +1182,9 @@ class seconds {
     return static_cast<T>(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<double>();
-  }
-
 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 <typename T> static constexpr T sec_factor() noexcept {
     return static_cast<T>(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,19 +1374,13 @@ class milliseconds {
     return static_cast<double>(m_sec) / static_cast<double>(max_in_day);
   }
 
-  /** Cast to fractional dso::seconds */
-  constexpr double to_fractional_seconds() const noexcept {
-    return static_cast<double>(m_sec) / sec_factor<double>();
-  }
-
   /** Cast to fractional hours */
   constexpr double to_fractional_hours() const noexcept {
     constexpr const underlying_type secinh =
         86400L * sec_factor<underlying_type>();
     return static_cast<double>(m_sec) / static_cast<double>(secinh);
   }
 
-private:
 /** Cast to any arithmetic type. */
 #if __cplusplus >= 202002L
   template <typename T>
@@ -1392,6 +1393,7 @@ class milliseconds {
     return static_cast<T>(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<underlying_type>::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 <typename T> static constexpr T sec_factor() noexcept {
     return static_cast<T>(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<double>(m_sec) / sec_factor<double>();
-  }
-
-private:
 /** Cast to any arithmetic type. */
 #if __cplusplus >= 202002L
   template <typename T>
@@ -1509,6 +1514,7 @@ class microseconds {
     return static_cast<T>(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<underlying_type>::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 <typename T> static constexpr T sec_factor() noexcept {
     return static_cast<T>(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<double>(m_sec) / static_cast<double>(secinh);
   }
 
-  /** Cast to fractional seconds */
-  constexpr double to_fractional_seconds() const noexcept {
-    return static_cast<double>(m_sec) / sec_factor<double>();
-  }
-
-private:
 /** Cast to any arithmetic type. */
 #if __cplusplus >= 202002L
   template <typename T>
@@ -1634,6 +1643,7 @@ class nanoseconds {
     return static_cast<T>(m_sec);
   }
 
+private:
   /** Nanoseconds as long ints. */
   underlying_type m_sec;
 }; /* class nanoseconds */

src/tpdate.hpp

@@ -78,8 +78,7 @@ class TwoPartDate {
 #endif
   TwoPartDate(const datetime<T> &d) noexcept
       : _mjd((double)d.imjd().as_underlying_type()),
-        //_fday(d.sec().fractional_days()) {
-        _fday(fractional_days<T>(d.sec())) {
+        _fday(to_fractional_days<T>(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 <gconcepts::is_sec_dt T>
+#else
+  template <typename T, typename = std::enable_if_t<T::is_of_sec_type>>
+#endif
+  double sec_of_day() const noexcept {
+    return fday() * static_cast<double>(T::max_in_day);
+  }
+
   /** Add seconds to instance.
    * @warning Does not take into account leap seconds.
    */

test/precision/tp2dt.cpp

@@ -0,0 +1,112 @@
+#include "dtcalendar.hpp"
+#include "tpdate.hpp"
+#include <cassert>
+#include <random>
+
+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<SecIntType>();
+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<SecIntType> nsdstr(
+      0, nsec::max_in_day); /* range for day of month */
+
+  int testnr = 0, ok;
+  datetime<nsec> 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<nsec>{year(ydstr(gen)), month(mdstr(gen)),
+                          day_of_month(ddstr(gen)), nsec(nsdstr(gen))};
+      d2 = datetime<nsec>{year(ydstr(gen)), month(mdstr(gen)),
+                          day_of_month(ddstr(gen)), nsec(nsdstr(gen))};
+      /* d3 on same day as d2 */
+      d3 = datetime<nsec>{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<seconds>() - to_fractional_seconds(d1.sec());
+      const double fs2 =
+          tpd2.sec_of_day<seconds>() - to_fractional_seconds(d2.sec());
+      const double fs3 =
+          tpd3.sec_of_day<seconds>() - 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<nanoseconds> 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<nanoseconds> and a TwoPartDate "
+         "instances is:\n%+.9e[sec]\n",
+         max_diff_sec);
+  printf("When transforming a datetime<nanoseconds> to a TwoPartDate instance, "
+         "the worst case scenario is an accuracy of about %.1e[nanosec]\n",
+         max_diff_nsec);
+  printf("When transforming a datetime<nanoseconds> to a TwoPartDate instance, "
+         "the worst case scenario is an accuracy of about %.1e[sec]\n",
+         max_diff_sec);
+
+  return 0;
+}