diff --git a/dwave/optimization/include/dwave-optimization/functional.hpp b/dwave/optimization/include/dwave-optimization/functional.hpp
index e050a7c11..9201d43e8 100644
--- a/dwave/optimization/include/dwave-optimization/functional.hpp
+++ b/dwave/optimization/include/dwave-optimization/functional.hpp
@@ -15,45 +15,185 @@
 #pragma once
 
 #include <algorithm>
+#include <cassert>
 #include <cmath>
 #include <concepts>
 #include <cstdlib>
+#include <limits>
+#include <utility>
+
+#include "dwave-optimization/interval.hpp"
+#include "dwave-optimization/typing.hpp"
 
 namespace dwave::optimization::functional {
 
-template <class T>
-struct abs {
-    static constexpr T operator()(const T& x) { return std::abs(x); }
+enum class Monotonicity { Decreasing = -1, None = 0, Increasing = 1 };
+
+template <typename UnaryOp>
+struct UnaryOpMixin {
+    template <DType T>
+    requires(UnaryOp::monotonic != Monotonicity::None)
+    static auto operator()(const interval<T>& domain) {
+        using return_type = interval<decltype(UnaryOp::operator()(T()))>;
+
+        // op(empty domain) -> empty domain
+        if (not static_cast<bool>(domain)) return return_type();
+
+        assert(
+            domain <= UnaryOp::template domain<T> and
+            "input domain must be a subset of the func's domain"
+        );
+
+        // We don't worry about outward rounding here because this overload is meant
+        // to reflect the behavior of the scalar overload, not necessarily to be
+        // mathematically correct.
+        // We *do* assume that UnaryOp (e.g., std::exp()) is monotonic, which
+        // is not always true, but I think it's an OK assumption for our purposes.
+        if constexpr (UnaryOp::monotonic == Monotonicity::Increasing) {
+            return return_type(
+                UnaryOp::operator()(domain.infimum), UnaryOp::operator()(domain.supremum)
+            );
+        } else if constexpr (UnaryOp::monotonic == Monotonicity::Decreasing) {
+            return return_type(
+                UnaryOp::operator()(domain.supremum), UnaryOp::operator()(domain.infimum)
+            );
+        } else {
+            assert(false and "unexpected monotonicity");
+            std::unreachable();
+        }
+    }
+
+    template <DType T>
+    static constexpr interval<T> domain = interval<T>::all();
 };
 
-template <class T>
-struct cos {
-    static auto operator()(const T& num) { return std::cos(num); }
+struct absolute : UnaryOpMixin<absolute> {
+    template <DType T>
+    static T operator()(const T& x) {
+        // Unlike NumPy/std, we define std::abs(INT_MIN) to equal INT_MAX under the reasoning
+        // that it's more important to us to preserve the sign than to preseve the correct value.
+        if constexpr (std::integral<T>) {
+            if (x == std::numeric_limits<T>::lowest()) return std::numeric_limits<T>::max();
+        }
+
+        // std::abs() is not defined for int8 or int16 so we static_cast to avoid widening.
+        return static_cast<T>(std::abs(x));
+    }
+    static bool operator()(const bool& x) { return x; }
+
+    template <DType T>
+    static interval<T> operator()(const interval<T>& domain) {
+        if (not static_cast<bool>(domain)) return {};  // op(empty domain) -> empty domain
+
+        assert(domain.infimum <= domain.supremum);  // implied by non-empty
+
+        // If the domain is non-negative, then absolute is identity
+        if (0 <= domain.infimum) return domain;
+
+        // If the domain is negative, then absolute is just the inverse
+        if (domain.supremum < 0) return -domain;
+
+        // Otherwise, the domain straddles 0
+
+        // Handle the -INT_MIN case. Again we treat abs(-INT_MIN) as INT_MAX under the reasoning
+        // that [INT_MIN, ...] is probably intended to mean unbounded.
+        if constexpr (std::integral<T>) {
+            if (domain.infimum == std::numeric_limits<T>::lowest()) {
+                return interval<T>(0, std::numeric_limits<T>::max());
+            }
+        }
+
+        return interval<T>(
+            0, -domain.infimum < domain.supremum ? domain.supremum : -domain.infimum
+        );
+    }
+    static interval<bool> operator()(const interval<bool>& domain) { return domain; }
+
+    static constexpr Monotonicity monotonic = Monotonicity::None;
 };
 
-template <class T>
-struct exp {
-    static constexpr auto operator()(const T& x) { return std::exp(x); }
+struct cos : UnaryOpMixin<cos> {
+    static auto operator()(const DType auto& x) { return std::cos(x); }
+
+    template <DType T>
+    static interval<decltype(std::cos(T()))> operator()(const interval<T>& domain) {
+        if (not static_cast<bool>(domain)) return {};  // op(empty domain) -> empty domain
+
+        // It is possible to be a lot more specific than this by checking whether
+        // our domain spans a full period or not, but I think this is of dubious
+        // benefit to the user so for now we just return [-1, +1]
+        return {-1, +1};
+    }
+
+    static constexpr Monotonicity monotonic = Monotonicity::None;
 };
 
-template <class T>
-struct expit {
-    static constexpr double operator()(const T& x) { return 1.0 / (1.0 + std::exp(-1. * x)); }
+struct exp : UnaryOpMixin<exp> {
+    static auto operator()(const DType auto& x) { return std::exp(x); }
+    using UnaryOpMixin::operator();
+
+    static constexpr Monotonicity monotonic = Monotonicity::Increasing;
 };
 
-template <class T>
-struct log {
-    static constexpr auto operator()(const T& x) { return std::log(x); }
+struct expit : UnaryOpMixin<expit> {
+    template <DType T>
+    static auto operator()(const T& x) {
+        return 1 / (1 + std::exp(-x));
+    }
+    using UnaryOpMixin::operator();
+
+    static constexpr Monotonicity monotonic = Monotonicity::Increasing;
 };
 
-template <class T>
-struct logical {
-    static constexpr bool operator()(const T& x) { return x; }
+struct log : UnaryOpMixin<log> {
+    template <DType T>
+    static auto operator()(const T& x) {
+        assert(domain<T>.contains(x) and "x must be non-negative");
+        return std::log(x);
+    }
+    using UnaryOpMixin::operator();
+
+    template <DType T>
+    static constexpr interval<T> domain = interval<T>::nonnegative();
+
+    static constexpr Monotonicity monotonic = Monotonicity::Increasing;
+};
+
+struct logical : UnaryOpMixin<logical> {
+    static bool operator()(const DType auto& x) { return x; }
+
+    static interval<bool> operator()(const interval<bool>& domain) { return domain; }
+    template <DType T>
+    static interval<bool> operator()(const interval<T>& domain) {
+        if (not static_cast<bool>(domain)) return {};  // op(empty domain) -> empty domain
+
+        if (domain.infimum == 0 and domain.supremum == 0) return interval(false, false);
+        if (domain.infimum <= 0 and domain.supremum >= 0) return interval(false, true);
+        return interval(true, true);
+    }
+
+    static constexpr Monotonicity monotonic = Monotonicity::None;
+};
+
+struct logical_not : UnaryOpMixin<logical_not> {
+    static bool operator()(const DType auto& x) { return not x; }
+
+    static interval<bool> operator()(const interval<bool>& domain) {
+        if (not static_cast<bool>(domain)) return {};  // op(empty domain) -> empty domain
+        return interval(not domain.supremum, not domain.infimum);
+    }
+    template <DType T>
+    static interval<bool> operator()(const interval<T>& domain) {
+        // Call the more specific interval<bool> overload
+        return operator()(logical{}(domain));
+    }
+
+    static constexpr Monotonicity monotonic = Monotonicity::None;
 };
 
 template <class T>
 struct logical_xor {
-    static constexpr bool operator()(const T& x, const T& y) {
+    static bool operator()(const T& x, const T& y) {
         return static_cast<bool>(x) != static_cast<bool>(y);
     }
 };
@@ -90,9 +230,28 @@ struct modulus {
     }
 };
 
-template <class T>
-struct rint {
-    static constexpr auto operator()(const T& x) { return std::rint(x); }
+struct negative : UnaryOpMixin<negative> {
+    template <class T>
+    requires(DType<T> and not std::same_as<T, bool>)  // not defined for bool
+    static auto operator()(const T& x) {
+        // We define -INT_MIN to equal INT_MAX under the reasoning that it's more
+        // important to us to preserve the sign than to preseve the correct value.
+        if constexpr (std::integral<T>) {
+            if (x == std::numeric_limits<T>::lowest()) return std::numeric_limits<T>::max();
+        }
+
+        return static_cast<T>(-x);  // so it doesn't widen e.g., int8_t->int
+    }
+    using UnaryOpMixin::operator();
+
+    static constexpr Monotonicity monotonic = Monotonicity::Decreasing;
+};
+
+struct rint : UnaryOpMixin<rint> {
+    static auto operator()(const DType auto& x) { return std::rint(x); }
+    using UnaryOpMixin::operator();
+
+    static constexpr Monotonicity monotonic = Monotonicity::Increasing;
 };
 
 template <class T>
@@ -103,24 +262,73 @@ struct safe_divides {
     }
 };
 
-template <class T>
-struct sin {
-    static auto operator()(const T& num) { return std::sin(num); }
+struct sin : UnaryOpMixin<sin> {
+    static auto operator()(const DType auto& x) { return std::sin(x); }
+
+    template <DType T>
+    static interval<decltype(std::sin(T()))> operator()(const interval<T>& domain) {
+        if (not static_cast<bool>(domain)) return {};  // op(empty domain) -> empty domain
+
+        // It is possible to be a lot more specific than this by checking whether
+        // our domain spans a full period or not, but I think this is of dubious
+        // benefit to the user so for now we just return [-1, +1]
+        return {-1, +1};
+    }
+
+    static constexpr Monotonicity monotonic = Monotonicity::None;
 };
 
-template <class T>
-struct square {
-    static constexpr T operator()(const T& x) { return x * x; }
+struct square : UnaryOpMixin<square> {
+    template <DType T>
+    static T operator()(const T& x) {
+        return x * x;
+    }
+    static bool operator()(const bool& x) { return x; }
+
+    template <DType T>
+    static interval<T> operator()(const interval<T>& domain) {
+        if (not static_cast<bool>(domain)) return {};  // op(empty domain) -> empty domain
+
+        assert(domain.infimum <= domain.supremum);  // implied by non-empty
+
+        square op{};
+        T inf_squared = op(domain.infimum);
+        T sup_squared = op(domain.supremum);
+
+        // Non-negative domain: square is increasing
+        if (0 <= domain.infimum) return interval<T>(inf_squared, sup_squared);
+
+        // Non-positive domain: square is decreasing
+        if (domain.supremum <= 0) return interval<T>(sup_squared, inf_squared);
+
+        // Otherwise the domain straddles 0: minimum is 0, maximum is the larger squared endpoint.
+
+        return interval<T>(0, inf_squared < sup_squared ? sup_squared : inf_squared);
+    }
+    static interval<bool> operator()(const interval<bool>& domain) { return domain; }
+
+    static constexpr Monotonicity monotonic = Monotonicity::None;
 };
 
-template <class T>
-struct square_root {
-    static constexpr auto operator()(const T& x) { return std::sqrt(x); }
+struct square_root : UnaryOpMixin<square_root> {
+    template <DType T>
+    static auto operator()(const T& x) {
+        assert(domain<T>.contains(x) and "x must be non-negative");
+        return std::sqrt(x);
+    }
+    using UnaryOpMixin::operator();
+
+    template <DType T>
+    static constexpr interval<T> domain = interval<T>::nonnegative();
+
+    static constexpr Monotonicity monotonic = Monotonicity::Increasing;
 };
 
-template <class T>
-struct tanh {
-    static auto operator()(const T& num) { return std::tanh(num); }
+struct tanh : UnaryOpMixin<tanh> {
+    static auto operator()(const DType auto& num) { return std::tanh(num); }
+    using UnaryOpMixin::operator();
+
+    static constexpr Monotonicity monotonic = Monotonicity::Increasing;
 };
 
 }  // namespace dwave::optimization::functional
diff --git a/dwave/optimization/include/dwave-optimization/nodes/unaryop.hpp b/dwave/optimization/include/dwave-optimization/nodes/unaryop.hpp
index 7e2795339..ea71d4af6 100644
--- a/dwave/optimization/include/dwave-optimization/nodes/unaryop.hpp
+++ b/dwave/optimization/include/dwave-optimization/nodes/unaryop.hpp
@@ -15,51 +15,82 @@
 #pragma once
 
 #include <cassert>
-#include <functional>
-#include <optional>
+#include <concepts>
+#include <memory>
+#include <ranges>
 #include <span>
+#include <stdexcept>
 #include <utility>
+#include <vector>
 
 #include "dwave-optimization/array.hpp"
 #include "dwave-optimization/functional.hpp"
 #include "dwave-optimization/graph.hpp"
+#include "dwave-optimization/interval.hpp"
 #include "dwave-optimization/state.hpp"
 
 namespace dwave::optimization {
 
+/// Apply a unary function to the predecessor node.
 template <class UnaryOp>
 class UnaryOpNode : public ArrayOutputMixin<ArrayNode> {
  public:
-    explicit UnaryOpNode(ArrayNode* node_ptr);
+    template <class... Args>
+    explicit UnaryOpNode(ArrayNode* node_ptr, Args&&... args) :
+        ArrayOutputMixin(node_ptr->shape()),
+        op(std::forward<Args>(args)...),
+        array_ptr_(node_ptr),
+        domain_(calculate_domain_(node_ptr)),
+        integral_(calculate_integral_(node_ptr)),
+        sizeinfo_(array_ptr_->sizeinfo()) {
+        add_predecessor_(node_ptr);
+    }
 
-    double const* buff(const State& state) const override;
-    std::span<const Update> diff(const State& state) const override;
+    /// @copydoc Array::buff()
+    double const* buff(const State& state) const override {
+        return data_ptr_<UnaryOpState_>(state)->buffer.data();
+    }
 
-    bool equal_to(const Node& rhs) const override;
-    bool equal_to(const UnaryOpNode& rhs) const;
+    /// @copydoc Node::commit()
+    void commit(State& state) const override {
+        auto* state_ptr = data_ptr_<UnaryOpState_>(state);
+        state_ptr->diff.clear();
+        state_ptr->previous_size = state_ptr->buffer.size();
+    }
 
-    /// @copydoc Array::integral()
-    bool integral() const override;
+    /// @copydoc Array::diff()
+    std::span<const Update> diff(const State& state) const override {
+        return data_ptr_<UnaryOpState_>(state)->diff;
+    }
 
-    /// @copydoc Array::min()
-    double min() const override;
+    /// @copydoc Node::equal_to()
+    bool equal_to(const Node& rhs) const override {
+        const UnaryOpNode* rhs_ptr = dynamic_cast<const UnaryOpNode*>(&rhs);
+        if (rhs_ptr == nullptr) return false;  // not same type so not equal
+        return this->equal_to(*rhs_ptr);       // use the equal_to(const UnaryOpNode&) overload
+    }
+    bool equal_to(const UnaryOpNode& rhs) const {
+        // If we're the same type, then we just need to check we have the same predecessor
+        return this->array_ptr_ == rhs.array_ptr_;
+    }
 
-    /// @copydoc Array::max()
-    double max() const override;
+    /// @copydoc Node::initialize_state()
+    void initialize_state(State& state) const override {
+        emplace_data_ptr_<UnaryOpState_>(
+            state, array_ptr_->view(state) | std::views::transform(op), array_ptr_->shape(state)
+        );
+    }
 
-    using ArrayOutputMixin::shape;
-    std::span<const ssize_t> shape(const State& state) const override;
-    using ArrayOutputMixin::size;
-    ssize_t size(const State& state) const override;
-    ssize_t size_diff(const State& state) const override;
-    SizeInfo sizeinfo() const override;
+    /// @copydoc Array::integral()
+    bool integral() const override { return integral_; }
 
-    void commit(State& state) const override;
-    void revert(State& state) const override;
-    void initialize_state(State& state) const override;
-    void propagate(State& state) const override;
+    /// @copydoc Array::max()
+    double max() const override { return domain_.supremum; }
 
-    // The predecessor of the operation, as an Array*.
+    /// @copydoc Array::min()
+    double min() const override { return domain_.infimum; }
+
+    /// The predecessor of the operation
     std::span<ArrayNode* const> operands() {
         assert(predecessors().size() == 1);
         return std::span<ArrayNode* const, 1>(&array_ptr_, 1);
@@ -69,8 +100,170 @@ class UnaryOpNode : public ArrayOutputMixin<ArrayNode> {
         return std::span<const ArrayNode* const, 1>(&array_ptr_, 1);
     }
 
+    /// @copydoc Node::propagate()
+    void propagate(State& state) const override {
+        const auto array_diff = array_ptr_->diff(state);
+
+        // If there are no updates the handle, return early.
+        if (array_diff.empty()) return;
+
+        auto* state_ptr = data_ptr_<UnaryOpState_>(state);
+        auto& buffer = state_ptr->buffer;
+        auto& diff = state_ptr->diff;
+
+        for (const auto& update : array_diff) {
+            assert(0 <= update.index);
+
+            if (update.removed()) {
+                assert(static_cast<size_t>(update.index) + 1 == buffer.size());
+
+                diff.emplace_back(Update::removal(update.index, buffer[update.index]));
+                buffer.pop_back();
+            } else if (update.placed()) {
+                assert(static_cast<size_t>(update.index) == buffer.size());
+
+                buffer.emplace_back(op(update.value));
+                diff.emplace_back(Update::placement(update.index, buffer.back()));
+            } else {
+                assert(static_cast<size_t>(update.index) < buffer.size());
+
+                double& old = buffer[update.index];
+                double value = op(update.value);
+
+                if (old == value) continue;  // no change to update
+
+                diff.emplace_back(update.index, old, value);
+                old = value;
+            }
+        }
+
+        if (ndim()) state_ptr->shape[0] = array_ptr_->shape(state)[0];
+
+        if (not diff.empty()) Node::propagate(state);
+    }
+
+    /// @copydoc Node::revert()
+    void revert(State& state) const override {
+        auto* state_ptr = data_ptr_<UnaryOpState_>(state);
+        std::vector<double>& buffer = state_ptr->buffer;
+        const ssize_t size = state_ptr->previous_size;
+        std::vector<Update>& diff = state_ptr->diff;
+        ssize_t* shape = state_ptr->shape.get();
+
+        const ssize_t propagated_size = buffer.size();
+
+        buffer.resize(size);
+        for (const auto& [index, old, _] : diff | std::views::reverse) {
+            assert(0 <= index);
+            if (size <= index) continue;
+            buffer[index] = old;
+        }
+        diff.clear();
+
+        // Adjust our shape back to what it should be, while avoiding reading our
+        // predecessor.
+        if (size != propagated_size) {
+            assert(0 < this->ndim());  // if our size changed this must be true
+            shape[0] = buffer.size();
+
+            // we could cache the divisor, but that feels like overkill in the context
+            // of a revert
+            for (ssize_t i = 1, ndim = this->ndim(); i < ndim; ++i) shape[0] /= shape[i];
+        }
+    }
+
+    /// @copydoc Array::shape()
+    std::span<const ssize_t> shape(const State& state) const override {
+        return std::span<const ssize_t>(data_ptr_<UnaryOpState_>(state)->shape.get(), ndim());
+    }
+    using ArrayOutputMixin::shape;
+
+    /// @copydoc Array::size()
+    ssize_t size(const State& state) const override {
+        return data_ptr_<UnaryOpState_>(state)->buffer.size();
+    }
+    using ArrayOutputMixin::size;
+
+    /// @copydoc Array::size_diff()
+    ssize_t size_diff(const State& state) const override {
+        const auto* state_ptr = data_ptr_<UnaryOpState_>(state);
+        return state_ptr->buffer.size() - state_ptr->previous_size;
+    }
+
+    /// @copydoc Array::sizeinfo()
+    SizeInfo sizeinfo() const override { return this->sizeinfo_; }
+
  private:
-    void replace_predecessor_(ssize_t index, Node* node_ptr) override;
+    // It would be more convenient to use ArrayStateData from _state.hpp
+    // but we don't currently have that in a public header.
+    // If that changes in the future, or we add default state support
+    // https://github.com/dwavesystems/dwave-optimization/issues/629 then
+    // we should use that.
+    struct UnaryOpState_ : NodeStateData {
+        UnaryOpState_() = delete;
+
+        template <std::ranges::range R>
+        UnaryOpState_(R&& values, std::span<const ssize_t> shape) :
+            buffer(std::ranges::to<std::vector<double>>(std::forward<R>(values))),
+            previous_size(buffer.size()),
+            shape(shape.size() ? std::make_unique<ssize_t[]>(shape.size()) : nullptr),
+            diff() {
+            for (ssize_t i = 0, ndim = shape.size(); i < ndim; ++i) {
+                this->shape[i] = shape[i];
+            }
+        }
+
+        std::vector<double> buffer;
+        ssize_t previous_size;
+
+        std::unique_ptr<ssize_t[]> shape;
+
+        std::vector<Update> diff;
+    };
+
+    // Calculate the min/max of the op applied to array_ptr
+    static interval<double> calculate_domain_(Array* array_ptr) {
+        interval<double> array_domain = interval<double>(array_ptr->min(), array_ptr->max());
+
+        if constexpr (std::same_as<UnaryOp, functional::square_root>) {
+            if (not(array_domain <= UnaryOp::template domain<double>)) {
+                throw std::invalid_argument("SquareRoot's predecessors must be non-negative");
+            }
+        }
+        if constexpr (std::same_as<UnaryOp, functional::log>) {
+            if (not(array_domain <= UnaryOp::template domain<double>)) {
+                throw std::invalid_argument("Log's predecessors must be non-negative");
+            }
+        }
+
+        // All the other ufuncs should have unbounded domains but as a sanity check...
+        assert(array_domain <= UnaryOp::template domain<double>);
+
+        // The range of our UnaryOp, i.e. the output min/max of our UnaryOpNode
+        return UnaryOp{}(array_domain);
+    }
+
+    // Determine whether the op applied to array_ptr will always result in an
+    // integral output.
+    static bool calculate_integral_(const Array* array_ptr) {
+        // rint() actually always returns a floating point. But it is an intergral
+        // floating point. This is a place where having proper dtypes would be
+        // very very nice.
+        if constexpr (std::same_as<UnaryOp, functional::rint>) return true;
+
+        // Otherwise, we ask our op what an integer input would result in.
+        if (array_ptr->integral()) {
+            return std::integral<decltype(UnaryOp{}(int()))>;
+        } else {
+            return std::integral<decltype(UnaryOp{}(double()))>;
+        }
+    }
+
+    void replace_predecessor_(ssize_t index, Node* node_ptr) override {
+        Node::replace_predecessor_(index, node_ptr);
+        array_ptr_ = dynamic_cast<ArrayNode*>(node_ptr);
+        assert(array_ptr_ != nullptr);
+    }
 
     UnaryOp op;
 
@@ -78,22 +271,49 @@ class UnaryOpNode : public ArrayOutputMixin<ArrayNode> {
     // predecessors(), but this is more performant
     ArrayNode* array_ptr_;
 
-    const ValuesInfo values_info_;
+    interval<double> domain_;  // the range of possible output values
+    bool integral_;            // whether the node will always output integral values
+
     const SizeInfo sizeinfo_;
 };
 
-using AbsoluteNode = UnaryOpNode<functional::abs<double>>;
-using CosNode = UnaryOpNode<functional::cos<double>>;
-using ExpitNode = UnaryOpNode<functional::expit<double>>;
-using ExpNode = UnaryOpNode<functional::exp<double>>;
-using LogNode = UnaryOpNode<functional::log<double>>;
-using LogicalNode = UnaryOpNode<functional::logical<double>>;
-using NegativeNode = UnaryOpNode<std::negate<double>>;
-using NotNode = UnaryOpNode<std::logical_not<double>>;
-using RintNode = UnaryOpNode<functional::rint<double>>;
-using SinNode = UnaryOpNode<functional::sin<double>>;
-using SquareNode = UnaryOpNode<functional::square<double>>;
-using SquareRootNode = UnaryOpNode<functional::square_root<double>>;
-using TanhNode = UnaryOpNode<functional::tanh<double>>;
+using AbsoluteNode = UnaryOpNode<functional::absolute>;
+extern template class UnaryOpNode<functional::absolute>;
+
+using CosNode = UnaryOpNode<functional::cos>;
+extern template class UnaryOpNode<functional::cos>;
+
+using ExpNode = UnaryOpNode<functional::exp>;
+extern template class UnaryOpNode<functional::exp>;
+
+using ExpitNode = UnaryOpNode<functional::expit>;
+extern template class UnaryOpNode<functional::expit>;
+
+using LogNode = UnaryOpNode<functional::log>;
+extern template class UnaryOpNode<functional::log>;
+
+using LogicalNode = UnaryOpNode<functional::logical>;
+extern template class UnaryOpNode<functional::logical>;
+
+using NegativeNode = UnaryOpNode<functional::negative>;
+extern template class UnaryOpNode<functional::negative>;
+
+using NotNode = UnaryOpNode<functional::logical_not>;
+extern template class UnaryOpNode<functional::logical_not>;
+
+using RintNode = UnaryOpNode<functional::rint>;
+extern template class UnaryOpNode<functional::rint>;
+
+using SinNode = UnaryOpNode<functional::sin>;
+extern template class UnaryOpNode<functional::sin>;
+
+using SquareNode = UnaryOpNode<functional::square>;
+extern template class UnaryOpNode<functional::square>;
+
+using SquareRootNode = UnaryOpNode<functional::square_root>;
+extern template class UnaryOpNode<functional::square_root>;
+
+using TanhNode = UnaryOpNode<functional::tanh>;
+extern template class UnaryOpNode<functional::tanh>;
 
 }  // namespace dwave::optimization
diff --git a/dwave/optimization/src/nodes/unaryop.cpp b/dwave/optimization/src/nodes/unaryop.cpp
index d871be910..71008a4f6 100644
--- a/dwave/optimization/src/nodes/unaryop.cpp
+++ b/dwave/optimization/src/nodes/unaryop.cpp
@@ -14,287 +14,20 @@
 
 #include "dwave-optimization/nodes/unaryop.hpp"
 
-#include <utility>
-
-#include "_state.hpp"
-
 namespace dwave::optimization {
 
-template <class UnaryOp>
-std::pair<double, double> calculate_values_minmax(const Array* array_ptr) {
-    // Do some checks to make sure the resulting domain/range will be valid
-    if constexpr (std::is_same<UnaryOp, functional::square_root<double>>::value) {
-        if (array_ptr->min() < 0) {
-            throw std::invalid_argument("SquareRoot's predecessors cannot take a negative value");
-        }
-    } else if constexpr (std::is_same<UnaryOp, functional::log<double>>::value) {
-        if (array_ptr->min() <= 0) {
-            throw std::invalid_argument("Log's predecessors cannot take a negative or zero value");
-        }
-    }
-
-    // If the output of the operation is boolean, then don't bother caching the result.
-    using result_type = typename std::invoke_result<UnaryOp, double&>::type;
-    if constexpr (std::same_as<result_type, bool>) {
-        return {false, true};
-    }
-
-    // Likewise for sin/cos/tanh the minmax is -1/+1. We could tighten it if the domain
-    // of our predecessor is smaller than 2pi, but let's keep it simple for now
-    if constexpr (
-        std::same_as<UnaryOp, functional::cos<double>> ||
-        std::same_as<UnaryOp, functional::sin<double>> ||
-        std::same_as<UnaryOp, functional::tanh<double>>
-    ) {
-        return {-1, +1};
-    }
-
-    // Otherwise the min and max depend on the predecessor
-
-    auto low = array_ptr->min();
-    auto high = array_ptr->max();
-    assert(low <= high);
-
-    if constexpr (std::same_as<UnaryOp, functional::abs<double>>) {
-        if (low >= 0 && high >= 0) {
-            return std::make_pair(low, high);
-        } else if (low >= 0) {
-            assert(false && "min > max");
-            std::unreachable();
-        } else if (high >= 0) {
-            return std::pair<double, double>(0.0, std::max<double>(-low, high));
-        } else {
-            return std::make_pair(-high, -low);
-        }
-    }
-    if constexpr (std::same_as<UnaryOp, functional::exp<double>>) {
-        return std::make_pair(std::exp(low), std::exp(high));
-    }
-    if constexpr (std::same_as<UnaryOp, functional::expit<double>>) {
-        double expit_low = 1.0 / (1.0 + std::exp(-low));
-        double expit_high = 1.0 / (1.0 + std::exp(-high));
-        return std::make_pair(expit_low, expit_high);
-    }
-    if constexpr (std::same_as<UnaryOp, functional::log<double>>) {
-        assert(low > 0);  // checked by constructor
-        return std::make_pair(std::log(low), std::log(high));
-    }
-    if constexpr (std::same_as<UnaryOp, functional::rint<double>>) {
-        return std::make_pair(std::rint(low), std::rint(high));
-    }
-    if constexpr (std::same_as<UnaryOp, functional::square<double>>) {
-        const auto highest = std::numeric_limits<double>::max();
-        return std::make_pair(
-            std::min({low * low, high * high, highest}),
-            std::min(
-                std::max({low * low, high * high}),
-                highest
-            )
-        );  // prevent inf
-    }
-    if constexpr (std::same_as<UnaryOp, functional::square_root<double>>) {
-        assert(low >= 0);  // checked by constructor
-        return std::make_pair(std::sqrt(low), std::sqrt(high));
-    }
-    if constexpr (std::same_as<UnaryOp, std::negate<double>>) {
-        return std::make_pair(-high, -low);
-    }
-
-    assert(false && "not implemeted yet");
-    std::unreachable();
-}
-
-template <class UnaryOp>
-bool calculate_integral(const Array*) {
-    using result_type = typename std::invoke_result<UnaryOp, double&>::type;
-    return std::is_integral<result_type>::value;
-}
-
-template <>
-bool calculate_integral<functional::abs<double>>(const Array* array_ptr) {
-    return array_ptr->integral();
-}
-
-template <>
-bool calculate_integral<functional::cos<double>>(const Array*) {
-    return false;
-}
-
-template <>
-bool calculate_integral<functional::exp<double>>(const Array*) {
-    return false;
-}
-
-template <>
-bool calculate_integral<functional::expit<double>>(const Array*) {
-    return false;
-}
-
-template <>
-bool calculate_integral<functional::log<double>>(const Array*) {
-    return false;
-}
-
-template <>
-bool calculate_integral<std::negate<double>>(const Array* array_ptr) {
-    return array_ptr->integral();
-}
-
-template <>
-bool calculate_integral<functional::rint<double>>(const Array*) {
-    return true;
-}
-
-template <>
-bool calculate_integral<functional::sin<double>>(const Array*) {
-    return false;
-}
-
-template <>
-bool calculate_integral<functional::square<double>>(const Array* array_ptr) {
-    return array_ptr->integral();
-}
-
-template <>
-bool calculate_integral<functional::tanh<double>>(const Array*) {
-    return false;
-}
-
-template <class UnaryOp>
-UnaryOpNode<UnaryOp>::UnaryOpNode(ArrayNode* node_ptr) :
-    ArrayOutputMixin(node_ptr->shape()),
-    array_ptr_(node_ptr),
-    values_info_(
-        calculate_values_minmax<UnaryOp>(array_ptr_),
-        calculate_integral<UnaryOp>(array_ptr_)
-    ),
-    sizeinfo_(array_ptr_->sizeinfo()) {
-    add_predecessor_(node_ptr);
-}
-
-template <class UnaryOp>
-void UnaryOpNode<UnaryOp>::commit(State& state) const {
-    data_ptr_<ArrayNodeStateData>(state)->commit();
-}
-
-template <class UnaryOp>
-double const* UnaryOpNode<UnaryOp>::buff(const State& state) const {
-    return data_ptr_<ArrayNodeStateData>(state)->buff();
-}
-
-template <class UnaryOp>
-std::span<const Update> UnaryOpNode<UnaryOp>::diff(const State& state) const {
-    return data_ptr_<ArrayNodeStateData>(state)->diff();
-}
-
-template <class UnaryOp>
-bool UnaryOpNode<UnaryOp>::equal_to(const Node& rhs) const {
-    const UnaryOpNode* rhs_ptr = dynamic_cast<const UnaryOpNode*>(&rhs);
-    if (rhs_ptr == nullptr) return false;  // not same type so not equal
-    return this->equal_to(*rhs_ptr);       // use the equal_to(const UnaryOpNode&) overload
-}
-
-template <class UnaryOp>
-bool UnaryOpNode<UnaryOp>::equal_to(const UnaryOpNode& rhs) const {
-    // If we're the same type, then we just need to check we have the same predecessor
-    return this->array_ptr_ == rhs.array_ptr_;
-}
-
-template <class UnaryOp>
-void UnaryOpNode<UnaryOp>::initialize_state(State& state) const {
-    std::vector<double> values;
-    values.reserve(array_ptr_->size(state));
-    for (const double val : array_ptr_->view(state)) {
-        values.emplace_back(op(val));
-    }
-
-    emplace_data_ptr_<ArrayNodeStateData>(state, std::move(values));
-}
-
-template <class UnaryOp>
-bool UnaryOpNode<UnaryOp>::integral() const {
-    return values_info_.integral;
-}
-
-template <class UnaryOp>
-double UnaryOpNode<UnaryOp>::min() const {
-    return this->values_info_.min;
-}
-
-template <class UnaryOp>
-double UnaryOpNode<UnaryOp>::max() const {
-    return this->values_info_.max;
-}
-
-template <class UnaryOp>
-void UnaryOpNode<UnaryOp>::propagate(State& state) const {
-    const auto diff = array_ptr_->diff(state);
-    // If there are no updates the handle, return early.
-    if (diff.empty()) return;
-
-    auto node_data = data_ptr_<ArrayNodeStateData>(state);
-
-    for (const auto& update : diff) {
-        const auto& [idx, _, value] = update;
-
-        if (update.placed()) {
-            assert(idx == static_cast<ssize_t>(node_data->size()));
-            node_data->emplace_back(op(value));
-        } else if (update.removed()) {
-            assert(idx == static_cast<ssize_t>(node_data->size()) - 1);
-            node_data->pop_back();
-        } else {
-            node_data->set(idx, op(value));
-        }
-    }
-
-    if (node_data->diff().size()) Node::propagate(state);
-}
-
-template <class UnaryOp>
-void UnaryOpNode<UnaryOp>::replace_predecessor_(ssize_t previous_index, Node* node_ptr) {
-    Node::replace_predecessor_(previous_index, node_ptr);
-    array_ptr_ = dynamic_cast<ArrayNode*>(node_ptr);
-    assert(array_ptr_ != nullptr);
-}
-
-template <class UnaryOp>
-void UnaryOpNode<UnaryOp>::revert(State& state) const {
-    data_ptr_<ArrayNodeStateData>(state)->revert();
-}
-
-template <class UnaryOp>
-std::span<const ssize_t> UnaryOpNode<UnaryOp>::shape(const State& state) const {
-    return array_ptr_->shape(state);
-}
-
-template <class UnaryOp>
-ssize_t UnaryOpNode<UnaryOp>::size(const State& state) const {
-    return data_ptr_<ArrayNodeStateData>(state)->size();
-}
-
-template <class UnaryOp>
-ssize_t UnaryOpNode<UnaryOp>::size_diff(const State& state) const {
-    return data_ptr_<ArrayNodeStateData>(state)->size_diff();
-}
-
-template <class UnaryOp>
-SizeInfo UnaryOpNode<UnaryOp>::sizeinfo() const {
-    return this->sizeinfo_;
-}
-
-template class UnaryOpNode<functional::abs<double>>;
-template class UnaryOpNode<functional::cos<double>>;
-template class UnaryOpNode<functional::exp<double>>;
-template class UnaryOpNode<functional::expit<double>>;
-template class UnaryOpNode<functional::log<double>>;
-template class UnaryOpNode<functional::logical<double>>;
-template class UnaryOpNode<functional::rint<double>>;
-template class UnaryOpNode<functional::sin<double>>;
-template class UnaryOpNode<functional::square<double>>;
-template class UnaryOpNode<functional::square_root<double>>;
-template class UnaryOpNode<functional::tanh<double>>;
-template class UnaryOpNode<std::negate<double>>;
-template class UnaryOpNode<std::logical_not<double>>;
+template class UnaryOpNode<functional::absolute>;
+template class UnaryOpNode<functional::cos>;
+template class UnaryOpNode<functional::exp>;
+template class UnaryOpNode<functional::expit>;
+template class UnaryOpNode<functional::log>;
+template class UnaryOpNode<functional::logical>;
+template class UnaryOpNode<functional::logical_not>;
+template class UnaryOpNode<functional::negative>;
+template class UnaryOpNode<functional::rint>;
+template class UnaryOpNode<functional::sin>;
+template class UnaryOpNode<functional::square>;
+template class UnaryOpNode<functional::square_root>;
+template class UnaryOpNode<functional::tanh>;
 
 }  // namespace dwave::optimization
diff --git a/releasenotes/notes/unaryop-rework-23877cf721e8f8a6.yaml b/releasenotes/notes/unaryop-rework-23877cf721e8f8a6.yaml
new file mode 100644
index 000000000..f4e8b7fdd
--- /dev/null
+++ b/releasenotes/notes/unaryop-rework-23877cf721e8f8a6.yaml
@@ -0,0 +1,9 @@
+---
+features:
+  - Make ``UnaryOpNode`` header-only.
+  - |
+    Rework C++ unary ops in ``dwave::optimization::functional`` to provide additional
+    information about the functions.
+upgrade:
+  - Rename C++ ``dwave::optimization::functional::abs()`` function to ``dwave::optimization::functional::absolute()``.
+  - Rename C++ ``dwave::optimization::functional::negate()`` function to ``dwave::optimization::functional::negative()``.
diff --git a/tests/cpp/nodes/test_unaryop.cpp b/tests/cpp/nodes/test_unaryop.cpp
index 8a960d279..373378a2a 100644
--- a/tests/cpp/nodes/test_unaryop.cpp
+++ b/tests/cpp/nodes/test_unaryop.cpp
@@ -33,26 +33,26 @@ namespace dwave::optimization {
 TEMPLATE_TEST_CASE(
     "UnaryOpNode",
     "",
-    functional::abs<double>,
-    functional::cos<double>,
-    functional::exp<double>,
-    functional::expit<double>,
-    functional::logical<double>,
-    functional::rint<double>,
-    functional::sin<double>,
-    functional::square<double>,
-    functional::tanh<double>,
-    std::negate<double>,
-    std::logical_not<double>
+    functional::absolute,
+    functional::cos,
+    functional::exp,
+    functional::expit,
+    functional::logical,
+    functional::logical_not,
+    functional::negative,
+    functional::rint,
+    functional::sin,
+    functional::square,
+    functional::tanh
 ) {
     auto graph = Graph();
 
     auto func = TestType();
 
     GIVEN("A constant scalar input") {
-        auto a_ptr = graph.emplace_node<ConstantNode>(-5);
+        auto* a_ptr = graph.emplace_node<ConstantNode>(-5);
 
-        auto p_ptr = graph.emplace_node<UnaryOpNode<TestType>>(a_ptr);
+        auto* p_ptr = graph.emplace_node<UnaryOpNode<TestType>>(a_ptr);
 
         THEN("The shape is also a scalar") {
             CHECK(p_ptr->ndim() == 0);
@@ -77,9 +77,9 @@ TEMPLATE_TEST_CASE(
     }
 
     GIVEN("A dynamic array input") {
-        auto a_ptr = graph.emplace_node<ListNode>(5, 0, 5);
+        auto* a_ptr = graph.emplace_node<ListNode>(5, 0, 5);
 
-        auto p_ptr = graph.emplace_node<UnaryOpNode<TestType>>(a_ptr);
+        auto* p_ptr = graph.emplace_node<UnaryOpNode<TestType>>(a_ptr);
 
         graph.emplace_node<ArrayValidationNode>(p_ptr);
 
@@ -139,13 +139,15 @@ TEMPLATE_TEST_CASE(
     }
 
     GIVEN("A 0-d integer decision input") {
-        auto a_ptr = graph.emplace_node<IntegerNode>(
+        auto* a_ptr = graph.emplace_node<IntegerNode>(
             std::span<const ssize_t>{},
             -100,
             100
         );  // Scalar output
 
-        auto p_ptr = graph.emplace_node<UnaryOpNode<TestType>>(a_ptr);
+        auto* p_ptr = graph.emplace_node<UnaryOpNode<TestType>>(a_ptr);
+
+        graph.emplace_node<ArrayValidationNode>(p_ptr);
 
         THEN("The integer is the operand") {
             CHECK(p_ptr->operands().size() == p_ptr->predecessors().size());
@@ -157,48 +159,42 @@ TEMPLATE_TEST_CASE(
             a_ptr->initialize_state(state, {-5});
             graph.initialize_state(state);
 
-            THEN("The output has the value and shape we expect") {
-                CHECK(p_ptr->size(state) == 1);
-                CHECK(p_ptr->shape(state).size() == 0);
-                CHECK(p_ptr->view(state)[0] == func(-5));
-            }
+            CHECK(p_ptr->size(state) == 1);
+            CHECK(p_ptr->shape(state).size() == 0);
+            CHECK_THAT(p_ptr->view(state), RangeEquals({func(-5)}));
 
             AND_WHEN("We change the integer's state and propagate") {
                 a_ptr->set_value(state, 0, 17);
-                a_ptr->propagate(state);
-                p_ptr->propagate(state);
+                graph.propagate(state);
+
+                CHECK_THAT(p_ptr->view(state), RangeEquals({func(17)}));
 
-                THEN("The output is what we expect") { CHECK(p_ptr->view(state)[0] == func(17)); }
 
                 AND_WHEN("We commit") {
-                    a_ptr->commit(state);
-                    p_ptr->commit(state);
+                    graph.commit(state);
 
-                    THEN("The value hasn't changed") { CHECK(p_ptr->view(state)[0] == func(17)); }
-                    THEN("The diffs are cleared") { CHECK(p_ptr->diff(state).size() == 0); }
+                    CHECK_THAT(p_ptr->view(state), RangeEquals({func(17)}));
+                    CHECK(p_ptr->diff(state).empty());
                 }
 
                 AND_WHEN("We revert") {
-                    a_ptr->revert(state);
-                    p_ptr->revert(state);
+                    graph.revert(state);
 
-                    THEN("The value reverts to the previous") {
-                        CHECK(p_ptr->view(state)[0] == func(-5));
-                    }
-                    THEN("The diffs are cleared") { CHECK(p_ptr->diff(state).size() == 0); }
+                    CHECK_THAT(p_ptr->view(state), RangeEquals({func(-5)}));
+                    CHECK(p_ptr->diff(state).empty());
                 }
             }
         }
     }
 
     GIVEN("A 3-d integer decision input") {
-        auto a_ptr = graph.emplace_node<IntegerNode>(
+        auto* a_ptr = graph.emplace_node<IntegerNode>(
             std::span<const ssize_t>({2, 3, 2}),
             -100,
             100
         );  // Scalar output
 
-        auto p_ptr = graph.emplace_node<UnaryOpNode<TestType>>(a_ptr);
+        auto* p_ptr = graph.emplace_node<UnaryOpNode<TestType>>(a_ptr);
 
         THEN("The integers node is the operand") {
             CHECK(p_ptr->operands().size() == p_ptr->predecessors().size());
@@ -281,8 +277,8 @@ TEST_CASE("UnaryOpNode - AbsoluteNode") {
     auto graph = Graph();
 
     GIVEN("An integer variable with domain [-3, 2]") {
-        auto i_ptr = graph.emplace_node<IntegerNode>(std::vector<ssize_t>{}, -3, 2);
-        auto abs_ptr = graph.emplace_node<AbsoluteNode>(i_ptr);
+        auto* i_ptr = graph.emplace_node<IntegerNode>(std::vector<ssize_t>{}, -3, 2);
+        auto* abs_ptr = graph.emplace_node<AbsoluteNode>(i_ptr);
 
         THEN("It has the min/max/integrality we expect") {
             CHECK(abs_ptr->min() == 0);
@@ -292,8 +288,8 @@ TEST_CASE("UnaryOpNode - AbsoluteNode") {
     }
 
     GIVEN("An integer variable with domain [-2, 4]") {
-        auto i_ptr = graph.emplace_node<IntegerNode>(std::vector<ssize_t>{}, -2, 4);
-        auto abs_ptr = graph.emplace_node<AbsoluteNode>(i_ptr);
+        auto* i_ptr = graph.emplace_node<IntegerNode>(std::vector<ssize_t>{}, -2, 4);
+        auto* abs_ptr = graph.emplace_node<AbsoluteNode>(i_ptr);
 
         THEN("It has the min/max we expect") {
             CHECK(abs_ptr->min() == 0);
@@ -308,14 +304,14 @@ TEST_CASE("UnaryOpNode - CosNode, SinNode, and TanhNode") {
     auto graph = Graph();
 
     GIVEN("x with min/max of -100/+100, y = cos(x), z = sin(x), a = tanh(x)") {
-        auto x = graph.emplace_node<ConstantNode>(std::vector<double>{-100, 1.5, +100});
-        auto y = graph.emplace_node<CosNode>(x);
-        auto z = graph.emplace_node<SinNode>(x);
-        auto a = graph.emplace_node<TanhNode>(x);
-
-        CHECK(!y->integral());
-        CHECK(!z->integral());
-        CHECK(!a->integral());
+        auto* x = graph.emplace_node<ConstantNode>(std::vector<double>{-100, 1.5, +100});
+        auto* y = graph.emplace_node<CosNode>(x);
+        auto* z = graph.emplace_node<SinNode>(x);
+        auto* a = graph.emplace_node<TanhNode>(x);
+
+        CHECK(not y->integral());
+        CHECK(not z->integral());
+        CHECK(not a->integral());
         CHECK(y->min() == -1);
         CHECK(y->max() == +1);
         CHECK(z->min() == -1);
@@ -329,8 +325,8 @@ TEST_CASE("UnaryOpNode - ExpitNode") {
     auto graph = Graph();
     GIVEN("An arbitrary number") {
         double c = 3.0;
-        auto c_ptr = graph.emplace_node<ConstantNode>(c);
-        auto expit_ptr = graph.emplace_node<ExpitNode>(c_ptr);
+        auto* c_ptr = graph.emplace_node<ConstantNode>(c);
+        auto* expit_ptr = graph.emplace_node<ExpitNode>(c_ptr);
         auto state = graph.initialize_state();
         CHECK(expit_ptr->min() == 1.0 / (1.0 + std::exp(-c)));
         CHECK(expit_ptr->max() == 1.0 / (1.0 + std::exp(-c)));
@@ -338,16 +334,16 @@ TEST_CASE("UnaryOpNode - ExpitNode") {
 
     GIVEN("A negative number") {
         double c = -4.0;
-        auto c_ptr = graph.emplace_node<ConstantNode>(c);
-        auto expit_ptr = graph.emplace_node<ExpitNode>(c_ptr);
+        auto* c_ptr = graph.emplace_node<ConstantNode>(c);
+        auto* expit_ptr = graph.emplace_node<ExpitNode>(c_ptr);
         auto state = graph.initialize_state();
         CHECK(expit_ptr->min() == 1.0 / (1.0 + std::exp(-c)));
         CHECK(expit_ptr->max() == 1.0 / (1.0 + std::exp(-c)));
     }
 
     GIVEN("A constant 1d array of doubles") {
-        auto c_ptr = graph.emplace_node<ConstantNode>(std::vector{-6.0, -0.3, 0.0, 1.2, 5.6});
-        auto expit_ptr = graph.emplace_node<ExpitNode>(c_ptr);
+        auto* c_ptr = graph.emplace_node<ConstantNode>(std::vector{-6.0, -0.3, 0.0, 1.2, 5.6});
+        auto* expit_ptr = graph.emplace_node<ExpitNode>(c_ptr);
 
         THEN("The min/max are expected") {
             THEN("expit(x) is not integral") { CHECK_FALSE(expit_ptr->integral()); }
@@ -357,8 +353,8 @@ TEST_CASE("UnaryOpNode - ExpitNode") {
     }
 
     GIVEN("An integer with max domain") {
-        auto i_ptr = graph.emplace_node<IntegerNode>(std::vector<ssize_t>{});
-        auto expit_ptr = graph.emplace_node<ExpitNode>(i_ptr);
+        auto* i_ptr = graph.emplace_node<IntegerNode>(std::vector<ssize_t>{});
+        auto* expit_ptr = graph.emplace_node<ExpitNode>(i_ptr);
         graph.emplace_node<ArrayValidationNode>(expit_ptr);
 
         THEN("The min/max are expected") {
@@ -372,8 +368,8 @@ TEST_CASE("UnaryOpNode - ExpNode") {
     auto graph = Graph();
     GIVEN("An arbitrary number") {
         double c = 3.0;
-        auto c_ptr = graph.emplace_node<ConstantNode>(c);
-        auto exp_ptr = graph.emplace_node<ExpNode>(c_ptr);
+        auto* c_ptr = graph.emplace_node<ConstantNode>(c);
+        auto* exp_ptr = graph.emplace_node<ExpNode>(c_ptr);
         auto state = graph.initialize_state();
         CHECK(exp_ptr->min() == std::exp(c));
         CHECK(exp_ptr->max() == std::exp(c));
@@ -381,16 +377,16 @@ TEST_CASE("UnaryOpNode - ExpNode") {
 
     GIVEN("A negative number") {
         double c = -4.0;
-        auto c_ptr = graph.emplace_node<ConstantNode>(c);
-        auto exp_ptr = graph.emplace_node<ExpNode>(c_ptr);
+        auto* c_ptr = graph.emplace_node<ConstantNode>(c);
+        auto* exp_ptr = graph.emplace_node<ExpNode>(c_ptr);
         auto state = graph.initialize_state();
         CHECK(exp_ptr->min() == std::exp(c));
         CHECK(exp_ptr->max() == std::exp(c));
     }
 
     GIVEN("A constant 1d array of doubles") {
-        auto c_ptr = graph.emplace_node<ConstantNode>(std::vector{-6.0, -0.3, 0.0, 1.2, 5.6});
-        auto exp_ptr = graph.emplace_node<ExpNode>(c_ptr);
+        auto* c_ptr = graph.emplace_node<ConstantNode>(std::vector{-6.0, -0.3, 0.0, 1.2, 5.6});
+        auto* exp_ptr = graph.emplace_node<ExpNode>(c_ptr);
 
         THEN("The min/max are expected") {
             THEN("exp(x) is not integral") { CHECK_FALSE(exp_ptr->integral()); }
@@ -400,8 +396,8 @@ TEST_CASE("UnaryOpNode - ExpNode") {
     }
 
     GIVEN("An integer with max domain") {
-        auto i_ptr = graph.emplace_node<IntegerNode>(std::vector<ssize_t>{});
-        auto exp_ptr = graph.emplace_node<ExpNode>(i_ptr);
+        auto* i_ptr = graph.emplace_node<IntegerNode>(std::vector<ssize_t>{});
+        auto* exp_ptr = graph.emplace_node<ExpNode>(i_ptr);
         graph.emplace_node<ArrayValidationNode>(exp_ptr);
 
         THEN("The min is expected") {
@@ -414,8 +410,8 @@ TEST_CASE("UnaryOpNode - ExpNode") {
 TEST_CASE("UnaryOpNode - LogNode") {
     auto graph = Graph();
     GIVEN("A constant 1d array of doubles") {
-        auto c_ptr = graph.emplace_node<ConstantNode>(std::vector{6.0, 0.3, 1.0, 1.2, 5.6});
-        auto log_ptr = graph.emplace_node<LogNode>(c_ptr);
+        auto* c_ptr = graph.emplace_node<ConstantNode>(std::vector{6.0, 0.3, 1.0, 1.2, 5.6});
+        auto* log_ptr = graph.emplace_node<LogNode>(c_ptr);
 
         THEN("The min/max are expected") {
             THEN("log(x) is not integral") { CHECK_FALSE(log_ptr->integral()); }
@@ -425,15 +421,15 @@ TEST_CASE("UnaryOpNode - LogNode") {
     }
     GIVEN("An arbitrary number") {
         double c = 10.0;
-        auto c_ptr = graph.emplace_node<ConstantNode>(c);
-        auto log_ptr = graph.emplace_node<LogNode>(c_ptr);
+        auto* c_ptr = graph.emplace_node<ConstantNode>(c);
+        auto* log_ptr = graph.emplace_node<LogNode>(c_ptr);
         auto state = graph.initialize_state();
         CHECK(log_ptr->min() == std::log(c));
         CHECK(log_ptr->max() == std::log(c));
     }
     GIVEN("A negative number") {
         double c = -10.0;
-        auto c_ptr = graph.emplace_node<ConstantNode>(c);
+        auto* c_ptr = graph.emplace_node<ConstantNode>(c);
         CHECK_THROWS(graph.emplace_node<LogNode>(c_ptr));
     }
 }
@@ -441,8 +437,8 @@ TEST_CASE("UnaryOpNode - LogNode") {
 TEST_CASE("UnaryOpNode - LogicalNode") {
     auto graph = Graph();
     GIVEN("A constant of mixed doubles and a negation of it") {
-        auto c_ptr = graph.emplace_node<ConstantNode>(std::vector{-2., -1., 0., 1., 2., -.5, .5});
-        auto logical_ptr = graph.emplace_node<LogicalNode>(c_ptr);
+        auto* c_ptr = graph.emplace_node<ConstantNode>(std::vector{-2., -1., 0., 1., 2., -.5, .5});
+        auto* logical_ptr = graph.emplace_node<LogicalNode>(c_ptr);
 
         THEN("NotNode is logical") {
             CHECK(logical_ptr->integral());
@@ -463,8 +459,8 @@ TEST_CASE("UnaryOpNode - LogicalNode") {
 TEST_CASE("UnaryOpNode - NegativeNode") {
     auto graph = Graph();
     GIVEN("An integer array and an asymmetric domain") {
-        auto i_ptr = graph.emplace_node<IntegerNode>(std::vector<ssize_t>{5}, -3, 8);
-        auto ni_ptr = graph.emplace_node<NegativeNode>(i_ptr);
+        auto* i_ptr = graph.emplace_node<IntegerNode>(std::vector<ssize_t>{5}, -3, 8);
+        auto* ni_ptr = graph.emplace_node<NegativeNode>(i_ptr);
 
         THEN("Negative has the min/max we expect") {
             CHECK(i_ptr->min() == -3);
@@ -480,8 +476,8 @@ TEST_CASE("UnaryOpNode - NegativeNode") {
 TEST_CASE("UnaryOpNode - NotNode") {
     auto graph = Graph();
     GIVEN("A constant of mixed doubles and a negation of it") {
-        auto c_ptr = graph.emplace_node<ConstantNode>(std::vector{-2., -1., 0., 1., 2., -.5, .5});
-        auto nc_ptr = graph.emplace_node<NotNode>(c_ptr);
+        auto* c_ptr = graph.emplace_node<ConstantNode>(std::vector{-2., -1., 0., 1., 2., -.5, .5});
+        auto* nc_ptr = graph.emplace_node<NotNode>(c_ptr);
 
         THEN("NotNode is logical") {
             CHECK(nc_ptr->integral());
@@ -503,8 +499,8 @@ TEST_CASE("UnaryOpNode - RintNode") {
     auto graph = Graph();
     GIVEN("An arbitrary number") {
         double c = 10.3;
-        auto c_ptr = graph.emplace_node<ConstantNode>(c);
-        auto rint_ptr = graph.emplace_node<RintNode>(c_ptr);
+        auto* c_ptr = graph.emplace_node<ConstantNode>(c);
+        auto* rint_ptr = graph.emplace_node<RintNode>(c_ptr);
         auto state = graph.initialize_state();
         CHECK(rint_ptr->min() == std::rint(c));
         CHECK(rint_ptr->max() == std::rint(c));
@@ -512,16 +508,16 @@ TEST_CASE("UnaryOpNode - RintNode") {
 
     GIVEN("A negative number") {
         double c = -10.5;
-        auto c_ptr = graph.emplace_node<ConstantNode>(c);
-        auto rint_ptr = graph.emplace_node<RintNode>(c_ptr);
+        auto* c_ptr = graph.emplace_node<ConstantNode>(c);
+        auto* rint_ptr = graph.emplace_node<RintNode>(c_ptr);
         auto state = graph.initialize_state();
         CHECK(rint_ptr->min() == std::rint(c));
         CHECK(rint_ptr->max() == std::rint(c));
     }
 
     GIVEN("A constant 1d array of doubles") {
-        auto c_ptr = graph.emplace_node<ConstantNode>(std::vector{-3.8, 0.3, 1.2, 5.6});
-        auto rint_ptr = graph.emplace_node<RintNode>(c_ptr);
+        auto* c_ptr = graph.emplace_node<ConstantNode>(std::vector{-3.8, 0.3, 1.2, 5.6});
+        auto* rint_ptr = graph.emplace_node<RintNode>(c_ptr);
 
         THEN("The min/max are expected") {
             CHECK(rint_ptr->integral());
@@ -532,18 +528,18 @@ TEST_CASE("UnaryOpNode - RintNode") {
 
     WHEN("We access a constant 1d array using integers from a RintNode") {
         double c0 = 0.3;
-        auto c0_ptr = graph.emplace_node<ConstantNode>(c0);
-        auto rint0_ptr = graph.emplace_node<RintNode>(c0_ptr);
+        auto* c0_ptr = graph.emplace_node<ConstantNode>(c0);
+        auto* rint0_ptr = graph.emplace_node<RintNode>(c0_ptr);
 
         double c3 = 2.8;
-        auto c3_ptr = graph.emplace_node<ConstantNode>(c3);
-        auto rint3_ptr = graph.emplace_node<RintNode>(c3_ptr);
+        auto* c3_ptr = graph.emplace_node<ConstantNode>(c3);
+        auto* rint3_ptr = graph.emplace_node<RintNode>(c3_ptr);
         ;
 
-        auto arr_ptr = graph.emplace_node<ConstantNode>(std::vector{0, 10, 20, 30});
+        auto* arr_ptr = graph.emplace_node<ConstantNode>(std::vector{0, 10, 20, 30});
 
-        auto a0_ptr = graph.emplace_node<AdvancedIndexingNode>(arr_ptr, rint0_ptr);
-        auto a3_ptr = graph.emplace_node<AdvancedIndexingNode>(arr_ptr, rint3_ptr);
+        auto* a0_ptr = graph.emplace_node<AdvancedIndexingNode>(arr_ptr, rint0_ptr);
+        auto* a3_ptr = graph.emplace_node<AdvancedIndexingNode>(arr_ptr, rint3_ptr);
 
         auto state = graph.initialize_state();
 
@@ -557,8 +553,8 @@ TEST_CASE("UnaryOpNode - RintNode") {
 TEST_CASE("UnaryOpNode - SquareNode") {
     auto graph = Graph();
     GIVEN("An integer with max domain") {
-        auto i_ptr = graph.emplace_node<IntegerNode>(std::vector<ssize_t>{});
-        auto square_ptr = graph.emplace_node<SquareNode>(i_ptr);
+        auto* i_ptr = graph.emplace_node<IntegerNode>(std::vector<ssize_t>{});
+        auto* square_ptr = graph.emplace_node<SquareNode>(i_ptr);
 
         THEN("The min/max are expected") {
             CHECK(square_ptr->min() == 0);
@@ -574,8 +570,8 @@ TEST_CASE("UnaryOpNode - SquareNode") {
 TEST_CASE("UnaryOpNode - SquareRootNode") {
     auto graph = Graph();
     GIVEN("An integer with max domain") {
-        auto i_ptr = graph.emplace_node<IntegerNode>(std::vector<ssize_t>{});
-        auto square_root_ptr = graph.emplace_node<SquareRootNode>(i_ptr);
+        auto* i_ptr = graph.emplace_node<IntegerNode>(std::vector<ssize_t>{});
+        auto* square_root_ptr = graph.emplace_node<SquareRootNode>(i_ptr);
         graph.emplace_node<ArrayValidationNode>(square_root_ptr);
 
         THEN("The min/max are expected") {
@@ -587,15 +583,15 @@ TEST_CASE("UnaryOpNode - SquareRootNode") {
     }
     GIVEN("An arbitrary number") {
         double c = 10.0;
-        auto c_ptr = graph.emplace_node<ConstantNode>(c);
-        auto square_root_ptr = graph.emplace_node<SquareRootNode>(c_ptr);
+        auto* c_ptr = graph.emplace_node<ConstantNode>(c);
+        auto* square_root_ptr = graph.emplace_node<SquareRootNode>(c_ptr);
         auto state = graph.initialize_state();
         CHECK(square_root_ptr->min() == std::sqrt(c));
         CHECK(square_root_ptr->max() == std::sqrt(c));
     }
     GIVEN("A negative number") {
         double c = -10.0;
-        auto c_ptr = graph.emplace_node<ConstantNode>(c);
+        auto* c_ptr = graph.emplace_node<ConstantNode>(c);
         CHECK_THROWS(graph.emplace_node<SquareRootNode>(c_ptr));
     }
 }
diff --git a/tests/cpp/test_functional.cpp b/tests/cpp/test_functional.cpp
index 931de059c..899a98ce9 100644
--- a/tests/cpp/test_functional.cpp
+++ b/tests/cpp/test_functional.cpp
@@ -12,24 +12,412 @@
 //    See the License for the specific language governing permissions and
 //    limitations under the License.
 
+#include <cmath>
+#include <concepts>
+#include <limits>
+
 #include <catch2/catch_template_test_macros.hpp>
 #include <catch2/catch_test_macros.hpp>
 
 #include "dwave-optimization/functional.hpp"
+#include "dwave-optimization/interval.hpp"
+#include "dwave-optimization/typing.hpp"
 
 namespace dwave::optimization::functional {
 
-TEMPLATE_TEST_CASE("modulus", "", double, int) {
+TEMPLATE_LIST_TEST_CASE("absolute", "", DTypes) {
+    constexpr absolute op{};
+
+    SECTION("absolute(scalar)") {
+        CHECK(op(TestType(0)) == 0);
+        CHECK(op(TestType(1)) == 1);
+
+        if constexpr (std::same_as<bool, TestType>) {
+            CHECK(op(true) == 1);  // abs(bool) is identity
+        } else if constexpr (std::integral<TestType>) {
+            CHECK(op(TestType(-1)) == 1);
+            CHECK(op(TestType(-10)) == 10);
+            CHECK(op(TestType(3)) == 3);
+            // We define abs(lowest) == max (see functional.hpp)
+            CHECK(
+                op(std::numeric_limits<TestType>::lowest()) == std::numeric_limits<TestType>::max()
+            );
+        } else {  // floating
+            CHECK(op(TestType(-1.5)) == 1.5);
+            CHECK(op(TestType(1.5)) == 1.5);
+        }
+    }
+
+    SECTION("absolute(interval)") {
+        CHECK(not op(interval<TestType>()));  // op(empty) -> empty
+
+        CHECK(op(interval<TestType>(0, 0)) == interval<TestType>(0, 0));
+        CHECK(op(interval<TestType>(0, 1)) == interval<TestType>(0, 1));
+
+        if constexpr (not std::same_as<TestType, bool>) {
+            CHECK(op(interval<TestType>(0, 5)) == interval<TestType>(0, 5));
+            CHECK(op(interval<TestType>(-7, -4)) == interval<TestType>(4, 7));
+            CHECK(op(interval<TestType>(-3, 1)) == interval<TestType>(0, 3));
+            CHECK(op(interval<TestType>(-1, 3)) == interval<TestType>(0, 3));
+            CHECK(op(interval<TestType>(-5, 5)) == interval<TestType>(0, 5));
+        }
+        if constexpr (std::floating_point<TestType>) {
+            CHECK(op(interval<TestType>(0.5, 5.2)) == interval<TestType>(0.5, 5.2));
+            CHECK(op(interval<TestType>(-5.2, -0.5)) == interval<TestType>(0.5, 5.2));
+        }
+    }
+}
+
+TEMPLATE_LIST_TEST_CASE("cos", "", DTypes) {
+    constexpr cos op{};
+
+    SECTION("cos(scalar)") {
+        CHECK(op(TestType(0)) == 1);  // cos(0) == 1 exactly
+        if constexpr (not std::same_as<bool, TestType>) {
+            CHECK(op(TestType(1)) == std::cos(TestType(1)));
+            CHECK(op(TestType(3)) == std::cos(TestType(3)));
+        }
+    }
+
+    SECTION("cos(interval)") {
+        CHECK(not op(interval<TestType>()));  // op(empty) -> empty
+        CHECK(op(interval<TestType>(0, 0)) == interval<double>(-1, +1));
+    }
+}
+
+TEMPLATE_LIST_TEST_CASE("exp", "", DTypes) {
+    constexpr exp op{};
+
+    SECTION("exp(scalar)") {
+        CHECK(op(TestType(0)) == 1);  // exp(0) == 1 exactly
+        if constexpr (not std::same_as<bool, TestType>) {
+            CHECK(op(TestType(1)) == std::exp(TestType(1)));
+            CHECK(op(TestType(-2)) == std::exp(TestType(-2)));
+        }
+    }
+
+    SECTION("exp(interval)") {
+        CHECK(not op(interval<TestType>()));  // op(empty) -> empty
+        CHECK(op(interval<TestType>(0, 1)) == interval(op(TestType(0)), op(TestType(1))));
+        if constexpr (not std::same_as<bool, TestType>) {
+            CHECK(op(interval<TestType>(-2, 3)) == interval(op(TestType(-2)), op(TestType(3))));
+        }
+    }
+
+    SECTION("exp domain is unrestricted") {
+        CHECK(exp::domain<TestType> == interval<TestType>::all());
+    }
+}
+
+TEMPLATE_LIST_TEST_CASE("expit", "", DTypes) {
+    constexpr expit op{};
+
+    SECTION("expit(scalar)") {
+        CHECK(op(TestType(0)) == 0.5);  // 1 / (1 + 1)
+        if constexpr (std::floating_point<TestType>) {
+            // no NaN at the extremes
+            CHECK(op(TestType(-1000)) == 0);
+            CHECK(op(TestType(1000)) == 1);
+        }
+    }
+
+    SECTION("expit(interval)") {
+        CHECK(not op(interval<TestType>()));  // op(empty) -> empty
+        CHECK(op(interval<TestType>(0, 1)) == interval(op(TestType(0)), op(TestType(1))));
+        if constexpr (not std::same_as<bool, TestType>) {
+            CHECK(op(interval<TestType>(-2, 3)) == interval(op(TestType(-2)), op(TestType(3))));
+        }
+    }
+}
+
+TEMPLATE_LIST_TEST_CASE("log", "", DTypes) {
+    constexpr log op{};
+
+    SECTION("log(scalar)") {
+        CHECK(op(TestType(1)) == 0);  // log(1) == 0 exactly
+        if constexpr (std::same_as<bool, TestType>) {
+        } else if constexpr (std::integral<TestType>) {
+            CHECK(op(TestType(2)) == std::log(TestType(2)));
+            CHECK(op(TestType(10)) == std::log(TestType(10)));
+        } else {  // floating
+            CHECK(op(TestType(2.5)) == std::log(TestType(2.5)));
+            CHECK(op(TestType(0.5)) == std::log(TestType(0.5)));
+        }
+    }
+
+    SECTION("log(interval)") {
+        CHECK(not op(interval<TestType>()));  // op(empty) -> empty
+        if constexpr (std::same_as<bool, TestType>) {
+            CHECK(op(interval<TestType>(1, 1)) == interval(op(TestType(1)), op(TestType(1))));
+        } else {
+            CHECK(op(interval<TestType>(1, 4)) == interval(op(TestType(1)), op(TestType(4))));
+            CHECK(op(interval<TestType>(2, 10)) == interval(op(TestType(2)), op(TestType(10))));
+        }
+    }
+
+    SECTION("log domain is non-negative") {
+        CHECK(log::domain<TestType> == interval<TestType>::nonnegative());
+    }
+}
+
+TEMPLATE_LIST_TEST_CASE("logical", "", DTypes) {
+    constexpr logical op{};
+
+    SECTION("logical(<scalar>)") {
+        CHECK(op(TestType(0)) == 0);
+
+        if constexpr (std::same_as<bool, TestType>) {
+            CHECK(op(true) == 1);
+        } else if constexpr (std::integral<TestType>) {
+            CHECK(op(TestType(-1)) == 1);
+            CHECK(op(TestType(1)) == 1);
+            CHECK(op(TestType(3)) == 1);
+        } else {  // floating
+            CHECK(op(TestType(-.000001)) == 1);
+            CHECK(op(TestType(.000001)) == 1);
+        }
+    }
+
+    SECTION("logical(<interval>)") {
+        CHECK(not op(interval<TestType>()));  // op(null) -> null
+
+        CHECK(op(interval<TestType>(0, 0)) == interval(false, false));
+        CHECK(op(interval<TestType>(1, 1)) == interval(true, true));
+        CHECK(op(interval<TestType>(0, 1)) == interval(false, true));
+
+        if constexpr (std::same_as<bool, TestType>) {
+            // already covered
+        } else if constexpr (std::integral<TestType>) {
+            CHECK(op(interval<TestType>(0, 5)) == interval(false, true));
+            CHECK(op(interval<TestType>(1, 5)) == interval(true, true));
+
+            CHECK(op(interval<TestType>(-3, 5)) == interval(false, true));
+
+            CHECK(op(interval<TestType>(-3, 0)) == interval(false, true));
+            CHECK(op(interval<TestType>(-3, -1)) == interval(true, true));
+        } else {  // floating
+            CHECK(op(interval<TestType>(0, .00001)) == interval(false, true));
+            CHECK(op(interval<TestType>(.000001, 5.5)) == interval(true, true));
+
+            CHECK(op(interval<TestType>(-3.4, 13.2)) == interval(false, true));
+
+            CHECK(op(interval<TestType>(-.00000001, 0)) == interval(false, true));
+            CHECK(op(interval<TestType>(-3.3, -.01)) == interval(true, true));
+        }
+    }
+}
+
+TEMPLATE_LIST_TEST_CASE("logical_not", "", DTypes) {
+    constexpr logical_not op{};
+
+    SECTION("logical_not(<scalar>)") {
+        CHECK(op(TestType(0)) == 1);
+
+        if constexpr (std::same_as<bool, TestType>) {
+            CHECK(op(true) == 0);
+        } else if constexpr (std::integral<TestType>) {
+            CHECK(op(TestType(-1)) == 0);
+            CHECK(op(TestType(1)) == 0);
+            CHECK(op(TestType(3)) == 0);
+        } else {  // floating
+            CHECK(op(TestType(-.000001)) == 0);
+            CHECK(op(TestType(.000001)) == 0);
+        }
+    }
+
+    SECTION("logical_not(<interval>)") {
+        CHECK(not op(interval<TestType>()));  // op(null) -> null
+
+        CHECK(op(interval<TestType>(0, 0)) == interval(true, true));
+        CHECK(op(interval<TestType>(1, 1)) == interval(false, false));
+        CHECK(op(interval<TestType>(0, 1)) == interval(false, true));
+
+        if constexpr (std::same_as<bool, TestType>) {
+            // already covered
+        } else if constexpr (std::integral<TestType>) {
+            CHECK(op(interval<TestType>(0, 5)) == interval(false, true));
+            CHECK(op(interval<TestType>(1, 5)) == interval(false, false));
+
+            CHECK(op(interval<TestType>(-3, 5)) == interval(false, true));
+
+            CHECK(op(interval<TestType>(-3, 0)) == interval(false, true));
+            CHECK(op(interval<TestType>(-3, -1)) == interval(false, false));
+        } else {  // floating
+            CHECK(op(interval<TestType>(0, .00001)) == interval(false, true));
+            CHECK(op(interval<TestType>(.000001, 5.5)) == interval(false, false));
+
+            CHECK(op(interval<TestType>(-3.4, 13.2)) == interval(false, true));
+
+            CHECK(op(interval<TestType>(-.00000001, 0)) == interval(false, true));
+            CHECK(op(interval<TestType>(-3.3, -.01)) == interval(false, false));
+        }
+    }
+}
+
+TEMPLATE_LIST_TEST_CASE("modulus", "", DTypes) {
     constexpr modulus<TestType> op;
 
-    // test for consistency with NumPy
     CHECK(op(1, 0) == 0);
     CHECK(op(0, 1) == 0);
-    CHECK(op(-1, 0) == 0);
-    CHECK(op(0, -1) == 0);
-    CHECK(op(-1, -10) == -1);
-    CHECK(op(-1, 10) == 9);
-    CHECK(op(1, -10) == -9);
+
+    if constexpr (not std::same_as<TestType, bool>) {
+        CHECK(op(-1, 0) == 0);
+        CHECK(op(0, -1) == 0);
+
+        CHECK(op(-1, -10) == -1);
+        CHECK(op(-1, 10) == 9);
+        CHECK(op(1, -10) == -9);
+    }
+}
+
+TEMPLATE_LIST_TEST_CASE("negative", "", DTypes) {
+    constexpr negative op{};
+    if constexpr (not std::same_as<TestType, bool>) {
+        SECTION("negative(scalar)") {
+            CHECK(op(TestType(0)) == 0);
+
+            if constexpr (std::integral<TestType>) {
+                CHECK(op(TestType(3)) == -3);
+                CHECK(op(TestType(-3)) == 3);
+            } else {  // floating
+                CHECK(op(TestType(1.5)) == -1.5);
+                CHECK(op(TestType(-1.5)) == 1.5);
+            }
+        }
+
+        SECTION("negative(interval)") {
+            CHECK(not op(interval<TestType>()));  // op(empty) -> empty
+
+            CHECK(op(interval<TestType>(0, 1)) == interval(op(TestType(1)), op(TestType(0))));
+            CHECK(op(interval<TestType>(-2, 3)) == interval(op(TestType(3)), op(TestType(-2))));
+            CHECK(op(interval<TestType>(-5, -1)) == interval(op(TestType(-1)), op(TestType(-5))));
+        }
+    }
+}
+
+TEMPLATE_LIST_TEST_CASE("rint", "", DTypes) {
+    constexpr rint op{};
+
+    SECTION("rint(scalar)") {
+        CHECK(op(TestType(0)) == 0);
+        if constexpr (std::same_as<bool, TestType>) {
+            CHECK(op(true) == 1);
+        } else if constexpr (std::integral<TestType>) {
+            CHECK(op(TestType(3)) == 3);
+            CHECK(op(TestType(-4)) == -4);
+        } else {  // floating: rounds half to even
+            CHECK(op(TestType(2.5)) == 2);
+            CHECK(op(TestType(3.5)) == 4);
+            CHECK(op(TestType(-2.5)) == -2);
+            CHECK(op(TestType(2.4)) == std::rint(TestType(2.4)));
+        }
+    }
+
+    SECTION("rint(interval)") {
+        CHECK(not op(interval<TestType>()));  // op(empty) -> empty
+        CHECK(op(interval<TestType>(0, 1)) == interval(op(TestType(0)), op(TestType(1))));
+        if constexpr (not std::same_as<bool, TestType>) {
+            CHECK(op(interval<TestType>(-3, 4)) == interval(op(TestType(-3)), op(TestType(4))));
+        }
+    }
+}
+
+TEMPLATE_LIST_TEST_CASE("sin", "", DTypes) {
+    constexpr sin op{};
+
+    SECTION("sin(scalar)") {
+        CHECK(op(TestType(0)) == 0);  // sin(0) == 0 exactly
+        if constexpr (not std::same_as<bool, TestType>) {
+            CHECK(op(TestType(1)) == std::sin(TestType(1)));
+            CHECK(op(TestType(2)) == std::sin(TestType(2)));
+        }
+    }
+
+    SECTION("sin(interval)") {
+        CHECK(not op(interval<TestType>()));  // op(empty) -> empty
+        CHECK(op(interval<TestType>(0, 0)) == interval<double>(-1, +1));
+    }
+}
+
+TEMPLATE_LIST_TEST_CASE("square", "", DTypes) {
+    constexpr square op{};
+
+    SECTION("square(scalar)") {
+        CHECK(op(TestType(0)) == 0);
+        CHECK(op(TestType(1)) == 1);
+        if constexpr (std::same_as<bool, TestType>) {
+            // square(bool) is identity
+        } else if constexpr (std::integral<TestType>) {
+            CHECK(op(TestType(3)) == 9);
+            CHECK(op(TestType(-3)) == 9);
+            CHECK(op(TestType(4)) == 16);
+        } else {  // floating
+            CHECK(op(TestType(2.5)) == 6.25);
+            CHECK(op(TestType(-1.5)) == 2.25);
+        }
+    }
+
+    SECTION("square(interval)") {
+        CHECK(not op(interval<TestType>()));  // op(empty) -> empty
+
+        if constexpr (std::same_as<bool, TestType>) {
+            CHECK(op(interval<bool>(0, 1)) == interval<bool>(0, 1));
+        } else {
+            CHECK(op(interval<TestType>(2, 3)) == interval(op(TestType(2)), op(TestType(3))));
+            CHECK(op(interval<TestType>(-3, -2)) == interval(op(TestType(-2)), op(TestType(-3))));
+            CHECK(op(interval<TestType>(-3, 2)) == interval(TestType(0), op(TestType(-3))));
+            CHECK(op(interval<TestType>(-2, 3)) == interval(TestType(0), op(TestType(3))));
+        }
+    }
+}
+
+TEMPLATE_LIST_TEST_CASE("square_root", "", DTypes) {
+    constexpr square_root op{};
+
+    SECTION("square_root(scalar)") {
+        CHECK(op(TestType(0)) == 0);
+        CHECK(op(TestType(1)) == 1);
+        if constexpr (not std::same_as<bool, TestType>) {
+            CHECK(op(TestType(4)) == 2);
+            CHECK(op(TestType(9)) == 3);
+            if constexpr (std::floating_point<TestType>) {
+                CHECK(op(TestType(2.0)) == std::sqrt(TestType(2.0)));
+            }
+        }
+    }
+
+    SECTION("square_root(interval)") {
+        CHECK(not op(interval<TestType>()));  // op(empty) -> empty
+        CHECK(op(interval<TestType>(0, 1)) == interval(op(TestType(0)), op(TestType(1))));
+        if constexpr (not std::same_as<bool, TestType>) {
+            CHECK(op(interval<TestType>(0, 4)) == interval(op(TestType(0)), op(TestType(4))));
+            CHECK(op(interval<TestType>(1, 9)) == interval(op(TestType(1)), op(TestType(9))));
+        }
+    }
+
+    SECTION("square_root domain is non-negative") {
+        CHECK(square_root::domain<TestType> == interval<TestType>::nonnegative());
+    }
+}
+
+TEMPLATE_LIST_TEST_CASE("tanh", "", DTypes) {
+    constexpr tanh op{};
+
+    SECTION("tanh(scalar)") {
+        CHECK(op(TestType(0)) == 0);  // tanh(0) == 0 exactly
+        if constexpr (not std::same_as<bool, TestType>) {
+            CHECK(op(TestType(1)) == std::tanh(TestType(1)));
+            CHECK(op(TestType(-2)) == std::tanh(TestType(-2)));
+        }
+    }
+
+    SECTION("tanh(interval)") {
+        CHECK(not op(interval<TestType>()));  // op(empty) -> empty
+        CHECK(op(interval<TestType>(0, 1)) == interval(op(TestType(0)), op(TestType(1))));
+        if constexpr (not std::same_as<bool, TestType>) {
+            CHECK(op(interval<TestType>(-2, 3)) == interval(op(TestType(-2)), op(TestType(3))));
+        }
+    }
 }
 
 }  // namespace dwave::optimization::functional