diff --git a/include/fesa/analysis/thermal_time_grid.h b/include/fesa/analysis/thermal_time_grid.h index 8a16f32..6a6ac90 100644 --- a/include/fesa/analysis/thermal_time_grid.h +++ b/include/fesa/analysis/thermal_time_grid.h @@ -1,11 +1,44 @@ #ifndef FESA_ANALYSIS_THERMAL_TIME_GRID_H_ #define FESA_ANALYSIS_THERMAL_TIME_GRID_H_ + #include + #include "fesa/core/status.h" + namespace fesa { -class ThermalTimeGrid final { public: - static Result Create(double period,double increment,std::uint32_t maximum); - std::uint32_t Count() const noexcept; double Time(std::uint32_t index) const noexcept; double Delta(std::uint32_t index) const noexcept; - private: ThermalTimeGrid(double p,double d,std::uint32_t n,bool close):period_(p),increment_(d),count_(n),close_(close){} double period_,increment_; std::uint32_t count_; bool close_; }; -} -#endif + +/// @brief Carries one stored time and its actual integration coefficient. +struct ThermalTimeFrame { + double time; + double delta; +}; + +/// @brief Stores an O(1) fixed-increment grid with exact input-cell endpoint tests. +class ThermalTimeGrid final { + public: + /// @brief Checks finite positive controls and the uint32 INC bound before use. + static Result Create(double period, double increment, + std::uint32_t maximum = 100U); + /// @brief Returns the number of completed increments, excluding initial time. + std::uint32_t Count() const noexcept; + /// @brief Returns fl(index * nominal increment), or the exact final input time. + /// @pre index is in [0, Count()]. + double Time(std::uint32_t index) const noexcept; + /// @brief Returns the actual coefficient, never a difference of stored times. + /// @pre index is in [1, Count()]. + double Delta(std::uint32_t index) const noexcept; + /// @brief Checks index, finite values and strict increase for one frame in O(1). + Result Frame(std::uint32_t index) const; + + private: + /// @brief Stores controls only after count and terminal interval validation. + ThermalTimeGrid(double period, double increment, std::uint32_t count, + double final_delta); + double period_; + double increment_; + std::uint32_t count_; + double final_delta_; +}; + +} // namespace fesa +#endif // FESA_ANALYSIS_THERMAL_TIME_GRID_H_ diff --git a/src/fesa/analysis/thermal_time_grid.cpp b/src/fesa/analysis/thermal_time_grid.cpp index 2f14576..a7499ab 100644 --- a/src/fesa/analysis/thermal_time_grid.cpp +++ b/src/fesa/analysis/thermal_time_grid.cpp @@ -1,4 +1,222 @@ #include "fesa/analysis/thermal_time_grid.h" + +#include #include +#include +#include #include -namespace fesa { Result ThermalTimeGrid::Create(double p,double d,std::uint32_t max) { if(!std::isfinite(p)||!std::isfinite(d)||p<=0||d<=0||max==0) return Result::Failure(Status::Failure(FailureCategory::kModel,{{Severity::kError,"invalid-thermal-time-grid",{},"STEP","","Invalid thermal time grid."}})); if(d>=p) return Result::Success(ThermalTimeGrid{p,d,1,true}); const double q=p/d; const auto m=static_cast(std::nearbyint(q)); const bool close=m>0 && std::abs(std::fma(-static_cast(m),d,p)) <= .5*(std::nextafter(p,std::numeric_limits::infinity())-p)+.5*static_cast(m)*(std::nextafter(d,std::numeric_limits::infinity())-d); const std::uint64_t n=close?m:static_cast(std::ceil(q)); if(n>max||n>std::numeric_limits::max()) return Result::Failure(Status::Failure(FailureCategory::kModel,{{Severity::kError,"thermal-increment-limit",{},"STEP","","Thermal increment count exceeds INC."}})); return Result::Success(ThermalTimeGrid{p,d,static_cast(n),close}); } std::uint32_t ThermalTimeGrid::Count()const noexcept{return count_;} double ThermalTimeGrid::Time(std::uint32_t i)const noexcept{return i==count_?period_:static_cast(i)*increment_;} double ThermalTimeGrid::Delta(std::uint32_t i)const noexcept{return (i==count_&&!close_)?std::fma(-static_cast(count_-1),increment_,period_):increment_;} } + +namespace fesa { +namespace { + +// Fixed four base-2^32 digits. Contract bounds all exact values to 122 bits. +struct UInt128 { + std::array digit{}; +}; + +/// @brief Orders unsigned integers from their most significant digit. +int Compare(const UInt128& a, const UInt128& b) { + for (int i = 3; i >= 0; --i) { + if (a.digit[i] != b.digit[i]) return a.digit[i] < b.digit[i] ? -1 : 1; + } + return 0; +} + +/// @brief Adds with explicit carry; false means the fixed storage overflowed. +bool Add(const UInt128& a, const UInt128& b, UInt128& result) { + std::uint64_t carry = 0; + for (std::size_t i = 0; i < 4; ++i) { + const std::uint64_t sum = std::uint64_t{a.digit[i]} + b.digit[i] + carry; + result.digit[i] = static_cast(sum); + carry = sum >> 32; + } + return carry == 0; +} + +/// @brief Subtracts nonnegative magnitudes with explicit borrow. +bool Subtract(const UInt128& a, const UInt128& b, UInt128& result) { + std::uint64_t borrow = 0; + for (std::size_t i = 0; i < 4; ++i) { + const std::uint64_t sub = std::uint64_t{b.digit[i]} + borrow; + result.digit[i] = static_cast(std::uint64_t{a.digit[i]} - sub); + borrow = std::uint64_t{a.digit[i]} < sub ? 1U : 0U; + } + return borrow == 0; +} + +/// @brief Places a significand on the exact 2^-88 scaled lattice. +bool Shift(std::uint64_t value, int shift, UInt128& result) { + result = {}; + if (shift < 0) return false; + for (int bit = 0; bit < 64; ++bit) { + if ((value & (std::uint64_t{1} << bit)) == 0) continue; + const int target = bit + shift; + if (target >= 128) return false; + result.digit[target / 32] |= std::uint32_t{1} << (target % 32); + } + return true; +} + +/// @brief Multiplies by at most 2^32+1 using bounded schoolbook products. +bool Multiply(const UInt128& a, std::uint64_t multiplier, UInt128& result) { + result = {}; + for (std::size_t j = 0; j < 2; ++j) { + const auto factor = static_cast(multiplier >> (32 * j)); + std::uint64_t carry = 0; + for (std::size_t i = 0; i < 4; ++i) { + if (i + j >= 4) { + if (a.digit[i] != 0 && factor != 0) return false; + continue; + } + const std::uint64_t product = std::uint64_t{a.digit[i]} * factor + + result.digit[i + j] + carry; + result.digit[i + j] = static_cast(product); + carry = product >> 32; + } + if (carry != 0) return false; + } + return true; +} + +struct Cell { + UInt128 center; + UInt128 lower; + UInt128 upper; + bool closed; +}; + +/// @brief Constructs ORIGINAL binary64 cells, including subnormal/maxfinite gaps. +bool MakeCell(double value, int period_exponent, Cell& cell) { + std::uint64_t bits = 0; + std::memcpy(&bits, &value, sizeof(bits)); + const int exponent_bits = static_cast((bits >> 52) & 0x7ffU); + const std::uint64_t fraction = bits & ((std::uint64_t{1} << 52) - 1); + const std::uint64_t significand = + exponent_bits == 0 ? fraction : fraction + (std::uint64_t{1} << 52); + const int exponent = exponent_bits == 0 ? -1074 : exponent_bits - 1075; + const int lattice_shift = exponent - period_exponent + 88; + UInt128 lower_gap; + UInt128 upper_gap; + // At a normal binade boundary only the predecessor gap is halved. The + // smallest normal has subnormal neighbors with the same original spacing. + const int boundary = fraction == 0 && exponent_bits > 1 ? 1 : 0; + cell.closed = (bits & 1) == 0; + return Shift(significand, lattice_shift, cell.center) && + Shift(1, lattice_shift - 1 - boundary, lower_gap) && + Shift(1, lattice_shift - 1, upper_gap) && + Subtract(cell.center, lower_gap, cell.lower) && + Add(cell.center, upper_gap, cell.upper); +} + +/// @brief Returns source-neutral utility failure for parser/procedure decoration. +Status GridFailure(const char* code, const char* message) { + return Status::Failure(FailureCategory::kModel, + {{Severity::kError, code, {}, "HEAT TRANSFER", "", + message}}); +} + +} // namespace + +Result ThermalTimeGrid::Create( + double period, double increment, std::uint32_t maximum) { + if (!std::isfinite(period) || !std::isfinite(increment) || period <= 0 || + increment <= 0 || maximum == 0) { + return Result::Failure(GridFailure( + "invalid-thermal-time-grid", "Time controls must be finite and positive.")); + } + if (increment >= period) { + return Result::Success( + ThermalTimeGrid{period, increment, 1U, period}); + } + const double quotient = period / increment; + const std::uint64_t limit = std::uint64_t{maximum} + 1U; + if (!std::isfinite(quotient) || quotient > static_cast(limit)) { + return Result::Failure(GridFailure( + "thermal-increment-limit", "Thermal increment count exceeds INC.")); + } + // No floating-to-integer cast precedes the finite, positive count bound. + std::uint64_t floor = static_cast(std::floor(quotient)); + const double fraction = quotient - static_cast(floor); + const std::uint64_t nearest = + floor + ((fraction > .5 || (fraction == .5 && (floor & 1) != 0)) ? 1 : 0); + int exponent = 0; + std::frexp(period, &exponent); + Cell period_cell{}; + Cell increment_cell{}; + UInt128 lower; + UInt128 upper; + bool exact_ok = MakeCell(period, exponent, period_cell) && + MakeCell(increment, exponent, increment_cell) && + Multiply(increment_cell.lower, nearest, lower) && + Multiply(increment_cell.upper, nearest, upper); + const int first = Compare(period_cell.upper, lower); + const int second = Compare(upper, period_cell.lower); + const bool closed = period_cell.closed && increment_cell.closed; + const bool close = nearest != 0 && + (first > 0 || (first == 0 && closed)) && + (second > 0 || (second == 0 && closed)); + std::uint64_t count = nearest; + if (!close && exact_ok) { + UInt128 product; + UInt128 next; + exact_ok = Multiply(increment_cell.center, floor, product) && + Multiply(increment_cell.center, floor + 1, next); + if (Compare(product, period_cell.center) > 0) { + --floor; + } else if (Compare(next, period_cell.center) <= 0) { + ++floor; + } + exact_ok = exact_ok && Multiply(increment_cell.center, floor, product) && + Multiply(increment_cell.center, floor + 1, next) && + Compare(product, period_cell.center) <= 0 && + Compare(next, period_cell.center) > 0; + count = floor + (Compare(product, period_cell.center) == 0 ? 0U : 1U); + } + if (!exact_ok) { + return Result::Failure(GridFailure( + "invalid-thermal-time-grid", "Exact time lattice bounds failed.")); + } + if (count == 0 || count > maximum) { + return Result::Failure(GridFailure( + "thermal-increment-limit", "Thermal increment count exceeds INC.")); + } + const double last = close ? increment : + std::fma(-static_cast(count - 1), increment, period); + ThermalTimeGrid grid{period, increment, static_cast(count), last}; + const auto frame = grid.Frame(grid.Count()); + if (!frame.HasValue()) return Result::Failure(frame.GetStatus()); + // For j <= 2^32, adjacent exact j*d are separated by more than 2^19 ULPs + // (subnormal products are exact). Thus the terminal check also suffices for + // O(1) preflight; every requested frame still checks its own stored interval. + return Result::Success(grid); +} + +ThermalTimeGrid::ThermalTimeGrid(double period, double increment, + std::uint32_t count, double final_delta) + : period_{period}, increment_{increment}, count_{count}, + final_delta_{final_delta} {} + +std::uint32_t ThermalTimeGrid::Count() const noexcept { return count_; } +double ThermalTimeGrid::Time(std::uint32_t index) const noexcept { + return index == count_ ? period_ : static_cast(index) * increment_; +} +double ThermalTimeGrid::Delta(std::uint32_t index) const noexcept { + return index == count_ ? final_delta_ : increment_; +} +Result ThermalTimeGrid::Frame(std::uint32_t index) const { + if (index == 0 || index > count_) { + return Result::Failure(GridFailure( + "invalid-thermal-time-frame", "Completed frame index is outside the grid.")); + } + const double time = Time(index); + const double delta = Delta(index); + const double previous = Time(index - 1); + if (!std::isfinite(time) || !std::isfinite(previous) || !std::isfinite(delta) || + !(time > previous) || !(delta > 0)) { + return Result::Failure(GridFailure( + "invalid-thermal-time-frame", "Stored times must be finite and strictly increasing.")); + } + return Result::Success({time, delta}); +} + +} // namespace fesa diff --git a/tests/unit/analysis/thermal_time_grid_test.cpp b/tests/unit/analysis/thermal_time_grid_test.cpp new file mode 100644 index 0000000..59448ee --- /dev/null +++ b/tests/unit/analysis/thermal_time_grid_test.cpp @@ -0,0 +1,103 @@ +#include "fesa/analysis/thermal_time_grid.h" + +#include + +#include +#include + +TEST(ThermalTimeGrid, SingleStepUsesPeriodAsActualIncrement) { + const auto grid = fesa::ThermalTimeGrid::Create(2.0, 3.0, 1U); + ASSERT_TRUE(grid.HasValue()); + EXPECT_EQ(grid.Value().Count(), 1U); + EXPECT_EQ(grid.Value().Delta(1U), 2.0); +} + +TEST(ThermalTimeGrid, PreservesGenuineRemainderOutsideInputCells) { + const auto grid = fesa::ThermalTimeGrid::Create( + std::nextafter(.9, 1.0), .3, 4U); + ASSERT_TRUE(grid.HasValue()); + EXPECT_EQ(grid.Value().Count(), 4U); + EXPECT_EQ(grid.Value().Delta(4U), 1.6653345369377348e-16); +} + +TEST(ThermalTimeGrid, OriginalSubnormalCellsRespectOpenAndClosedContact) { + const double tiny = std::numeric_limits::denorm_min(); + const auto even = fesa::ThermalTimeGrid::Create(20 * tiny, 6 * tiny, 4U); + const auto odd = fesa::ThermalTimeGrid::Create(17 * tiny, 5 * tiny, 4U); + ASSERT_TRUE(even.HasValue()); + ASSERT_TRUE(odd.HasValue()); + EXPECT_EQ(even.Value().Count(), 3U); + EXPECT_EQ(odd.Value().Count(), 4U); + EXPECT_EQ(odd.Value().Delta(4U), 2 * tiny); +} + +TEST(ThermalTimeGrid, DecimalCellsCounterexampleAndScale) { + for (const int exponent : {-1000, 0, 1000}) { + const double period = std::ldexp(.9, exponent); + const double delta = std::ldexp(.3, exponent); + const auto grid = fesa::ThermalTimeGrid::Create(period, delta, 3U); + ASSERT_TRUE(grid.HasValue()); + EXPECT_EQ(grid.Value().Count(), 3U); + EXPECT_EQ(grid.Value().Delta(3U), delta); + EXPECT_EQ(grid.Value().Time(3U), period); + } + const auto old_counterexample = fesa::ThermalTimeGrid::Create( + 5.351455173733616, .17262758624947144, 31U); + ASSERT_TRUE(old_counterexample.HasValue()); + EXPECT_EQ(old_counterexample.Value().Count(), 31U); + EXPECT_EQ(old_counterexample.Value().Delta(31U), .17262758624947144); + EXPECT_FALSE(fesa::ThermalTimeGrid::Create( + 5.351455173733616, .17262758624947144, 30U).HasValue()); +} + +TEST(ThermalTimeGrid, BinadeLimitsAndMaximumFinitePreserveOriginalCells) { + const double max = std::numeric_limits::max(); + const auto grid = fesa::ThermalTimeGrid::Create(max, max / 3, 3U); + ASSERT_TRUE(grid.HasValue()); + EXPECT_EQ(grid.Value().Count(), 3U); + EXPECT_EQ(grid.Value().Delta(3U), max / 3); + const auto upper = fesa::ThermalTimeGrid::Create(4294967295., 1., 4294967295U); + ASSERT_TRUE(upper.HasValue()); + EXPECT_EQ(upper.Value().Count(), 4294967295U); + EXPECT_TRUE(upper.Value().Frame(4294967295U).HasValue()); + EXPECT_FALSE(fesa::ThermalTimeGrid::Create(4294967296., 1., 4294967295U).HasValue()); + EXPECT_FALSE(fesa::ThermalTimeGrid::Create(max, + std::numeric_limits::denorm_min(), 4294967295U).HasValue()); +} + +TEST(ThermalTimeGrid, ShortenedIncrementAndCheckedFrames) { + const auto grid = fesa::ThermalTimeGrid::Create(10., 3., 4U); + ASSERT_TRUE(grid.HasValue()); + EXPECT_EQ(grid.Value().Count(), 4U); + EXPECT_EQ(grid.Value().Time(0U), 0.); + for (std::uint32_t i = 1; i <= 4; ++i) { + ASSERT_TRUE(grid.Value().Frame(i).HasValue()); + EXPECT_EQ(grid.Value().Time(i), i == 4 ? 10. : 3. * i); + EXPECT_EQ(grid.Value().Delta(i), i == 4 ? 1. : 3.); + } + EXPECT_FALSE(grid.Value().Frame(0U).HasValue()); + EXPECT_FALSE(grid.Value().Frame(5U).HasValue()); + const auto decimals = fesa::ThermalTimeGrid::Create(1., .1); + ASSERT_TRUE(decimals.HasValue()); + EXPECT_EQ(decimals.Value().Time(3U), .30000000000000004); + EXPECT_EQ(decimals.Value().Delta(10U), .1); + const double below = std::nextafter(.9, 0.); + const auto close = fesa::ThermalTimeGrid::Create(below, .3, 3U); + const auto short_last = fesa::ThermalTimeGrid::Create(std::nextafter(below, 0.), .3, 3U); + ASSERT_TRUE(close.HasValue()); + ASSERT_TRUE(short_last.HasValue()); + EXPECT_EQ(close.Value().Delta(3U), .3); + EXPECT_EQ(short_last.Value().Delta(3U), .2999999999999998); +} + +TEST(ThermalTimeGrid, RejectsInvalidControlsBeforeIntegerConversion) { + for (double bad : {0., -1., std::numeric_limits::infinity(), + std::numeric_limits::quiet_NaN()}) { + EXPECT_FALSE(fesa::ThermalTimeGrid::Create(bad, 1.).HasValue()); + EXPECT_FALSE(fesa::ThermalTimeGrid::Create(1., bad).HasValue()); + } + EXPECT_FALSE(fesa::ThermalTimeGrid::Create(1., 1., 0U).HasValue()); + const auto equal = fesa::ThermalTimeGrid::Create(2., 2., 1U); + ASSERT_TRUE(equal.HasValue()); + EXPECT_EQ(equal.Value().Delta(1U), 2.); +}