long double not working good

DSOlab · Oct 28, 2023 · c7e5cb4 · c7e5cb4
1 parent 332a6e9
commit c7e5cb4
Show file tree

Hide file tree

Showing 7 changed files with 136 additions and 40 deletions.
diff --git a/SConstruct b/SConstruct
@@ -89,6 +89,7 @@ if test:
     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")
+    test_sources += glob.glob(r"test/timer/*.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'))

diff --git a/src/datetime_tops.hpp b/src/datetime_tops.hpp
@@ -7,6 +7,7 @@
 #include "dtconcepts.hpp"
 #include <limits>
 #include <type_traits>
+#include <cstdio> // only for debugging
 
 namespace dso {
 
@@ -385,20 +386,26 @@ constexpr typename S::underlying_type max_days_allowed() {
 #if __cplusplus >= 202002L
 template <gconcepts::is_sec_dt S>
 #else
-template <typename S, typename = std::enable_if_t<S::is_of_sec_type>>
+template <typename S, typename T = double,
+          typename = std::enable_if_t<S::is_of_sec_type>,
+          typename = std::enable_if_t<std::is_floating_point<T>::value>>
 #endif
-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);
+T to_fractional_days(S nsec) noexcept {
+  printf("\tto_fractional_days got %ld\n", nsec.as_underlying_type());
+  const T sec = static_cast<T>(nsec.__member_ref__());
+  //return sec / static_cast<T>(S::max_in_day);
+  return sec / 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>>
+template <typename S, typename T = double,
+          typename = std::enable_if_t<S::is_of_sec_type>,
+          typename = std::enable_if_t<std::is_floating_point<T>::value>>
 #endif
-double to_fractional_seconds(S nsec) noexcept {
-  return nsec. S::template cast_to<double>() * S::sec_inv_factor();
+T to_fractional_seconds(S nsec) noexcept {
+  return nsec.S::template cast_to<T>() * S::sec_inv_factor();
 }
 
 /** Explicit cast of any second type to another second type.

diff --git a/src/datetime_write.hpp b/src/datetime_write.hpp
@@ -98,6 +98,7 @@ class SpitTime<S, HMSFormat::HHMMSSF> {
   static int spit(const hms_time<S> &hms, char *buffer) noexcept {
     /* seconds of minute (real) */
     double sec = to_fractional_seconds(hms.nsec());
+    printf("FUNCTION  %.15f\n", sec);
     return std::sprintf(
         buffer, "%02d:%02d:%012.9f", hms.hr().as_underlying_type(),
         hms.mn().as_underlying_type(), sec);

diff --git a/src/hms_time.hpp b/src/hms_time.hpp
@@ -45,7 +45,8 @@ class hms_time {
     assert(remaining < secInHour);
 #endif
     /* seconds in minute */
-    const SecIntType secInMin = 60L * S::template sec_factor<SecIntType>();
+    constexpr const SecIntType secInMin =
+        60L * S::template sec_factor<SecIntType>();
     /* compute/remove minutes */
     const SecIntType mn = remaining / secInMin;
     _minutes = mn;
@@ -55,9 +56,10 @@ class hms_time {
 #endif
     /* remaining S seconds */
     _sec = remaining;
+    printf("FUNCTION2 %15ld\n", _sec.as_underlying_type());
 #ifdef DEBUG
-    assert(_sec.as_underlying_type() + secInMin * minutes +
-               secInHour * _hours ==
+    assert(_sec.as_underlying_type() + secInMin * _minutes.as_underlying_type() +
+               secInHour * _hours.as_underlying_type() ==
            seconds.as_underlying_type());
 #endif
   }

diff --git a/src/tpdate.hpp b/src/tpdate.hpp
@@ -68,7 +68,7 @@ class TwoPartDate {
 private:
   using FDOUBLE = long double;
   double _mjd;  /** Mjd */
-  double _fday; /** fractional days */
+  FDOUBLE _fday; /** fractional days */
 
 public:
   /** Constructor from datetime<T> */
@@ -79,29 +79,29 @@ class TwoPartDate {
 #endif
   TwoPartDate(const datetime<T> &d) noexcept
       : _mjd((double)d.imjd().as_underlying_type()),
-        _fday(to_fractional_days<T>(d.sec())) {
+        _fday(to_fractional_days<T,FDOUBLE>(d.sec())) {
     this->normalize();
   }
 
   /** Constructor from a pair of doubles, such that MJD = a + b */
-  explicit constexpr TwoPartDate(double b = 0, double s = 0) noexcept
+  explicit TwoPartDate(double b = 0, FDOUBLE s = 0) noexcept
       : _mjd(b), _fday(s) {
     this->normalize();
   }
 
   /** Get the MJD as an intgral number, i.e. no fractional part */
   double imjd() const noexcept { return _mjd; }
 
-  /** Get the fractional part of the MJD*/
-  double fday() const noexcept { return _fday; }
-
+  /** Get the fractional part of the MJD */
+  FDOUBLE 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);
+  FDOUBLE sec_of_day() const noexcept {
+    return fday() * static_cast<FDOUBLE>(T::max_in_day);
   }
 
   /** @brief Transform the (integral part of the) date to Year Month Day */
@@ -110,13 +110,13 @@ class TwoPartDate {
   /** Add seconds to instance.
    * @warning Does not take into account leap seconds.
    */
-  void add_seconds(double sec) noexcept {
+  void add_seconds(FDOUBLE sec) noexcept {
     // first approach
     // _fday += sec / SEC_PER_DAY;
     // this->normalize();
     //
     // second approach
-    double dsec = _fday * SEC_PER_DAY;
+    FDOUBLE dsec = _fday * SEC_PER_DAY;
     dsec += sec;
     _fday = dsec / SEC_PER_DAY;
     this->normalize();
@@ -161,7 +161,7 @@ class TwoPartDate {
    * @warning Does not take into account leap seconds.
    */
   template <DateTimeDifferenceType DT>
-  double diff(const TwoPartDate &d) const noexcept {
+  FDOUBLE diff(const TwoPartDate &d) const noexcept {
     if constexpr (DT == DateTimeDifferenceType::FractionalDays) {
       /* difference as fractional days */
       return (_mjd - d._mjd) + (_fday - d._fday);
@@ -172,15 +172,15 @@ class TwoPartDate {
   }
 
   /** Get the date as (fractional) Julian Date */
-  double julian_date() const noexcept { return _fday + (_mjd + dso::MJD0_JD); }
+  FDOUBLE julian_date() const noexcept { return _fday + (_mjd + dso::MJD0_JD); }
 
   /** Transform instance to TT, assuming it is in TAI
    *
    * The two time scales are connected by the formula:
    * \f$ TT = TAI + ΔT \$ where \f$ ΔT = TT - TAI = 32.184 [sec] \f$
    */
   TwoPartDate tai2tt() const noexcept {
-    constexpr const double dtat = TT_MINUS_TAI / SEC_PER_DAY;
+    constexpr const FDOUBLE dtat = TT_MINUS_TAI / SEC_PER_DAY;
     return TwoPartDate(_mjd, _fday + dtat);
   }
 
@@ -190,14 +190,15 @@ class TwoPartDate {
    * \f$ UTC = TAI + ΔAT \$ where ΔAT are the leap seconds.
    */
   TwoPartDate utc2tai() const noexcept {
+    constexpr const FDOUBLE SPD = SEC_PER_DAY;
     auto utc = *this;
     /* Get TAI-UTC at 0h today and extra seconds in day (if any) */
     int extra;
     const int leap = dat(modified_julian_day((int)utc._mjd), extra);
     /* Remove any scaling applied to spread leap into preceding day */
-    utc._fday *= (SEC_PER_DAY + extra) / SEC_PER_DAY;
+    utc._fday *= (SPD + extra) / SPD;
     /* Assemble the TAI result, preserving the UTC split and order */
-    return TwoPartDate(utc._mjd, utc._fday + leap / SEC_PER_DAY);
+    return TwoPartDate(utc._mjd, utc._fday + leap / SPD);
   }
 
   /** Transform an instance to TT assuming it is in UTC */
@@ -210,8 +211,8 @@ class TwoPartDate {
   TwoPartDate tai2utc() const noexcept {
     // do it the SOFA way ...
     auto utc1(*this);
-    double small = utc1._fday;
-    double big = utc1._mjd;
+    FDOUBLE small = utc1._fday;
+    FDOUBLE big = utc1._mjd;
     for (int i = 0; i < 3; i++) {
       TwoPartDate guess = TwoPartDate(big, small).utc2tai();
       small += utc1._mjd - guess._mjd;
@@ -228,7 +229,8 @@ class TwoPartDate {
 
   /**  Transform an instance to TAI assuming it is in TT */
   TwoPartDate tt2tai() const noexcept {
-    constexpr const double dtat = TT_MINUS_TAI / SEC_PER_DAY;
+    constexpr const FDOUBLE SPD = SEC_PER_DAY;
+    constexpr const FDOUBLE dtat = TT_MINUS_TAI / SPD;
     return TwoPartDate(_mjd, _fday - dtat);
   }
 
@@ -242,14 +244,14 @@ class TwoPartDate {
    * @return The corresponding UT1 MJD, computed using:
    *         ∆T = TT − UT1 = 32.184[sec] + ∆AT − ∆UT1
    */
-  TwoPartDate tt2ut1(double dut1) const noexcept {
+  TwoPartDate tt2ut1(FDOUBLE dut1) const noexcept {
     /* note that ΔUT1 = UT1 − UTC hence UT1 = ΔUT1 + UTC */
     return tt2utc() + TwoPartDate(0e0, dut1 / SEC_PER_DAY);
   }
 
-  double as_mjd() const noexcept { return _fday + _mjd; }
+  FDOUBLE as_mjd() const noexcept { return _fday + _mjd; }
 
-  double jcenturies_sinceJ2000() const noexcept {
+  FDOUBLE jcenturies_sinceJ2000() const noexcept {
     return (_mjd - J2000_MJD) / DAYS_IN_JULIAN_CENT +
            _fday / DAYS_IN_JULIAN_CENT;
   }
@@ -267,29 +269,33 @@ class TwoPartDate {
     return (_mjd < d._mjd) || ((_mjd == d._mjd) && (_fday <= d._fday));
   }
 
-  /// @brief Keep _fday < 1e0 and _mjd integral only
-  constexpr void normalize() noexcept {
-    // fractional part should NOT be negative
+  /** @brief Normalize an instance. 
+   *
+   * Normalize here is meant in the sense that the fractional part of day is 
+   * stored in the _fday part, while _mjd holds the integral part of date
+   */
+  void normalize() noexcept {
+    /* fractional part should NOT be negative */
     while (_fday < 0e0) {
       _fday = 1 - _fday;
       _mjd -= 1e0;
     }
-    double fmore(0e0), extra(0e0);
-    // check if _mjd part has a fractional part
-    if ((fmore = std::modf(_mjd, &extra)) != 0e0) {
-      // assign fractional part to _fday and keep integral part to _mjd
+    FDOUBLE fmore(0e0), extra(0e0);
+    /* check if _mjd part has a fractional part */
+    if ((fmore = std::modf((FDOUBLE)_mjd, &extra)) != 0e0) {
+      /* assign fractional part to _fday and keep integral part to _mjd */
       _fday += fmore;
       _mjd = extra;
     }
-    // check if fractional part is >= 1e0
+    /* check if fractional part is >= 1e0 */
     if (_fday >= 1e0) {
       _fday = std::modf(_fday, &extra);
       _mjd += extra;
     }
 #ifdef DEBUG
     assert(_fday >= 0e0 && _fday < 1e0);
 #endif
-    // all done
+    /* all done */
     return;
   }
 

diff --git a/test/timer/test_tpdate_normalization.cpp b/test/timer/test_tpdate_normalization.cpp
@@ -0,0 +1,74 @@
+#include "dtcalendar.hpp"
+#include "tpdate.hpp"
+#include <cassert>
+#include <random>
+#include <vector>
+#include <chrono>
+#include <iostream>
+
+using namespace std::chrono;
+
+constexpr const long num_tests = 1'000'000;
+using nsec = dso::nanoseconds;
+using dso::datetime;
+typedef nsec::underlying_type SecIntType;
+using dso::TwoPartDate;
+using dso::SEC_PER_DAY;
+//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 */
+  std::uniform_real_distribution<double> rnds(-86400e0, 86400e0);
+
+  for (int Y=0; Y<5; Y++) {
+  auto start = high_resolution_clock::now();
+  int testnr = 0, ok, donotoptimize=0;
+  datetime<nsec> d1;
+  TwoPartDate tpd;
+  while (testnr < num_tests) {
+    /* do we have a valid date ? */
+    try {
+      d1 = datetime<nsec>{dso::year(ydstr(gen)), dso::month(mdstr(gen)),
+                          dso::day_of_month(ddstr(gen)), nsec(nsdstr(gen))};
+      ok = 1;
+    } catch (std::exception &) {
+      ok = 0;
+    }
+    if (ok) {
+      std::vector<double> dv;
+      /* construct a TwoPartDate from a datetime */
+      TwoPartDate tpd1(d1);
+      for (int i = 0; i < 10; i++) {
+        double s = rnds(gen);
+        /* do smthng studip to empty cache */
+        for (int j = 0; j < 10; j++) {
+          dv.push_back(s + j);
+        }
+        tpd1.add_seconds(s / SEC_PER_DAY);
+        for (int j = 0; j < 10; j++) {
+          if (dv[j] < 0e0) ++ok;
+        }
+        donotoptimize += ok - 100;
+        if (donotoptimize>1000) tpd = tpd1;
+      }
+    }
+    ++testnr;
+  }
+  auto stop = high_resolution_clock::now();
+  auto duration = duration_cast<microseconds>(stop - start);
+
+  std::cout << "Time spent on fucntion: " << duration.count() << "microsec\n";
+  printf("Here is smthng irrelevant, dummy=%d, mjd=%.1f\n", donotoptimize,
+         tpd.imjd());
+  }
+
+  return 0;
+}
diff --git a/test/unit_tests/dwrite9.cpp b/test/unit_tests/dwrite9.cpp
@@ -36,8 +36,13 @@ int main() {
   sz = dso::SpitDate<YMDFormat::YYYYMMDD>::numChars +
        dso::SpitTime<nanoseconds, HMSFormat::HHMMSSF>::numChars + 1;
   assert(!std::strncmp(buffer, "2023/10/24 00:00:00.000000059", sz));
+  printf("--------------------------------------------------\n");
   td1 = TwoPartDate(d1);
   to_char<YMDFormat::YYYYMMDD, HMSFormat::HHMMSSF>(td1, buffer);
+  printf("Expected: 2023/10/24 00:00:00.000000059\n");
+  printf("Got:      %.29s\n", buffer);
+  printf("Fr. Day:  %.15Lf\n", td1.sec_of_day<nanoseconds>());
+  printf("Fr. Day:  %.15Lf\n", td1.sec_of_day<dso::seconds>());
   assert(!std::strncmp(buffer, "2023/10/24 00:00:00.000000059", sz));
 
   d1 = datetime<nanoseconds>(year(2023), month(10), day_of_month(24),