From 21d6c5ab8490495b8566f1571b2fc2b477d13cf3 Mon Sep 17 00:00:00 2001
From: Alexander Condello <arcondello@gmail.com>
Date: Wed, 20 May 2026 15:56:58 -0700
Subject: [PATCH 1/4] Add minimal interval implementation

---
 .../include/dwave-optimization/interval.hpp   | 90 +++++++++++++++++++
 dwave/optimization/src/interval.cpp           | 19 ++++
 meson.build                                   |  1 +
 tests/cpp/meson.build                         |  1 +
 tests/cpp/test_interval.cpp                   | 69 ++++++++++++++
 5 files changed, 180 insertions(+)
 create mode 100644 dwave/optimization/include/dwave-optimization/interval.hpp
 create mode 100644 dwave/optimization/src/interval.cpp
 create mode 100644 tests/cpp/test_interval.cpp

diff --git a/dwave/optimization/include/dwave-optimization/interval.hpp b/dwave/optimization/include/dwave-optimization/interval.hpp
new file mode 100644
index 000000000..840b4e304
--- /dev/null
+++ b/dwave/optimization/include/dwave-optimization/interval.hpp
@@ -0,0 +1,90 @@
+// Copyright 2026 D-Wave
+//
+//    Licensed under the Apache License, Version 2.0 (the "License");
+//    you may not use this file except in compliance with the License.
+//    You may obtain a copy of the License at
+//
+//        http://www.apache.org/licenses/LICENSE-2.0
+//
+//    Unless required by applicable law or agreed to in writing, software
+//    distributed under the License is distributed on an "AS IS" BASIS,
+//    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+//    See the License for the specific language governing permissions and
+//    limitations under the License.
+
+#pragma once
+
+#include <algorithm>
+#include <cassert>
+#include <cmath>
+
+#include "dwave-optimization/common.hpp"  // for ssize_t
+#include "dwave-optimization/typing.hpp"
+
+namespace dwave::optimization {
+
+// dev note: true interval arithmetic requires "outward" rounding. For right
+// now I am ignoring this in the interest of getting an API solidified, but
+// we should support it ASAP. I have marked the places we need to make changes
+// with a "numeric issues" comment for easy future reference.
+
+template <DType T>
+class interval {
+ public:
+    interval() noexcept : infimum(1), supremum(0) {};
+
+    // todo: note
+    interval(T value) noexcept : interval(value, value) {}
+
+    interval(T infimum, T supremum) noexcept : infimum(infimum), supremum(supremum) {}
+
+    static interval unbounded() {
+        using limits = std::numeric_limits<T>;
+        if constexpr (limits::has_infinity) {
+            return interval(-limits::infinity(), limits::infinity());
+        } else {
+            return interval(limits::lowest(), limits::max());
+        }
+    }
+
+    /// An interval evalutes to `true` if it is not empty.
+    explicit operator bool() const noexcept { return infimum <= supremum; }
+
+    bool operator==(const interval& rhs) const {
+        return infimum == rhs.infimum and supremum == rhs.supremum;
+    }
+
+    interval& operator|=(DType auto scalar) { return *this |= interval<decltype(scalar)>(scalar); }
+    template <DType U>
+    interval& operator|=(const interval<U>& rhs) {
+        // dev note: numeric issues
+        if constexpr (std::floating_point<U> and std::integral<T>) {
+            infimum = std::min<T>(infimum, std::floor(rhs.infimum));
+            supremum = std::max<T>(supremum, std::ceil(rhs.supremum));
+        } else {
+            infimum = std::min<T>(infimum, rhs.infimum);
+            supremum = std::max<T>(supremum, rhs.supremum);
+        }
+        return *this;
+    }
+
+    bool contains(auto&& value) const {
+        // dev note: numeric issues
+        return infimum <= value and value <= supremum;
+    }
+
+    /// The size of the interval is the max(supremum - infimum, 0)
+    double size() const noexcept requires(std::floating_point<T>) {
+        if (supremum >= infimum) return static_cast<double>(supremum) - infimum;
+        return 0;
+    }
+    std::ptrdiff_t size() const noexcept requires(std::integral<T>) {
+        if (supremum >= infimum) return static_cast<std::ptrdiff_t>(supremum) - infimum;
+        return 0;
+    }
+
+    T infimum;
+    T supremum;
+};
+
+}  // namespace dwave::optimization
diff --git a/dwave/optimization/src/interval.cpp b/dwave/optimization/src/interval.cpp
new file mode 100644
index 000000000..c7e962fc4
--- /dev/null
+++ b/dwave/optimization/src/interval.cpp
@@ -0,0 +1,19 @@
+// Copyright 2026 D-Wave
+//
+//    Licensed under the Apache License, Version 2.0 (the "License");
+//    you may not use this file except in compliance with the License.
+//    You may obtain a copy of the License at
+//
+//        http://www.apache.org/licenses/LICENSE-2.0
+//
+//    Unless required by applicable law or agreed to in writing, software
+//    distributed under the License is distributed on an "AS IS" BASIS,
+//    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+//    See the License for the specific language governing permissions and
+//    limitations under the License.
+
+#include "dwave-optimization/interval.hpp"
+
+namespace dwave::optimization {
+
+}  // namespace dwave::optimization
diff --git a/meson.build b/meson.build
index 690d71652..eb2d252d2 100644
--- a/meson.build
+++ b/meson.build
@@ -52,6 +52,7 @@ dwave_optimization_src = [
     'dwave/optimization/src/array.cpp',
     'dwave/optimization/src/fraction.cpp',
     'dwave/optimization/src/graph.cpp',
+    'dwave/optimization/src/interval.cpp',
     'dwave/optimization/src/simplex.cpp',
 ]
 
diff --git a/tests/cpp/meson.build b/tests/cpp/meson.build
index ef23f7641..f9cb39251 100644
--- a/tests/cpp/meson.build
+++ b/tests/cpp/meson.build
@@ -35,6 +35,7 @@ tests_all = executable(
         'test_functional.cpp',
         'test_functional_.cpp',
         'test_graph.cpp',
+        'test_interval.cpp',
         'test_iterators.cpp',
         'test_simplex.cpp',
         'test_type_list.cpp',
diff --git a/tests/cpp/test_interval.cpp b/tests/cpp/test_interval.cpp
new file mode 100644
index 000000000..6f5644b94
--- /dev/null
+++ b/tests/cpp/test_interval.cpp
@@ -0,0 +1,69 @@
+// Copyright 2026 D-Wave
+//
+//    Licensed under the Apache License, Version 2.0 (the "License");
+//    you may not use this file except in compliance with the License.
+//    You may obtain a copy of the License at
+//
+//        http://www.apache.org/licenses/LICENSE-2.0
+//
+//    Unless required by applicable law or agreed to in writing, software
+//    distributed under the License is distributed on an "AS IS" BASIS,
+//    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+//    See the License for the specific language governing permissions and
+//    limitations under the License.
+
+#include <catch2/catch_template_test_macros.hpp>
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include "dwave-optimization/interval.hpp"
+
+namespace dwave::optimization {
+
+TEMPLATE_TEST_CASE(
+    "interval",
+    "",
+    float,
+    double,
+    bool,
+    std::int8_t,
+    std::int16_t,
+    std::int32_t,
+    std::int64_t
+) {
+    GIVEN("a default-constructed interval") {
+        auto in = interval<TestType>();
+
+        CHECK(not static_cast<bool>(in));
+        CHECK(in.size() == 0);
+    }
+
+    GIVEN("an interval containing only the value 0") {
+        auto in = interval<TestType>(0);
+
+        CHECK(static_cast<bool>(in));
+        CHECK(in.size() == 0);
+        CHECK(in.contains(0));
+        CHECK(not in.contains(1));
+        CHECK(not in.contains(std::numeric_limits<double>::epsilon()));
+
+        auto [inf, sup] = in;
+        STATIC_REQUIRE(std::same_as<decltype(inf), TestType>);
+        STATIC_REQUIRE(std::same_as<decltype(sup), TestType>);
+        CHECK(inf == 0);
+        CHECK(sup == 0);
+    }
+}
+
+TEST_CASE("interval") {
+    SECTION("union operations") {
+        auto in = interval<int>();
+        in |= .5;
+        CHECK(in == interval<int>(0, 1));
+        in |= -5;
+        CHECK(in == interval<int>(-5, 1));
+    }
+    // CHECK((interval<int>() | static_cast<float>(.5)) == interval<float>(0, 1));
+}
+
+}  // namespace dwave::optimization

From 317fae808a832e6bcac8b8a384f30f2ae66d94c6 Mon Sep 17 00:00:00 2001
From: Alexander Condello <arcondello@gmail.com>
Date: Thu, 21 May 2026 16:21:53 -0700
Subject: [PATCH 2/4] Expand interval implementation

---
 .../include/dwave-optimization/interval.hpp   | 180 +++++++++++++++---
 dwave/optimization/src/interval.cpp           |  11 ++
 ...feature-cpp-interval-232c1ed8d06f7fec.yaml |   3 +
 tests/cpp/test_interval.cpp                   |  66 ++++++-
 4 files changed, 227 insertions(+), 33 deletions(-)
 create mode 100644 releasenotes/notes/feature-cpp-interval-232c1ed8d06f7fec.yaml

diff --git a/dwave/optimization/include/dwave-optimization/interval.hpp b/dwave/optimization/include/dwave-optimization/interval.hpp
index 840b4e304..36c9365ca 100644
--- a/dwave/optimization/include/dwave-optimization/interval.hpp
+++ b/dwave/optimization/include/dwave-optimization/interval.hpp
@@ -17,6 +17,8 @@
 #include <algorithm>
 #include <cassert>
 #include <cmath>
+#include <iosfwd>
+#include <type_traits>
 
 #include "dwave-optimization/common.hpp"  // for ssize_t
 #include "dwave-optimization/typing.hpp"
@@ -24,67 +26,193 @@
 namespace dwave::optimization {
 
 // dev note: true interval arithmetic requires "outward" rounding. For right
-// now I am ignoring this in the interest of getting an API solidified, but
+// now I am ignoring this and focusing on getting the API solidified, but
 // we should support it ASAP. I have marked the places we need to make changes
 // with a "numeric issues" comment for easy future reference.
 
 template <DType T>
-class interval {
+requires(not std::is_const_v<T>) class interval {
  public:
-    interval() noexcept : infimum(1), supremum(0) {};
-
-    // todo: note
-    interval(T value) noexcept : interval(value, value) {}
+    // by default we construct an empty interval
+    interval() = default;
+
+    // rule of five stuff all works the way you'd expect
+    interval(const interval&) = default;
+    interval(interval&&) = default;
+    ~interval() = default;
+    interval& operator=(const interval&) = default;
+    interval& operator=(interval&&) = default;
+
+    // todo: note the difference between interval<int>(.5, .5) and interval<int>(interval<float>(.5,
+    // .5))
+    explicit interval(T value) noexcept : interval(value, value) {}
+    interval(T infimum, T supremum) noexcept : infimum_(infimum), supremum_(supremum) {
+        assert(not std::isnan(infimum_) and "infimum must not be nan");
+        assert(not std::isnan(supremum_) and "supremum must not be nan");
+    }
 
-    interval(T infimum, T supremum) noexcept : infimum(infimum), supremum(supremum) {}
+    template <DType U>
+    explicit(std::integral<T> != std::integral<U>) interval(const interval<U>& rhs) : interval() {
+        *this = rhs;
+    }
+    template <DType U>
+    explicit(std::integral<T> != std::integral<U>) interval(interval<U>&& rhs) : interval() {
+        *this = rhs;
+    }
 
-    static interval unbounded() {
-        using limits = std::numeric_limits<T>;
-        if constexpr (limits::has_infinity) {
-            return interval(-limits::infinity(), limits::infinity());
+    template <DType U>
+    interval& operator=(const interval<U>& rhs) {
+        // dev note: numeric issues
+        if constexpr (std::floating_point<U> and std::integral<T>) {
+            infimum_ = std::floor(rhs.infimum());
+            supremum_ = std::ceil(rhs.supremum());
         } else {
-            return interval(limits::lowest(), limits::max());
+            infimum_ = rhs.infimum();
+            supremum_ = rhs.supremum();
         }
+        return *this;
     }
 
     /// An interval evalutes to `true` if it is not empty.
-    explicit operator bool() const noexcept { return infimum <= supremum; }
+    explicit operator bool() const noexcept { return infimum_ <= supremum_; }
 
     bool operator==(const interval& rhs) const {
-        return infimum == rhs.infimum and supremum == rhs.supremum;
+        return infimum_ == rhs.infimum_ and supremum_ == rhs.supremum_;
     }
 
     interval& operator|=(DType auto scalar) { return *this |= interval<decltype(scalar)>(scalar); }
+
     template <DType U>
     interval& operator|=(const interval<U>& rhs) {
+        // if rhs is an empty interval, do nothing
+        if (not static_cast<bool>(rhs)) return *this;
+
+        // if lhs is an empty interval, then set it to rhs (which we know is not empty)
+        if (not static_cast<bool>(*this)) return *this = rhs;
+
+        // otherwise, adjust our interval to also contain the other.
         // dev note: numeric issues
-        if constexpr (std::floating_point<U> and std::integral<T>) {
-            infimum = std::min<T>(infimum, std::floor(rhs.infimum));
-            supremum = std::max<T>(supremum, std::ceil(rhs.supremum));
+        if constexpr (std::integral<T> and std::floating_point<U>) {
+            infimum_ = std::min<T>(infimum_, std::floor(rhs.infimum()));
+            supremum_ = std::max<T>(supremum_, std::ceil(rhs.supremum()));
         } else {
-            infimum = std::min<T>(infimum, rhs.infimum);
-            supremum = std::max<T>(supremum, rhs.supremum);
+            infimum_ = std::min<T>(infimum_, rhs.infimum());
+            supremum_ = std::max<T>(supremum_, rhs.supremum());
         }
         return *this;
     }
 
-    bool contains(auto&& value) const {
+    template <DType U>
+    friend auto operator|(interval lhs, const interval<U>& rhs) {
+        interval<decltype(T() + U())> out(std::move(lhs));
+        out |= rhs;
+        return out;
+    }
+
+    template <DType U>
+    interval& operator+=(const interval<U>& rhs) {
+        // if lhs is an empty interval than the result is also empty
+        if (not static_cast<bool>(*this)) return *this;
+
+        // if rhs is an empty interval, then set lhs to rhs (making lhs empty)
+        if (not static_cast<bool>(rhs)) return *this = rhs;
+
+        // neither are empty! So add everything
         // dev note: numeric issues
-        return infimum <= value and value <= supremum;
+        if constexpr (std::integral<T> and std::floating_point<U>) {
+            infimum_ += std::floor(rhs.infimum());
+            supremum_ += std::ceil(rhs.supremum());
+        } else {
+            infimum_ += rhs.infimum();
+            supremum_ += rhs.supremum();
+        }
+        return *this;
     }
 
-    /// The size of the interval is the max(supremum - infimum, 0)
+    template <DType U>
+    friend auto operator+(interval lhs, const interval<U>& rhs) {
+        interval<decltype(T() + U())> out(std::move(lhs));
+        out += rhs;
+        return out;
+    }
+
+    bool contains(DType auto&& value) const {
+        // dev note: numeric issues
+        return infimum_ <= value and value <= supremum_;
+    }
+
+    // Implement the get() to support the structured binding. It's not OK to
+    // overload std functions so this needs to be found via ADL.
+    // Unfortunately that prevents interval from satisfying the tuple-like
+    // concept. Why can't we have nice things C++?? See note below for
+    // more details.
+    template <std::size_t I>
+    friend const T& get(const interval& in) noexcept requires(I <= 1) {
+        if constexpr (I == 0) return in.infimum_;
+        if constexpr (I == 1) return in.supremum_;
+    }
+    template <std::size_t I>
+    friend T&& get(interval&& in) noexcept requires(I <= 1) {
+        if constexpr (I == 0) return std::move(in.infimum_);
+        if constexpr (I == 1) return std::move(in.supremum_);
+    }
+    template <std::size_t I>
+    friend const T&& get(const interval&& in) noexcept requires(I <= 1) {
+        if constexpr (I == 0) return std::move(in.infimum_);
+        if constexpr (I == 1) return std::move(in.supremum_);
+    }
+
+    const T& infimum() const noexcept { return infimum_; }
+
+    /// The size of the interval is the max(supremum_ - infimum_, 0)
     double size() const noexcept requires(std::floating_point<T>) {
-        if (supremum >= infimum) return static_cast<double>(supremum) - infimum;
+        if (supremum_ >= infimum_) return static_cast<double>(supremum_) - infimum_;
         return 0;
     }
     std::ptrdiff_t size() const noexcept requires(std::integral<T>) {
-        if (supremum >= infimum) return static_cast<std::ptrdiff_t>(supremum) - infimum;
+        if (supremum_ >= infimum_) return static_cast<std::ptrdiff_t>(supremum_) - infimum_;
         return 0;
     }
 
-    T infimum;
-    T supremum;
+    const T& supremum() const noexcept { return supremum_; }
+
+    static interval unbounded() {
+        using limits = std::numeric_limits<T>;
+        if constexpr (limits::has_infinity) {
+            return interval(-limits::infinity(), limits::infinity());
+        } else {
+            return interval(limits::lowest(), limits::max());
+        }
+    }
+
+ private:
+    T infimum_ = 1;
+    T supremum_ = 0;
 };
 
+std::ostream& operator<<(std::ostream& os, const interval<std::int64_t>& rhs);
+std::ostream& operator<<(std::ostream& os, const interval<double>& rhs);
+
 }  // namespace dwave::optimization
+
+// Do a lot of convoluted stuff to support structured binding.
+// Unfortunately this doesn't support the full tuple-like concept that would let us work
+// out of the box with a lot of nice methods.
+// See https://stackoverflow.com/questions/79660771/stdapply-to-a-custom-tuple-like-type
+// See https://lists.isocpp.org/std-discussion/2021/09/1438.php
+namespace std {
+
+template <dwave::optimization::DType T>
+struct tuple_size<dwave::optimization::interval<T>> : integral_constant<size_t, 2> {};
+
+template <dwave::optimization::DType T>
+struct tuple_element<0, dwave::optimization::interval<T>> {
+    using type = T;
+};
+
+template <dwave::optimization::DType T>
+struct tuple_element<1, dwave::optimization::interval<T>> {
+    using type = T;
+};
+
+}  // namespace std
diff --git a/dwave/optimization/src/interval.cpp b/dwave/optimization/src/interval.cpp
index c7e962fc4..cbfb04fa5 100644
--- a/dwave/optimization/src/interval.cpp
+++ b/dwave/optimization/src/interval.cpp
@@ -13,7 +13,18 @@
 //    limitations under the License.
 
 #include "dwave-optimization/interval.hpp"
+#include "dwave-optimization/common.hpp"  // for ssize_t
+#include "dwave-optimization/typing.hpp"
+
+#include <iostream>
 
 namespace dwave::optimization {
 
+std::ostream& operator<<(std::ostream& os, const interval<std::int64_t>& rhs) {
+    return os << "[" << rhs.infimum() << ", " << rhs.supremum() << "]";
+}
+std::ostream& operator<<(std::ostream& os, const interval<double>& rhs) {
+    return os << "[" << rhs.infimum() << ", " << rhs.supremum() << "]";
+}
+
 }  // namespace dwave::optimization
diff --git a/releasenotes/notes/feature-cpp-interval-232c1ed8d06f7fec.yaml b/releasenotes/notes/feature-cpp-interval-232c1ed8d06f7fec.yaml
new file mode 100644
index 000000000..d0369db2a
--- /dev/null
+++ b/releasenotes/notes/feature-cpp-interval-232c1ed8d06f7fec.yaml
@@ -0,0 +1,3 @@
+---
+features:
+  - Add C++ ``interval<T>()`` class for interval arithmetic.
diff --git a/tests/cpp/test_interval.cpp b/tests/cpp/test_interval.cpp
index 6f5644b94..95dcd9dd0 100644
--- a/tests/cpp/test_interval.cpp
+++ b/tests/cpp/test_interval.cpp
@@ -15,9 +15,12 @@
 #include <catch2/catch_template_test_macros.hpp>
 #include <catch2/catch_test_macros.hpp>
 #include <catch2/matchers/catch_matchers_all.hpp>
+#include <catch2/matchers/catch_matchers_range_equals.hpp>
 
 #include "dwave-optimization/interval.hpp"
 
+using Catch::Matchers::RangeEquals;
+
 namespace dwave::optimization {
 
 TEMPLATE_TEST_CASE(
@@ -46,24 +49,73 @@ TEMPLATE_TEST_CASE(
         CHECK(in.contains(0));
         CHECK(not in.contains(1));
         CHECK(not in.contains(std::numeric_limits<double>::epsilon()));
+    }
+
+    SECTION("get<...>(interval)") {
+        SECTION("const rvalue") {
+            const interval<TestType> in{0, 1};
+            const TestType&& inf = get<0>(std::move(in));
+            CHECK(inf == 0);
+        }
+    }
+
+    SECTION("printing") {
+        std::stringstream ss;
+        ss << interval<TestType>(0, 1);
+        CHECK(ss.str() == "[0, 1]");
+    }
+
+    SECTION("structured binding") {
+        SECTION("const reference") {
+            auto in = interval<TestType>(0, 1);
+            const auto& [inf, sup] = in;
+            STATIC_REQUIRE(std::same_as<decltype(inf), const TestType>);
+            STATIC_REQUIRE(std::same_as<decltype(sup), const TestType>);
+            CHECK(inf == 0);
+            CHECK(sup == 1);
+        }
 
-        auto [inf, sup] = in;
-        STATIC_REQUIRE(std::same_as<decltype(inf), TestType>);
-        STATIC_REQUIRE(std::same_as<decltype(sup), TestType>);
-        CHECK(inf == 0);
-        CHECK(sup == 0);
+        SECTION("rvalue") {
+            auto in = interval<TestType>(0, 1);
+            auto [inf, sup] = in;
+            STATIC_REQUIRE(std::same_as<decltype(inf), TestType>);
+            STATIC_REQUIRE(std::same_as<decltype(sup), TestType>);
+            CHECK(inf == 0);
+            CHECK(sup == 1);
+        }
+
+        SECTION("const rvalue") {
+            const auto in = interval<TestType>(0, 1);
+
+            auto&& [inf, sup] = in;
+            STATIC_REQUIRE(std::same_as<decltype(inf), const TestType>);
+            STATIC_REQUIRE(std::same_as<decltype(sup), const TestType>);
+            CHECK(inf == 0);
+            CHECK(sup == 1);
+        }
     }
 }
 
 TEST_CASE("interval") {
-    SECTION("union operations") {
+    SECTION("addition") {
+        CHECK((interval<int>(0, 3) + interval<float>(.5, .5)) == interval<float>(.5, 3.5));
+    }
+
+    SECTION("union") {
         auto in = interval<int>();
         in |= .5;
         CHECK(in == interval<int>(0, 1));
         in |= -5;
         CHECK(in == interval<int>(-5, 1));
+
+        CHECK((interval<int>(0, 1) | interval<double>(.5, 1.5)) == interval<double>(0, 1.5));
+        CHECK((interval<float>(0, 1) | interval<double>(.5, 1.5)) == interval<double>(0, 1.5));
+        CHECK((interval<float>(0, 1.5) | interval<bool>(0, 1)) == interval<float>(0, 1.5));
+    }
+
+    SECTION("Copy construction") {
+        CHECK(interval<int>(interval<float>(.2, .8)) == interval<int>(0, 1));
     }
-    // CHECK((interval<int>() | static_cast<float>(.5)) == interval<float>(0, 1));
 }
 
 }  // namespace dwave::optimization

From fc5de9afe2082ae017191ddcd7ae8dac7a208f11 Mon Sep 17 00:00:00 2001
From: Alexander Condello <arcondello@gmail.com>
Date: Tue, 26 May 2026 12:35:29 -0700
Subject: [PATCH 3/4] Add a proof-of-concept UnaryOp that uses the new interval
 class

---
 .../include/dwave-optimization/interval.hpp   | 96 ++++++++++++++-----
 dwave/optimization/src/functional_.hpp        | 21 ++++
 tests/cpp/test_functional_.cpp                | 36 ++++++-
 tests/cpp/test_interval.cpp                   |  3 +-
 4 files changed, 128 insertions(+), 28 deletions(-)

diff --git a/dwave/optimization/include/dwave-optimization/interval.hpp b/dwave/optimization/include/dwave-optimization/interval.hpp
index 36c9365ca..8c1eccd5c 100644
--- a/dwave/optimization/include/dwave-optimization/interval.hpp
+++ b/dwave/optimization/include/dwave-optimization/interval.hpp
@@ -43,8 +43,7 @@ requires(not std::is_const_v<T>) class interval {
     interval& operator=(const interval&) = default;
     interval& operator=(interval&&) = default;
 
-    // todo: note the difference between interval<int>(.5, .5) and interval<int>(interval<float>(.5,
-    // .5))
+    // todo: note the difference between interval<int>(.5, .5) and union operations
     explicit interval(T value) noexcept : interval(value, value) {}
     interval(T infimum, T supremum) noexcept : infimum_(infimum), supremum_(supremum) {
         assert(not std::isnan(infimum_) and "infimum must not be nan");
@@ -52,24 +51,17 @@ requires(not std::is_const_v<T>) class interval {
     }
 
     template <DType U>
-    explicit(std::integral<T> != std::integral<U>) interval(const interval<U>& rhs) : interval() {
-        *this = rhs;
-    }
+    explicit(std::integral<T> != std::integral<U>) interval(const interval<U>& rhs) :
+        infimum_(rhs.infimum()), supremum_(rhs.supremum()) {}
     template <DType U>
-    explicit(std::integral<T> != std::integral<U>) interval(interval<U>&& rhs) : interval() {
-        *this = rhs;
-    }
+    explicit(std::integral<T> != std::integral<U>) interval(interval<U>&& rhs) :
+        infimum_(rhs.infimum()), supremum_(rhs.supremum()) {}
 
     template <DType U>
     interval& operator=(const interval<U>& rhs) {
         // dev note: numeric issues
-        if constexpr (std::floating_point<U> and std::integral<T>) {
-            infimum_ = std::floor(rhs.infimum());
-            supremum_ = std::ceil(rhs.supremum());
-        } else {
-            infimum_ = rhs.infimum();
-            supremum_ = rhs.supremum();
-        }
+        infimum_ = rhs.infimum();
+        supremum_ = rhs.supremum();
         return *this;
     }
 
@@ -80,6 +72,44 @@ requires(not std::is_const_v<T>) class interval {
         return infimum_ == rhs.infimum_ and supremum_ == rhs.supremum_;
     }
 
+    interval operator-() const { return interval(-supremum_, -infimum_); }
+
+    template <DType U>
+    interval& operator&=(const interval<U>& rhs) {
+        // if lhs is an empty interval than the result is also empty
+        if (not static_cast<bool>(*this)) return *this;
+
+        // if rhs is an empty interval, then set lhs to rhs (making lhs empty)
+        if (not static_cast<bool>(rhs)) {
+            // if lhs is an empty interval, then set it to rhs (which we know is not empty)
+            if constexpr (std::integral<T> and std::floating_point<U>) {
+                infimum_ = std::floor(rhs.infimum());
+                supremum_ = std::ceil(rhs.supremum());
+            } else {
+                infimum_ = rhs.infimum();
+                supremum_ = rhs.supremum();
+            }
+        } else {
+            // otherwise, adjust our interval to be the union
+            if constexpr (std::integral<T> and std::floating_point<U>) {
+                infimum_ = std::max<T>(infimum_, std::floor(rhs.infimum()));
+                supremum_ = std::min<T>(supremum_, std::ceil(rhs.supremum()));
+            } else {
+                infimum_ = std::max<T>(infimum_, rhs.infimum());
+                supremum_ = std::min<T>(supremum_, rhs.supremum());
+            }
+        }
+
+        return *this;
+    }
+
+    template <DType U>
+    friend auto operator&(interval lhs, const interval<U>& rhs) {
+        interval<decltype(T() + U())> out(std::move(lhs));
+        out &= rhs;
+        return out;
+    }
+
     interval& operator|=(DType auto scalar) { return *this |= interval<decltype(scalar)>(scalar); }
 
     template <DType U>
@@ -87,18 +117,27 @@ requires(not std::is_const_v<T>) class interval {
         // if rhs is an empty interval, do nothing
         if (not static_cast<bool>(rhs)) return *this;
 
-        // if lhs is an empty interval, then set it to rhs (which we know is not empty)
-        if (not static_cast<bool>(*this)) return *this = rhs;
-
-        // otherwise, adjust our interval to also contain the other.
         // dev note: numeric issues
-        if constexpr (std::integral<T> and std::floating_point<U>) {
-            infimum_ = std::min<T>(infimum_, std::floor(rhs.infimum()));
-            supremum_ = std::max<T>(supremum_, std::ceil(rhs.supremum()));
+        if (not static_cast<bool>(*this)) {
+            // if lhs is an empty interval, then set it to rhs (which we know is not empty)
+            if constexpr (std::integral<T> and std::floating_point<U>) {
+                infimum_ = std::floor(rhs.infimum());
+                supremum_ = std::ceil(rhs.supremum());
+            } else {
+                infimum_ = rhs.infimum();
+                supremum_ = rhs.supremum();
+            }
         } else {
-            infimum_ = std::min<T>(infimum_, rhs.infimum());
-            supremum_ = std::max<T>(supremum_, rhs.supremum());
+            // otherwise, adjust our interval to also contain the other.
+            if constexpr (std::integral<T> and std::floating_point<U>) {
+                infimum_ = std::min<T>(infimum_, std::floor(rhs.infimum()));
+                supremum_ = std::max<T>(supremum_, std::ceil(rhs.supremum()));
+            } else {
+                infimum_ = std::min<T>(infimum_, rhs.infimum());
+                supremum_ = std::max<T>(supremum_, rhs.supremum());
+            }
         }
+
         return *this;
     }
 
@@ -164,6 +203,15 @@ requires(not std::is_const_v<T>) class interval {
 
     const T& infimum() const noexcept { return infimum_; }
 
+    static interval nonnegative() {
+        using limits = std::numeric_limits<T>;
+        if constexpr (limits::has_infinity) {
+            return interval(0, limits::infinity());
+        } else {
+            return interval(0, limits::max());
+        }
+    }
+
     /// The size of the interval is the max(supremum_ - infimum_, 0)
     double size() const noexcept requires(std::floating_point<T>) {
         if (supremum_ >= infimum_) return static_cast<double>(supremum_) - infimum_;
diff --git a/dwave/optimization/src/functional_.hpp b/dwave/optimization/src/functional_.hpp
index e1ad70821..57fb4aea7 100644
--- a/dwave/optimization/src/functional_.hpp
+++ b/dwave/optimization/src/functional_.hpp
@@ -24,8 +24,11 @@
 
 #pragma once
 
+#include <cmath>
+
 #include "dwave-optimization/array.hpp"
 #include "dwave-optimization/functional.hpp"
+#include "dwave-optimization/interval.hpp"
 #include "dwave-optimization/typing.hpp"
 
 namespace dwave::optimization::functional_ {
@@ -87,6 +90,24 @@ struct BinaryFunctionMixin {
     }
 };
 
+template <typename UnaryFunction>
+struct UnaryFunctionMixin {};
+
+template <DType T>
+struct Abs : UnaryFunctionMixin<Abs<T>> {
+    using result_type = T;
+
+    constexpr T operator()(DType auto&& x) const noexcept {
+        return std::abs(x);
+    }
+
+    template<DType U>
+    interval<T> operator()(const interval<U>& in) const noexcept {
+        // first get the absolute value, using the type of the input interval
+        return static_cast<interval<T>>((-in | in) & interval<U>::nonnegative());
+    }
+};
+
 /// Add two elements.
 template <DType T>
 struct Add : BinaryFunctionMixin<Add<double>> {
diff --git a/tests/cpp/test_functional_.cpp b/tests/cpp/test_functional_.cpp
index 773981d1d..cedad9f58 100644
--- a/tests/cpp/test_functional_.cpp
+++ b/tests/cpp/test_functional_.cpp
@@ -12,13 +12,43 @@
 //    See the License for the specific language governing permissions and
 //    limitations under the License.
 
-#include "catch2/catch_template_test_macros.hpp"
-#include "catch2/catch_test_macros.hpp"
-#include "catch2/generators/catch_generators.hpp"
+#include <catch2/catch_template_test_macros.hpp>
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/generators/catch_generators.hpp>
+
+#include "dwave-optimization/interval.hpp"
 #include "functional_.hpp"
 
 namespace dwave::optimization::functional_ {
 
+TEST_CASE("Abs") {
+    SECTION(".operator(...)") {
+        SECTION("Abs<double>") {
+            const Abs<double> op;
+
+            CHECK(op(1.5) == 1.5);
+            CHECK(op(-1) == 1);
+            CHECK(op(static_cast<const double>(1.5)) == 1.5);
+
+            CHECK(op(interval(-1, 1)) == interval<double>(0, 1));
+            CHECK(op(interval(-12.5, 1.)) == interval<double>(0, 12.5));
+        }
+        SECTION("Abs<int>") {
+            const Abs<int> op;
+
+            CHECK(op(1.999) == 1);  // going to int truncates
+            CHECK(op(-1.999) == 1);  // going to int truncates
+            CHECK(op(-1) == 1);
+            CHECK(op(static_cast<const double>(-1.9999999)) == 1);
+
+            CHECK(op(interval(-1, 1)) == interval<int>(0, 1));
+
+            // because abs itself truncates, we want our bounds to also be truncated
+            CHECK(op(interval(-12.5, 1.)) == interval<int>(0, 12));
+        }
+    }
+}
+
 TEST_CASE("Add") {
     static_assert(Add<double>::associative);
     static_assert(Add<double>::commutative);
diff --git a/tests/cpp/test_interval.cpp b/tests/cpp/test_interval.cpp
index 95dcd9dd0..823a618ac 100644
--- a/tests/cpp/test_interval.cpp
+++ b/tests/cpp/test_interval.cpp
@@ -114,7 +114,8 @@ TEST_CASE("interval") {
     }
 
     SECTION("Copy construction") {
-        CHECK(interval<int>(interval<float>(.2, .8)) == interval<int>(0, 1));
+        // truncates
+        CHECK(interval<int>(interval<float>(.2, .8)) == interval<int>(0, 0));
     }
 }
 

From 0bc57442212b473b5dcc86d2506f75c91b85b555 Mon Sep 17 00:00:00 2001
From: William Bernoudy <wbernoudy@dwavesys.com>
Date: Fri, 5 Jun 2026 14:03:34 -0700
Subject: [PATCH 4/4] Fix contains and size for intervals

and add Abs::inverse()
---
 .../include/dwave-optimization/interval.hpp       |  4 ++--
 dwave/optimization/src/functional_.hpp            | 15 +++++++++++----
 2 files changed, 13 insertions(+), 6 deletions(-)

diff --git a/dwave/optimization/include/dwave-optimization/interval.hpp b/dwave/optimization/include/dwave-optimization/interval.hpp
index 8c1eccd5c..4bcb4cc56 100644
--- a/dwave/optimization/include/dwave-optimization/interval.hpp
+++ b/dwave/optimization/include/dwave-optimization/interval.hpp
@@ -175,7 +175,7 @@ requires(not std::is_const_v<T>) class interval {
         return out;
     }
 
-    bool contains(DType auto&& value) const {
+    bool contains(auto&& value) const {
         // dev note: numeric issues
         return infimum_ <= value and value <= supremum_;
     }
@@ -218,7 +218,7 @@ requires(not std::is_const_v<T>) class interval {
         return 0;
     }
     std::ptrdiff_t size() const noexcept requires(std::integral<T>) {
-        if (supremum_ >= infimum_) return static_cast<std::ptrdiff_t>(supremum_) - infimum_;
+        if (supremum_ >= infimum_) return static_cast<std::ptrdiff_t>(supremum_) - infimum_ + 1;
         return 0;
     }
 
diff --git a/dwave/optimization/src/functional_.hpp b/dwave/optimization/src/functional_.hpp
index 57fb4aea7..fcf452ece 100644
--- a/dwave/optimization/src/functional_.hpp
+++ b/dwave/optimization/src/functional_.hpp
@@ -97,15 +97,22 @@ template <DType T>
 struct Abs : UnaryFunctionMixin<Abs<T>> {
     using result_type = T;
 
-    constexpr T operator()(DType auto&& x) const noexcept {
-        return std::abs(x);
-    }
+    constexpr T operator()(DType auto&& x) const noexcept { return std::abs(x); }
 
-    template<DType U>
+    template <DType U>
     interval<T> operator()(const interval<U>& in) const noexcept {
         // first get the absolute value, using the type of the input interval
         return static_cast<interval<T>>((-in | in) & interval<U>::nonnegative());
     }
+
+    template <DType U>
+    interval<T> inverse(const interval<U>& in) const noexcept {
+        assert(
+            (interval<U>::nonnegative() | in) == in and
+            "argument for inverse of abs should be nonnegative"
+        );
+        return in | -in;
+    }
 };
 
 /// Add two elements.