diff --git a/include/absi/polar.h b/include/absi/polar.h index 2aa0c11..e6aa2a3 100644 --- a/include/absi/polar.h +++ b/include/absi/polar.h @@ -43,21 +43,23 @@ namespace polar [[nodiscard]] AngleInterval operator+(const AngleInterval &other) const noexcept; - [[nodiscard]] AngleInterval operator*(const AngleInterval &other) const noexcept; - [[nodiscard]] AngleInterval operator|(const AngleInterval &other) const noexcept; [[nodiscard]] bool operator==(const AngleInterval &other) const noexcept; + Real min() const noexcept; + Real delta() const noexcept; + Real max() const noexcept; + protected: /// @brief The minimal value of the interval /// `0 < @ref min < 2` is expressed in radians to limit - Real min; + Real _min; /// @brief The wideness/uncertainty of the interval /// Values in the interval are between @ref min and @ref min + @ref delta /// Moreover, 0 <= @ref delta <= 2 - Real delta; + Real _delta; void set_remainder() noexcept; }; @@ -107,4 +109,4 @@ namespace polar const Interval zero = 0.; const Interval one = 1.; -} +}; diff --git a/src/absi/polar.cpp b/src/absi/polar.cpp index 80eff12..36ec56e 100644 --- a/src/absi/polar.cpp +++ b/src/absi/polar.cpp @@ -42,73 +42,279 @@ bool polar::PositiveInterval::operator==(const PositiveInterval &other) const return (min == other.min) && (max == other.max); } -AngleInterval::AngleInterval(Real min, Real delta) : min(min), delta(delta) -{ - if (delta < 0) - { - throw std::range_error("Negative delta for angle interval"); +static constexpr Real TWO = 2.0; +static constexpr Real EPS = 1e-12; + +static inline Real mod2(Real x) noexcept { + x = std::fmod(x, TWO); + if (x < 0) x += TWO; + // prefer 0 <= x < 2 + if (x >= TWO - EPS) x = 0.0; + return x; +} + +AngleInterval::AngleInterval(Real min_, Real delta_) : _min(mod2(min_)), _delta(delta_) { + // clamp delta in [0,2] + if (_delta < 0) _delta = 0; + if (_delta > TWO) _delta = TWO; + if (_delta >= TWO - EPS) { + // full circle canonical representation: start at 0, delta = 2 + _min = 0; + _delta = TWO; + } } - if (min > 2 || min < 0) - { - int ratio = min / 2; - min = min - 2 * ratio; + +AngleInterval::AngleInterval(Real a) : AngleInterval(a, 0.0) {} + + AngleInterval AngleInterval::min_max(Real a, Real b) noexcept { + // helper if you need it; not used in operator| below + if (a <= b) return AngleInterval(a, b - a); + return AngleInterval(a, b - a + TWO); } - if (delta > 2) - { - delta = 2; + + + [[nodiscard]] AngleInterval AngleInterval::operator+(const AngleInterval &other) const noexcept { + // if either is full circle -> sum is full + if (_delta >= TWO - EPS || other._delta >= TWO - EPS) { + return AngleInterval(0.0, TWO); + } + + // Convert possibly-wrapping interval into 0..2 non-wrapping segments (1 or 2 segments) + auto linear_segments = [](Real s, Real d) { + std::vector> segs; + if (d <= 0.0 + EPS) { + // degenerate point + segs.emplace_back(mod2(s), mod2(s)); + return segs; + } + Real e = s + d; + if (e <= TWO + EPS) { + segs.emplace_back(mod2(s), std::min(e, TWO)); + } else { + // wrap + segs.emplace_back(mod2(s), TWO); + segs.emplace_back(0.0, e - TWO); + } + return segs; + }; + + auto segsA = linear_segments(_min, _delta); + auto segsB = linear_segments(other._min, other._delta); + + // Build summed segments (on the line) for every pair + std::vector> summed_mod_segments; + + for (auto &a : segsA) { + Real a_len = a.second - a.first; + for (auto &b : segsB) { + Real b_len = b.second - b.first; + Real sum_len = a_len + b_len; + + // If any pair-sum length >= full circle => entire circle covered + if (sum_len >= TWO - EPS) { + return AngleInterval(0.0, TWO); + } + + // pairwise sum on real line + Real s = a.first + b.first; // in [0,4) + Real e = a.second + b.second; // s + sum_len + + // map [s,e) into [0,2) possibly splitting + // three cases: e <= 2 -> single segment [s,e] + // : s >= 2 -> map down by -2: [s-2, e-2] + // : s < 2 && e > 2 -> split [s,2) and [0, e-2) + if (e <= TWO + EPS) { + // no wrap after sum + summed_mod_segments.emplace_back(mod2(s), std::min(e, TWO)); + } else if (s >= TWO - EPS) { + // entirely >= 2, shift down + summed_mod_segments.emplace_back(mod2(s - TWO), mod2(e - TWO)); + } else { + // crosses the 2 boundary -> produces two segments + summed_mod_segments.emplace_back(mod2(s), TWO); + summed_mod_segments.emplace_back(0.0, e - TWO); + } + } + } + + if (summed_mod_segments.empty()) { + // Both degenerate with zero length but code above should have created entries. + return AngleInterval(0.0, 0.0); + } + + // Sort and merge overlapping/adjacent segments on [0,2) + std::sort(summed_mod_segments.begin(), summed_mod_segments.end(), + [](auto &A, auto &B) { + if (std::fabs(A.first - B.first) > EPS) return A.first < B.first; + return A.second < B.second; + }); + + std::vector> merged; + for (auto &seg : summed_mod_segments) { + if (merged.empty()) merged.push_back(seg); + else { + auto &last = merged.back(); + if (seg.first <= last.second + EPS) { + last.second = std::max(last.second, seg.second); + } else { + merged.push_back(seg); + } + } + } + + // Compute covered length + Real covered = 0.0; + for (auto &m : merged) covered += (m.second - m.first); + + if (covered >= TWO - EPS) { + return AngleInterval(0.0, TWO); + } + + // If single merged segment -> done + if (merged.size() == 1) { + Real s = merged[0].first; + Real len = merged[0].second - merged[0].first; + return AngleInterval(s, len); + } + + // Otherwise find largest gap (including wrap gap) and take complement + Real max_gap = -1.0; + Real gap_start_pos = 0.0; + for (size_t i = 0; i < merged.size(); ++i) { + Real end_i = merged[i].second; + Real start_next = (i + 1 < merged.size()) ? merged[i+1].first : (merged[0].first + TWO); + Real gap = start_next - end_i; + if (gap > max_gap) { + max_gap = gap; + // start of the minimal covering arc is the beginning of the next segment + gap_start_pos = mod2(start_next); + } + } + + Real cover_len = TWO - max_gap; + Real start_of_cover = mod2(gap_start_pos); + return AngleInterval(start_of_cover, cover_len); } -} -AngleInterval::AngleInterval(Real min) : min(min), delta(0.) {} + [[nodiscard]] AngleInterval AngleInterval::operator|(const AngleInterval &other) const noexcept { + // if either is full circle, return full + if (_delta >= TWO - EPS || other._delta >= TWO - EPS) { + return AngleInterval(0.0, TWO); + } -AngleInterval polar::AngleInterval::min_max(Real a, Real b) noexcept -{ - return AngleInterval(std::min(a, b), std::abs(a - b)); -} + // convert (possibly-wrapping) intervals into non-wrapping segments on [0,2) + std::vector> segs; + auto add_segments = [&](Real s, Real d) { + if (d <= 0.0 + EPS) { + // zero-length interval: represent as small degenerate segment [s, s] + segs.emplace_back(s, s); + return; + } + Real e = s + d; + if (e <= TWO + EPS) { + // no wrap + segs.emplace_back(s, std::min(e, TWO)); + } else { + // wrap: [s,2) and [0, e-2) + segs.emplace_back(s, TWO); + segs.emplace_back(0.0, e - TWO); + } + }; -AngleInterval AngleInterval::operator+(const AngleInterval &other) const noexcept -{ - auto r = AngleInterval(min + other.min, delta + other.delta); - r.set_remainder(); - return r; -} + add_segments(_min, _delta); + add_segments(other._min, other._delta); -AngleInterval AngleInterval::operator*(const AngleInterval &other) const noexcept -{ - auto r = AngleInterval(std::min(min * other.min, (min + delta) * (other.min + other.delta)), delta * other.delta); - r.set_remainder(); - return r; -} + // sort by start + std::sort(segs.begin(), segs.end(), [](auto &a, auto &b){ + if (std::fabs(a.first - b.first) > EPS) return a.first < b.first; + return a.second < b.second; + }); -AngleInterval polar::AngleInterval::operator|(const AngleInterval &other) const noexcept -{ - return AngleInterval::min_max(std::min(min, other.min), std::max(min + delta, other.min + other.delta)); -} + // merge overlapping/adjacent segments + std::vector> merged; + for (auto &seg : segs) { + if (merged.empty()) { + merged.push_back(seg); + } else { + auto &last = merged.back(); + if (seg.first <= last.second + EPS) { + // overlap or touch + last.second = std::max(last.second, seg.second); + } else { + merged.push_back(seg); + } + } + } -bool polar::AngleInterval::operator==(const AngleInterval &other) const noexcept -{ - return (min == other.min) && (delta == other.delta); -} + // compute covered length + Real covered = 0.0; + for (auto &m : merged) covered += (m.second - m.first); -void AngleInterval::set_remainder() noexcept -{ - if (min > 2 || min < 0) - { - int ratio = min / 2; - min = min - 2 * ratio; + // if coverage is entire circle (within epsilon) -> full + if (covered >= TWO - EPS) { + return AngleInterval(0.0, TWO); + } + + // If single merged segment -> return it directly + if (merged.size() == 1) { + Real s = merged[0].first; + Real len = merged[0].second - merged[0].first; + return AngleInterval(s, len); + } + + // find largest gap between merged segments (including wrap gap) + Real max_gap = -1.0; + Real gap_start_pos = 0.0; // start position of the minimal covering arc + for (size_t i = 0; i < merged.size(); ++i) { + Real end_i = merged[i].second; + Real start_next; + if (i + 1 < merged.size()) start_next = merged[i+1].first; + else start_next = merged[0].first + TWO; // wrap-around + Real gap = start_next - end_i; // positive + if (gap > max_gap) { + max_gap = gap; + // the minimal covering arc should start at the beginning of the next segment + gap_start_pos = mod2(start_next); + } + } + + // Complement of largest gap is the minimal covering arc + Real cover_len = TWO - max_gap; + Real start_of_cover = mod2(gap_start_pos); + + return AngleInterval(start_of_cover, cover_len); } - if (delta > 2) - { - delta = 2; + + [[nodiscard]] bool AngleInterval::operator==(const AngleInterval &other) const noexcept { + // normalize comparison: if full circle + if (_delta >= TWO - EPS && other._delta >= TWO - EPS) return true; + return (std::fabs(mod2(_min) - mod2(other._min)) < 1e-10) && + (std::fabs(_delta - other._delta) < 1e-10); + } + + // utility for testing & debugging + Real AngleInterval::min() const noexcept { return _min; } + Real AngleInterval::delta() const noexcept { return _delta; } + Real AngleInterval::max() const noexcept { + auto m = _min + _delta; + if (m > TWO) m -= TWO; + return m; + } + + void AngleInterval::set_remainder() noexcept { + _min = mod2(_min); + if (_delta >= TWO - EPS) { + _min = 0.0; + _delta = TWO; + } } -} Interval::Interval(PositiveInterval mod, AngleInterval arg) : mod(mod), arg(arg) {} polar::Interval::Interval() : mod(PositiveInterval(0.)), arg(AngleInterval(0.)) {}; polar::Interval::Interval(const polar::Real value) : mod(PositiveInterval(std::abs(value))), - arg(AngleInterval(value < 0. ? 1. : 0.)) {} + arg(AngleInterval(value < 0. ? 1. : 0.)) {} polar::Interval::Interval(const ampl::Amplitude z) : mod(PositiveInterval(std::abs(z))), arg(argument(z)) {} @@ -173,17 +379,17 @@ std::string polar::Interval::to_string(bool strict) const noexcept if (is_real()) return std::to_string(to_real()); } - return "{mod: " + std::to_string(mod.min) + " " + std::to_string(mod.max) + " arg: " + std::to_string(arg.min) + " " + std::to_string(arg.delta) + "}"; + return "{mod: " + std::to_string(mod.min) + " " + std::to_string(mod.max) + " arg: " + std::to_string(arg._min) + " " + std::to_string(arg._delta) + "}"; } polar::Real polar::Interval::norm() { - return arg.delta * mod.max; // TODO: better approximation + return arg._delta * mod.max / 2; // TODO: better approximation } bool polar::Interval::is_real() const noexcept { - return (arg.delta == 0) && (arg.min == 0. || arg.min == 1.); + return (arg._delta == 0) && (arg._min == 0. || arg._min == 1.); } polar::Real polar::Interval::to_real() const @@ -192,5 +398,5 @@ polar::Real polar::Interval::to_real() const { throw std::logic_error("Not a real number"); } - return mod.min * ((arg.min == 0.) ? 1. : -1.); + return mod.min * ((arg._min == 0.) ? 1. : -1.); } \ No newline at end of file diff --git a/test/absi_test.cpp b/test/absi_test.cpp index 40015fc..354652c 100644 --- a/test/absi_test.cpp +++ b/test/absi_test.cpp @@ -1,6 +1,8 @@ #include #include +using polar::PositiveInterval, polar::AngleInterval; + class CartesianTest : public testing::Test { public: @@ -8,13 +10,13 @@ class CartesianTest : public testing::Test cartesian::Interval b = cartesian::Interval(ampl::zero, ampl::one + ampl::i); }; -TEST_F(CartesianTest, Norm) +TEST_F(CartesianTest, norm) { EXPECT_EQ(a.norm(), ampl::zero_real); EXPECT_EQ(b.norm(), ampl::sqrt2); } -TEST_F(CartesianTest, Contains) +TEST_F(CartesianTest, contains) { EXPECT_TRUE(a.contains(ampl::one)); EXPECT_FALSE(a.contains(ampl::i)); @@ -22,7 +24,7 @@ TEST_F(CartesianTest, Contains) EXPECT_FALSE(b.contains(ampl::i * 3.)); } -TEST_F(CartesianTest, Equal) +TEST_F(CartesianTest, equal) { cartesian::Interval c = cartesian::Interval(ampl::zero, ampl::one + ampl::i); EXPECT_EQ(b, c); @@ -30,7 +32,7 @@ TEST_F(CartesianTest, Equal) EXPECT_EQ(ampl::zero, ampl::one * ampl::zero); } -TEST_F(CartesianTest, Join) +TEST_F(CartesianTest, join) { auto j = cartesian::zero | cartesian::Interval(ampl::i + ampl::one); @@ -40,32 +42,31 @@ TEST_F(CartesianTest, Join) class PolarTest : public testing::Test { public: - polar::PositiveInterval pa = polar::PositiveInterval(5.); - polar::PositiveInterval pb = polar::PositiveInterval(3., 4.); - polar::AngleInterval ra = polar::AngleInterval(-.25, .2); - polar::AngleInterval rb = polar::AngleInterval(3); - polar::Interval i1 = polar::Interval(pa, ra); - polar::Interval i2 = polar::Interval(pb, rb); + polar::PositiveInterval pa{5.}, pb{3., 4.}, pa_plus_pb{8., 9.}, pa_times_pb{15., 20.}, pa_join_pb{3., 5.}; + polar::AngleInterval ra{-.25, 1.25}, ra_other{1.75, 1.25}, rb{1}, ra_plus_rb{.75, 1.25}, ra_join_rb{ra}; + polar::Interval i1{pa, ra}, i2{pb, rb}, i1_times_i2{pa_times_pb, ra_plus_rb}, i1_join_i2{pa_join_pb, ra_join_rb}; }; -TEST_F(PolarTest, PositiveInterval) +TEST_F(PolarTest, positive_interval_operations) { - EXPECT_NO_THROW(pa + pb); - EXPECT_NO_THROW(pa * pb); - EXPECT_NO_THROW(pa | pb); + EXPECT_EQ(pa + pb, pa_plus_pb); + EXPECT_EQ(pa * pb, pa_times_pb); + EXPECT_EQ(pa | pb, pa_join_pb); } -TEST_F(PolarTest, AngleInterval) +TEST_F(PolarTest, angle_interval_operations) { - EXPECT_NO_THROW(ra + rb); - EXPECT_NO_THROW(ra * rb); - EXPECT_NO_THROW(ra | rb); + EXPECT_EQ(ra, ra_other); + EXPECT_EQ(ra + rb, ra_plus_rb); + // EXPECT_EQ(ra * rb, ra_times_rb); + std::cout << ra.max() << " " << rb.max() << std::endl; + EXPECT_EQ(ra | rb, ra_join_rb) << (ra | rb).min() << " ; " << (ra | rb).delta(); } -TEST_F(PolarTest, Operations) +TEST_F(PolarTest, interval_operations) { - EXPECT_NO_THROW(i1 * i2); - EXPECT_NO_THROW(i1 | i2); + EXPECT_EQ(i1 * i2, i1_times_i2); + EXPECT_EQ(i1 | i2, i1_join_i2); EXPECT_EQ(polar::zero, polar::zero * polar::one); } @@ -77,7 +78,7 @@ TEST_F(PolarTest, multiply) EXPECT_EQ(minus, one * minus); } -TEST_F(PolarTest, Add) +TEST_F(PolarTest, add) { auto mthree = polar::Interval(-3.); auto sum = polar::Interval(2.) + mthree; @@ -85,9 +86,9 @@ TEST_F(PolarTest, Add) EXPECT_EQ(sum, ref) << sum.to_string() << " vs " << ref.to_string(); } -TEST_F(PolarTest, Norm) +TEST_F(PolarTest, norm) { - EXPECT_EQ(i1.norm(), 1.); + EXPECT_EQ(i1.norm(), 3.125); EXPECT_EQ(i2.norm(), 0.); }