From 38ad64deeeaf3e9cb452af4a6c5acf9b935af19b Mon Sep 17 00:00:00 2001
From: Alexander Condello <arcondello@gmail.com>
Date: Tue, 18 Aug 2026 22:11:36 -0700
Subject: [PATCH 01/10] Build Windows wheels with LLVM

---
 .circleci/config.yml                           | 9 +++++++--
 meson.build                                    | 8 +++++---
 releasenotes/notes/c++23-60fb26a150a004e8.yaml | 3 +++
 3 files changed, 15 insertions(+), 5 deletions(-)
 create mode 100644 releasenotes/notes/c++23-60fb26a150a004e8.yaml

diff --git a/.circleci/config.yml b/.circleci/config.yml
index 3693ac0e9..bc6d90179 100644
--- a/.circleci/config.yml
+++ b/.circleci/config.yml
@@ -7,6 +7,9 @@ orbs:
 commands:
   run-cibuildwheel:
     parameters:
+      cibw-config:  # environment variables
+        type: string
+        default: ""
       cibw-version:
         type: string
         default: 4.2.0  # latest as of August 2026
@@ -16,7 +19,7 @@ commands:
           shell: bash -eo pipefail
           command: |
             python -m pip install --user cibuildwheel==<< parameters.cibw-version >>
-            python -m cibuildwheel --output-dir dist
+            << parameters.cibw-config >> python -m cibuildwheel --output-dir dist
 
       - store_artifacts: &store-artifacts
           path: ./dist
@@ -122,7 +125,9 @@ jobs:
 
     steps:
       - checkout
-      - run-cibuildwheel
+      - run: choco install llvm -y
+      - run-cibuildwheel:
+          cibw-config: CIBW_ENVIRONMENT_PASS=CXX CXX=clang-cl
 
   python-sdist:
     docker:
diff --git a/meson.build b/meson.build
index b8882ec82..cc2c65c1a 100644
--- a/meson.build
+++ b/meson.build
@@ -19,12 +19,14 @@ project(
     ).stdout(),
 )
 
-# We want some debugging symbols
-add_project_arguments(['-g1'], language: ['cpp'])
-
 cpp = meson.get_compiler('cpp')
 py = import('python').find_installation(pure: false)
 
+# We want some debugging symbols
+if cpp.has_argument('-g1')
+  add_project_arguments('-g1', language: 'cpp')
+endif
+
 dwave_optimization_include = include_directories('dwave/optimization/include/')
 dwave_optimization_src = [
     'dwave/optimization/src/nodes/_checkpoints.cpp',
diff --git a/releasenotes/notes/c++23-60fb26a150a004e8.yaml b/releasenotes/notes/c++23-60fb26a150a004e8.yaml
new file mode 100644
index 000000000..c2b893b07
--- /dev/null
+++ b/releasenotes/notes/c++23-60fb26a150a004e8.yaml
@@ -0,0 +1,3 @@
+---
+other:
+  - Build Windows wheels with LLVM rather than MSVC.

From 4d548108cc392a27b677c445d888b4f0319e9322 Mon Sep 17 00:00:00 2001
From: Alexander Condello <arcondello@gmail.com>
Date: Tue, 18 Aug 2026 22:30:16 -0700
Subject: [PATCH 02/10] Compile C++ library with C++23

---
 meson.build                                    | 2 +-
 releasenotes/notes/c++23-60fb26a150a004e8.yaml | 2 ++
 2 files changed, 3 insertions(+), 1 deletion(-)

diff --git a/meson.build b/meson.build
index cc2c65c1a..1e64bf33c 100644
--- a/meson.build
+++ b/meson.build
@@ -2,7 +2,7 @@ project(
     'dwave-optimization',
     'c', 'cpp', 'cython',
     default_options: [
-        'cpp_std=c++20',
+        'cpp_std=c++23',
         'buildtype=release',
         'b_ndebug=if-release',  # add -DNDEBUG when building for "release", which we do by default
         'optimization=3',
diff --git a/releasenotes/notes/c++23-60fb26a150a004e8.yaml b/releasenotes/notes/c++23-60fb26a150a004e8.yaml
index c2b893b07..a76532cb7 100644
--- a/releasenotes/notes/c++23-60fb26a150a004e8.yaml
+++ b/releasenotes/notes/c++23-60fb26a150a004e8.yaml
@@ -1,3 +1,5 @@
 ---
+features:
+  - Compile C++ library with C++23.
 other:
   - Build Windows wheels with LLVM rather than MSVC.

From 2c93f656a7e6a2cf0986ebe039bc92a6a9dae27d Mon Sep 17 00:00:00 2001
From: Alexander Condello <arcondello@gmail.com>
Date: Tue, 18 Aug 2026 22:32:24 -0700
Subject: [PATCH 03/10] Use std::unreachable() rather than backport

---
 .../include/dwave-optimization/common.hpp           | 13 -------------
 .../include/dwave-optimization/iterators.hpp        |  3 ++-
 .../include/dwave-optimization/nodes/testing.hpp    |  9 +++++----
 dwave/optimization/src/nodes/binaryop.cpp           | 11 ++++++-----
 dwave/optimization/src/nodes/collections.cpp        |  2 +-
 dwave/optimization/src/nodes/creation.cpp           |  9 +++++----
 dwave/optimization/src/nodes/indexing.cpp           |  7 ++++---
 dwave/optimization/src/nodes/lambda.cpp             |  4 +++-
 dwave/optimization/src/nodes/lp.cpp                 |  3 ++-
 dwave/optimization/src/nodes/naryop.cpp             |  5 +++--
 dwave/optimization/src/nodes/numbers.cpp            |  8 ++++----
 dwave/optimization/src/nodes/unaryop.cpp            |  6 ++++--
 12 files changed, 39 insertions(+), 41 deletions(-)

diff --git a/dwave/optimization/include/dwave-optimization/common.hpp b/dwave/optimization/include/dwave-optimization/common.hpp
index bcbf78981..07344982e 100644
--- a/dwave/optimization/include/dwave-optimization/common.hpp
+++ b/dwave/optimization/include/dwave-optimization/common.hpp
@@ -32,16 +32,3 @@ typedef std::int64_t ssize_t;
 #include <sys/types.h>  // for ssize_t
 
 #endif
-
-namespace dwave::optimization {
-
-// backport unreachable from c++23
-[[noreturn]] inline void unreachable() {
-#if defined(_MSC_VER) && !defined(__clang__)  // MSVC
-    __assume(false);
-#else  // GCC, Clang
-    __builtin_unreachable();
-#endif
-}
-
-}  // namespace dwave::optimization
diff --git a/dwave/optimization/include/dwave-optimization/iterators.hpp b/dwave/optimization/include/dwave-optimization/iterators.hpp
index f547459fb..4ec2b545b 100644
--- a/dwave/optimization/include/dwave-optimization/iterators.hpp
+++ b/dwave/optimization/include/dwave-optimization/iterators.hpp
@@ -19,6 +19,7 @@
 #include <concepts>
 #include <memory>
 #include <span>
+#include <utility>
 
 #include "dwave-optimization/common.hpp"  // for ssize_t
 #include "dwave-optimization/typing.hpp"
@@ -461,7 +462,7 @@ requires(
                 case FormatCharacter::signedlonglong_:
                     return *static_cast<const signed long long*>(ptr);
             }
-            unreachable();
+            std::unreachable();
         }
     }
 
diff --git a/dwave/optimization/include/dwave-optimization/nodes/testing.hpp b/dwave/optimization/include/dwave-optimization/nodes/testing.hpp
index 3954ff67d..d079d00e2 100644
--- a/dwave/optimization/include/dwave-optimization/nodes/testing.hpp
+++ b/dwave/optimization/include/dwave-optimization/nodes/testing.hpp
@@ -16,6 +16,7 @@
 
 #include <optional>
 #include <random>
+#include <utility>
 #include <vector>
 
 #include "dwave-optimization/array.hpp"
@@ -113,18 +114,18 @@ class DynamicArrayTestingNode : public ArrayOutputMixin<ArrayNode>, public Decis
         checkpoint_type& checkpoint
     ) const override {
         assert(false and "not implemented");
-        unreachable();
+        std::unreachable();
     }
     [[noreturn]] void assign_from_checkpoint(
         State& state,
         checkpoint_type&& checkpoint
     ) const override {
         assert(false and "not implemented");
-        unreachable();
+        std::unreachable();
     }
     [[noreturn]] virtual checkpoint_type checkpoint(State& state) const override {
         assert(false and "not implemented");
-        unreachable();
+        std::unreachable();
     }
 
     // State mutation methods *************************************************
@@ -185,7 +186,7 @@ class DynamicArrayTestingNode : public ArrayOutputMixin<ArrayNode>, public Decis
             case 2:
                 return set(state, rng);
             default:
-                unreachable();
+                std::unreachable();
         }
     }
 
diff --git a/dwave/optimization/src/nodes/binaryop.cpp b/dwave/optimization/src/nodes/binaryop.cpp
index 78ee6535e..5f3279b37 100644
--- a/dwave/optimization/src/nodes/binaryop.cpp
+++ b/dwave/optimization/src/nodes/binaryop.cpp
@@ -19,6 +19,7 @@
 #include <ranges>
 #include <stdexcept>
 #include <type_traits>
+#include <utility>
 #include <vector>
 
 #include "_state.hpp"
@@ -151,7 +152,7 @@ std::pair<double, double> calculate_values_minmax(const Array* lhs_ptr, const Ar
     }
 
     assert(false && "not implemeted yet");
-    unreachable();
+    std::unreachable();
 }
 
 template <class BinaryOp>
@@ -182,7 +183,7 @@ bool calculate_integral(const Array* lhs_ptr, const Array* rhs_ptr) {
     }
 
     assert(false && "not implemeted yet");
-    unreachable();
+    std::unreachable();
 }
 
 template <class BinaryOp>
@@ -294,7 +295,7 @@ void BinaryOpNode<BinaryOp>::initialize_state(State& state) const {
     } else {
         // this case is complicated we need to "stretch" dimensions into each other
         assert(false && "not yet implemented");
-        unreachable();
+        std::unreachable();
     }
 
     this->template emplace_data_ptr_<ArrayNodeStateData>(state, std::move(values));
@@ -427,7 +428,7 @@ void BinaryOpNode<BinaryOp>::propagate(State& state) const {
     } else {
         // this case is complicated we need to "stretch" dimensions into eachother
         assert(false && "not yet implemented");
-        unreachable();
+        std::unreachable();
     }
 
     if (ptr->diff().size()) Node::propagate(state);
@@ -512,7 +513,7 @@ SizeInfo binaryop_calculate_sizeinfo(
 
     // not possible for us to be dynamic and none of our predecessors to be
     assert(false && "not implemeted");
-    unreachable();
+    std::unreachable();
 }
 
 template <class BinaryOp>
diff --git a/dwave/optimization/src/nodes/collections.cpp b/dwave/optimization/src/nodes/collections.cpp
index bee598ca5..099553ff7 100644
--- a/dwave/optimization/src/nodes/collections.cpp
+++ b/dwave/optimization/src/nodes/collections.cpp
@@ -532,7 +532,7 @@ struct DisjointBitSetsNodeData_ : CheckpointableState, NodeStateData {
         }
 
         assert(false and "disjoint set elements must be in exactly one bit-set once");
-        unreachable();
+        std::unreachable();
     }
 
     void commit() {
diff --git a/dwave/optimization/src/nodes/creation.cpp b/dwave/optimization/src/nodes/creation.cpp
index 492ece3ab..9d7cf54ee 100644
--- a/dwave/optimization/src/nodes/creation.cpp
+++ b/dwave/optimization/src/nodes/creation.cpp
@@ -14,6 +14,7 @@
 
 #include "dwave-optimization/nodes/creation.hpp"
 
+#include <utility>
 #include <variant>
 
 #include "_state.hpp"
@@ -109,7 +110,7 @@ std::vector<double> arange(const ssize_t start, const ssize_t stop, const ssize_
         }
     } else {
         assert(false && "0 step not allowed");
-        unreachable();
+        std::unreachable();
     }
 
     return arange;
@@ -178,7 +179,7 @@ std::pair<double, double> calculate_values_minmax(
     }
 
     assert(false && "zero step not allowed");
-    unreachable();
+    std::unreachable();
 }
 
 const SizeInfo calculate_arange_sizeinfo(
@@ -226,7 +227,7 @@ const SizeInfo calculate_arange_sizeinfo(
         );
     } else {
         assert(false && "unreachable");
-        unreachable();
+        std::unreachable();
     }
 
     // Handles all cases EXCEPT the following: "Exactly one predecessor, it
@@ -423,7 +424,7 @@ void ARangeNode::propagate(State& state) const {
         }
     } else {
         assert(false && "zero step not allowed");
-        unreachable();
+        std::unreachable();
     }
 
     if (ptr->diff().size()) Node::propagate(state);
diff --git a/dwave/optimization/src/nodes/indexing.cpp b/dwave/optimization/src/nodes/indexing.cpp
index 879880749..47542a91b 100644
--- a/dwave/optimization/src/nodes/indexing.cpp
+++ b/dwave/optimization/src/nodes/indexing.cpp
@@ -17,6 +17,7 @@
 #include <algorithm>
 #include <ranges>
 #include <unordered_set>
+#include <utility>
 #include <variant>
 
 #include "_state.hpp"
@@ -781,7 +782,7 @@ void AdvancedIndexingNode::replace_predecessor_(ssize_t index, Node* node_ptr) {
     }
 
     assert(false and "should not be able to get here");
-    unreachable();
+    std::unreachable();
 }
 
 void AdvancedIndexingNode::revert(State& state) const {
@@ -1405,13 +1406,13 @@ SizeInfo basicindexing_calculate_sizeinfo(
         } else if (start < 0 && stop >= 0) {
             // -start:stop:step - imposes a nonlinear maximum size
             assert(false && "not linear, handled above");
-            unreachable();
+            std::unreachable();
         } else if (start < 0 && stop < 0) {
             // -start:-stop:step - imposes a maximum size
             sizeinfo.max = num_per_row * (std::max<ssize_t>(stop - start, 0) + step - 1) / step;
         } else {
             assert(false && "shouldn't be reachable");
-            unreachable();
+            std::unreachable();
         }
     }
 
diff --git a/dwave/optimization/src/nodes/lambda.cpp b/dwave/optimization/src/nodes/lambda.cpp
index b40ec9a55..75c875611 100644
--- a/dwave/optimization/src/nodes/lambda.cpp
+++ b/dwave/optimization/src/nodes/lambda.cpp
@@ -14,6 +14,8 @@
 
 #include "dwave-optimization/nodes/lambda.hpp"
 
+#include <utility>
+
 #include "_state.hpp"
 #include "dwave-optimization/array.hpp"
 #include "dwave-optimization/graph.hpp"
@@ -352,7 +354,7 @@ void AccumulateZipNode::propagate(State& state) const {
             data->emplace_back(val);
         } else {
             assert(false && "index is too large for current buffer");
-            unreachable();
+            std::unreachable();
         }
     }
 
diff --git a/dwave/optimization/src/nodes/lp.cpp b/dwave/optimization/src/nodes/lp.cpp
index b1f06f2c9..1ab60b839 100644
--- a/dwave/optimization/src/nodes/lp.cpp
+++ b/dwave/optimization/src/nodes/lp.cpp
@@ -15,6 +15,7 @@
 #include "dwave-optimization/nodes/lp.hpp"
 
 #include <string>
+#include <utility>
 
 #include "../simplex.hpp"
 #include "_state.hpp"
@@ -433,7 +434,7 @@ void LinearProgramNode::replace_predecessor_(ssize_t index, Node* node_ptr) {
     if (check_and_replace(lb_ptr_)) return;
     if (check_and_replace(ub_ptr_)) return;
 
-    unreachable();
+    std::unreachable();
     assert(false and "should never get here");
 }
 
diff --git a/dwave/optimization/src/nodes/naryop.cpp b/dwave/optimization/src/nodes/naryop.cpp
index 5dc47e482..150337f11 100644
--- a/dwave/optimization/src/nodes/naryop.cpp
+++ b/dwave/optimization/src/nodes/naryop.cpp
@@ -19,6 +19,7 @@
 #include <cmath>
 #include <stdexcept>
 #include <type_traits>
+#include <utility>
 
 #include "_state.hpp"
 
@@ -83,7 +84,7 @@ bool calculate_integral(const std::vector<Array*>& operands) {
     }
 
     assert(false && "not implemeted yet");
-    unreachable();
+    std::unreachable();
 }
 
 template <class BinaryOp>
@@ -142,7 +143,7 @@ ValuesInfo calculate_values_info(const std::vector<Array*>& operands) {
     }
 
     assert(false && "not implemeted yet");
-    unreachable();
+    std::unreachable();
 }
 
 template <class BinaryOp>
diff --git a/dwave/optimization/src/nodes/numbers.cpp b/dwave/optimization/src/nodes/numbers.cpp
index bed900251..a12d667cc 100644
--- a/dwave/optimization/src/nodes/numbers.cpp
+++ b/dwave/optimization/src/nodes/numbers.cpp
@@ -502,7 +502,7 @@ bool satisfies_sum_constraint(
                     break;
                 default:
                     assert(false && "Unexpected operator type.");
-                    unreachable();
+                    std::unreachable();
             }
         }
     }
@@ -593,7 +593,7 @@ double sum_constraint_delta(
             return (lhs < bound) ? (bound - lhs) : 0.0;
         default:
             assert(false && "Unexpected operator type.");
-            unreachable();
+            std::unreachable();
     }
 }
 
@@ -717,7 +717,7 @@ void NumberNode::initialize_state(State& state) const {
         initialize_state(state, std::move(values));
     } else {
         assert(false && "Multiple sum constraints not yet supported.");
-        unreachable();
+        std::unreachable();
     }
 }
 
@@ -1608,7 +1608,7 @@ void BinaryNode::initialize_state(State& state) const {
         initialize_state(state, std::move(values));
     } else {
         assert(false && "Multiple sum constraints not yet supported.");
-        unreachable();
+        std::unreachable();
     }
 }
 
diff --git a/dwave/optimization/src/nodes/unaryop.cpp b/dwave/optimization/src/nodes/unaryop.cpp
index 739027618..d871be910 100644
--- a/dwave/optimization/src/nodes/unaryop.cpp
+++ b/dwave/optimization/src/nodes/unaryop.cpp
@@ -14,6 +14,8 @@
 
 #include "dwave-optimization/nodes/unaryop.hpp"
 
+#include <utility>
+
 #include "_state.hpp"
 
 namespace dwave::optimization {
@@ -58,7 +60,7 @@ std::pair<double, double> calculate_values_minmax(const Array* array_ptr) {
             return std::make_pair(low, high);
         } else if (low >= 0) {
             assert(false && "min > max");
-            unreachable();
+            std::unreachable();
         } else if (high >= 0) {
             return std::pair<double, double>(0.0, std::max<double>(-low, high));
         } else {
@@ -99,7 +101,7 @@ std::pair<double, double> calculate_values_minmax(const Array* array_ptr) {
     }
 
     assert(false && "not implemeted yet");
-    unreachable();
+    std::unreachable();
 }
 
 template <class UnaryOp>

From 4d2a199a9840d32b45c8104ff74e336f1a9e2e34 Mon Sep 17 00:00:00 2001
From: Alexander Condello <arcondello@gmail.com>
Date: Tue, 18 Aug 2026 22:57:47 -0700
Subject: [PATCH 04/10] Use c++23 features in a few places

Specifically places where we have comments wishing for C++23
---
 .circleci/config.yml                           | 2 +-
 dwave/optimization/src/array.cpp               | 2 --
 dwave/optimization/src/nodes/flow.cpp          | 9 ++-------
 dwave/optimization/src/nodes/manipulation.cpp  | 9 +--------
 releasenotes/notes/c++23-60fb26a150a004e8.yaml | 2 ++
 5 files changed, 6 insertions(+), 18 deletions(-)

diff --git a/.circleci/config.yml b/.circleci/config.yml
index bc6d90179..a325f692b 100644
--- a/.circleci/config.yml
+++ b/.circleci/config.yml
@@ -412,7 +412,7 @@ workflows:
           name: cpp-gcc-<< matrix.gcc-version >>
           matrix:
             parameters:
-              gcc-version: ["12", "latest"]
+              gcc-version: ["14", "latest"]
       - cpp-macOS
       - serialization:
           requires:
diff --git a/dwave/optimization/src/array.cpp b/dwave/optimization/src/array.cpp
index 3484d30ca..5f897f4df 100644
--- a/dwave/optimization/src/array.cpp
+++ b/dwave/optimization/src/array.cpp
@@ -241,8 +241,6 @@ std::vector<ssize_t> broadcast_shapes(
     std::vector<ssize_t> shape(std::max(lhs.size(), rhs.size()));
 
     // Walk backwards through the shapes, checking for dimension compatibility.
-    // Technically span::rbegin() etc are c++23 features but it seems to work on
-    // all the compilers we care about. Whereas the c++20 ranges::rbegin() etc do not.
     auto lit = lhs.rbegin();
     const auto lend = lhs.rend();
     auto rit = rhs.rbegin();
diff --git a/dwave/optimization/src/nodes/flow.cpp b/dwave/optimization/src/nodes/flow.cpp
index 52dc5847c..0776410b6 100644
--- a/dwave/optimization/src/nodes/flow.cpp
+++ b/dwave/optimization/src/nodes/flow.cpp
@@ -301,13 +301,8 @@ void WhereNode::initialize_state(State& state) const {
         std::vector<double> values;
         values.reserve(condition.size());
 
-        // zip would be very nice here...
-        for (
-            auto cit = condition.begin(), xit = x.begin(), yit = y.begin();
-            cit != std::default_sentinel;
-            ++cit, ++xit, ++yit
-        ) {
-            values.emplace_back((*cit) ? *xit : *yit);
+        for (const auto& [ci, xi, yi] : std::views::zip(condition, x, y)) {
+            values.emplace_back(ci ? xi : yi);
         }
 
         emplace_data_ptr_<WhereNodeData>(state, std::move(values));
diff --git a/dwave/optimization/src/nodes/manipulation.cpp b/dwave/optimization/src/nodes/manipulation.cpp
index 1c45621ff..b4b2f0631 100644
--- a/dwave/optimization/src/nodes/manipulation.cpp
+++ b/dwave/optimization/src/nodes/manipulation.cpp
@@ -1193,14 +1193,7 @@ void ResizeNode::initialize_state(State& state) const {
     const ssize_t size = this->size();  // the desired size of our state
     assert(size >= 0);                  // we're never dynamic
 
-    std::vector<double> values;
-    values.reserve(size);
-
-    // Fill in from our predecessor, up to our size.
-    // In c++23 we could use append_range(...) which would be nicer.
-    for (const auto& v : array_ptr_->view(state) | std::views::take(size)) {
-        values.emplace_back(v);
-    }
+    auto values = std::ranges::to<std::vector>(array_ptr_->view(state) | std::views::take(size));
 
     // Now fill in everything else with our fill value
     assert(
diff --git a/releasenotes/notes/c++23-60fb26a150a004e8.yaml b/releasenotes/notes/c++23-60fb26a150a004e8.yaml
index a76532cb7..01b4b7bd2 100644
--- a/releasenotes/notes/c++23-60fb26a150a004e8.yaml
+++ b/releasenotes/notes/c++23-60fb26a150a004e8.yaml
@@ -1,5 +1,7 @@
 ---
 features:
   - Compile C++ library with C++23.
+upgrade:
+  - Drop support for GCC<14.
 other:
   - Build Windows wheels with LLVM rather than MSVC.

From dc664e052c5d559bf0d5215ced37d31361806814 Mon Sep 17 00:00:00 2001
From: Alexander Condello <arcondello@gmail.com>
Date: Thu, 20 Aug 2026 14:37:23 -0700
Subject: [PATCH 05/10] Update note about C++23 feature

---
 dwave/optimization/src/nodes/_checkpoints.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/dwave/optimization/src/nodes/_checkpoints.cpp b/dwave/optimization/src/nodes/_checkpoints.cpp
index c45bbfe80..cc3432088 100644
--- a/dwave/optimization/src/nodes/_checkpoints.cpp
+++ b/dwave/optimization/src/nodes/_checkpoints.cpp
@@ -80,7 +80,7 @@ void DiffCheckpoint::revert_updates(std::vector<Update> updates) {
 
     // We want to track the updates that would revert the changes from the
     // current state.
-    // In C++23 we could use assign_range() which would be nicer
+    // We'd like to use append_range(), but unfortunately gcc14 doesn't support it.
     auto relevant = std::move(updates) | std::views::take(drop_) | std::views::reverse |
                     std::views::transform([](const Update& up) { return up.inverse(); });
     updates_.emplace_back(relevant.begin(), relevant.end());

From 9196c6b897ba667859576ee64b5599c93046b0f4 Mon Sep 17 00:00:00 2001
From: Alexander Condello <arcondello@gmail.com>
Date: Thu, 20 Aug 2026 14:37:38 -0700
Subject: [PATCH 06/10] Mark additional Update methods as constexpr

---
 .../include/dwave-optimization/array.hpp      | 20 +++++++------------
 1 file changed, 7 insertions(+), 13 deletions(-)

diff --git a/dwave/optimization/include/dwave-optimization/array.hpp b/dwave/optimization/include/dwave-optimization/array.hpp
index dbfbacef5..e1251c726 100644
--- a/dwave/optimization/include/dwave-optimization/array.hpp
+++ b/dwave/optimization/include/dwave-optimization/array.hpp
@@ -326,29 +326,23 @@ struct Update {
     friend std::ostream& operator<<(std::ostream& os, const Update& update);
 
     // Whether the given index was placed when the state was grown
-    bool placed() const {
-        // We'd like to constexpr this, but std::isnan is not constexpr in C++20
-        return std::isnan(old);
-    }
+    constexpr bool placed() const { return std::isnan(old); }
 
     // Whether the given index was removed when the state was resized
-    bool removed() const {
-        // We'd like to constexpr this, but std::isnan is not constexpr in C++20
-        return std::isnan(value);
-    }
+    constexpr bool removed() const { return std::isnan(value); }
 
     // Returns true if the Update's goes from nothing to nothing (index can be anything)
-    bool null() const { return std::isnan(old) && std::isnan(value); }
+    constexpr bool null() const { return std::isnan(old) && std::isnan(value); }
 
-    double old_or(double val) const { return std::isnan(old) ? val : old; }
+    constexpr double old_or(double val) const { return std::isnan(old) ? val : old; }
 
-    double value_or(double val) const { return std::isnan(value) ? val : value; }
+    constexpr double value_or(double val) const { return std::isnan(value) ? val : value; }
 
     // Return true if the update does nothing - that is old and value are the same.
-    bool identity() const { return null() || old == value; }
+    constexpr bool identity() const { return null() || old == value; }
 
     // Return the update that would undo the current update
-    Update inverse() const { return Update(index, value, old); }
+    constexpr Update inverse() const { return Update(index, value, old); }
 
     // Use NaN to represent the "nothing" value used in placements/removals
     static constexpr double nothing = std::numeric_limits<double>::signaling_NaN();

From 0c58aeed1b8ca47bf5062c9446d61cec6b2a67b8 Mon Sep 17 00:00:00 2001
From: Alexander Condello <arcondello@gmail.com>
Date: Thu, 20 Aug 2026 15:06:51 -0700
Subject: [PATCH 07/10] Mark functional operator() as static

---
 .../include/dwave-optimization/functional.hpp | 32 +++++++++----------
 dwave/optimization/src/functional_.hpp        | 18 +++++------
 2 files changed, 25 insertions(+), 25 deletions(-)

diff --git a/dwave/optimization/include/dwave-optimization/functional.hpp b/dwave/optimization/include/dwave-optimization/functional.hpp
index 06928b5fc..e050a7c11 100644
--- a/dwave/optimization/include/dwave-optimization/functional.hpp
+++ b/dwave/optimization/include/dwave-optimization/functional.hpp
@@ -23,54 +23,54 @@ namespace dwave::optimization::functional {
 
 template <class T>
 struct abs {
-    constexpr T operator()(const T& x) const { return std::abs(x); }
+    static constexpr T operator()(const T& x) { return std::abs(x); }
 };
 
 template <class T>
 struct cos {
-    auto operator()(const T& num) const { return std::cos(num); }
+    static auto operator()(const T& num) { return std::cos(num); }
 };
 
 template <class T>
 struct exp {
-    constexpr auto operator()(const T& x) const { return std::exp(x); }
+    static constexpr auto operator()(const T& x) { return std::exp(x); }
 };
 
 template <class T>
 struct expit {
-    constexpr double operator()(const T& x) const { return 1.0 / (1.0 + std::exp(-1. * x)); }
+    static constexpr double operator()(const T& x) { return 1.0 / (1.0 + std::exp(-1. * x)); }
 };
 
 template <class T>
 struct log {
-    constexpr auto operator()(const T& x) const { return std::log(x); }
+    static constexpr auto operator()(const T& x) { return std::log(x); }
 };
 
 template <class T>
 struct logical {
-    constexpr bool operator()(const T& x) const { return x; }
+    static constexpr bool operator()(const T& x) { return x; }
 };
 
 template <class T>
 struct logical_xor {
-    constexpr bool operator()(const T& x, const T& y) const {
+    static constexpr bool operator()(const T& x, const T& y) {
         return static_cast<bool>(x) != static_cast<bool>(y);
     }
 };
 
 template <class T>
 struct max {
-    constexpr T operator()(const T& x, const T& y) const { return std::max(x, y); }
+    static constexpr T operator()(const T& x, const T& y) { return std::max(x, y); }
 };
 
 template <class T>
 struct min {
-    constexpr T operator()(const T& x, const T& y) const { return std::min(x, y); }
+    static constexpr T operator()(const T& x, const T& y) { return std::min(x, y); }
 };
 
 template <class T>
 struct modulus {
-    constexpr T operator()(const T& x, const T& y) const {
+    static constexpr T operator()(const T& x, const T& y) {
         // Copy numpy behavior and return 0 for `x % 0`
         if (y == 0) return 0;
 
@@ -92,12 +92,12 @@ struct modulus {
 
 template <class T>
 struct rint {
-    constexpr auto operator()(const T& x) const { return std::rint(x); }
+    static constexpr auto operator()(const T& x) { return std::rint(x); }
 };
 
 template <class T>
 struct safe_divides {
-    constexpr T operator()(const T& lhs, const T& rhs) const {
+    static constexpr T operator()(const T& lhs, const T& rhs) {
         if (!rhs) return 0;
         return lhs / rhs;
     }
@@ -105,22 +105,22 @@ struct safe_divides {
 
 template <class T>
 struct sin {
-    auto operator()(const T& num) const { return std::sin(num); }
+    static auto operator()(const T& num) { return std::sin(num); }
 };
 
 template <class T>
 struct square {
-    constexpr T operator()(const T& x) const { return x * x; }
+    static constexpr T operator()(const T& x) { return x * x; }
 };
 
 template <class T>
 struct square_root {
-    constexpr auto operator()(const T& x) const { return std::sqrt(x); }
+    static constexpr auto operator()(const T& x) { return std::sqrt(x); }
 };
 
 template <class T>
 struct tanh {
-    auto operator()(const T& num) const { return std::tanh(num); }
+    static auto operator()(const T& num) { return std::tanh(num); }
 };
 
 }  // namespace dwave::optimization::functional
diff --git a/dwave/optimization/src/functional_.hpp b/dwave/optimization/src/functional_.hpp
index e1ad70821..18d311299 100644
--- a/dwave/optimization/src/functional_.hpp
+++ b/dwave/optimization/src/functional_.hpp
@@ -107,7 +107,7 @@ struct Add : BinaryFunctionMixin<Add<double>> {
     ///
     /// The `rhs` may be `result_type`. Some binary operations might support
     /// additional types.
-    result_type operator()(const DType auto& lhs, const DType auto& rhs) const noexcept {
+    static result_type operator()(const DType auto& lhs, const DType auto& rhs) noexcept {
         return lhs + rhs;
     }
 
@@ -199,10 +199,10 @@ struct LogicalAnd : BinaryFunctionMixin<LogicalAnd> {
 
     /// @brief Return the logical and of `lhs` and `rhs`.
     /// @copydetails Add::operator()
-    result_type operator()(const DType auto& lhs, const DType auto& rhs) const noexcept {
+    static result_type operator()(const DType auto& lhs, const DType auto& rhs) noexcept {
         return lhs and rhs;
     }
-    reduction_type operator()(reduction_type lhs, const DType auto& rhs) const noexcept {
+    static reduction_type operator()(reduction_type lhs, const DType auto& rhs) noexcept {
         if (rhs == 0) lhs.num_falsy_ += 1;
         return lhs;
     }
@@ -263,10 +263,10 @@ struct LogicalOr : BinaryFunctionMixin<LogicalOr> {
 
     /// @brief Return the logical or of `lhs` and `rhs`.
     /// @copydetails Add::operator()
-    result_type operator()(const DType auto& lhs, const DType auto& rhs) const noexcept {
+    static result_type operator()(const DType auto& lhs, const DType auto& rhs) noexcept {
         return lhs or rhs;
     }
-    reduction_type operator()(reduction_type lhs, const DType auto& rhs) const noexcept {
+    static reduction_type operator()(reduction_type lhs, const DType auto& rhs) noexcept {
         lhs.num_truthy_ += (rhs != 0);
         return lhs;
     }
@@ -315,7 +315,7 @@ struct Maximum : BinaryFunctionMixin<Maximum<T>> {
 
     /// @brief Return the max of `lhs` and `rhs`.
     /// @copydetails Add::operator()
-    result_type operator()(const DType auto& lhs, const DType auto& rhs) const noexcept {
+    static result_type operator()(const DType auto& lhs, const DType auto& rhs) noexcept {
         return std::max<result_type>(lhs, rhs);
     }
 
@@ -359,7 +359,7 @@ struct Minimum : BinaryFunctionMixin<Minimum<T>> {
 
     /// @brief Return the min of `lhs` and `rhs`.
     /// @copydetails Add::operator()
-    result_type operator()(const DType auto& lhs, const DType auto& rhs) const noexcept {
+    static result_type operator()(const DType auto& lhs, const DType auto& rhs) noexcept {
         return std::min<result_type>(lhs, rhs);
     }
 
@@ -419,10 +419,10 @@ struct Multiply : BinaryFunctionMixin<Multiply<T>> {
 
     /// @brief Return the product of `lhs` and `rhs`.
     /// @copydetails Add::operator()
-    result_type operator()(const DType auto& lhs, const DType auto& rhs) const noexcept {
+    static result_type operator()(const DType auto& lhs, const DType auto& rhs) noexcept {
         return lhs * rhs;
     }
-    reduction_type operator()(reduction_type lhs, const DType auto& rhs) const noexcept {
+    static reduction_type operator()(reduction_type lhs, const DType auto& rhs) noexcept {
         if (rhs == 0) {
             lhs.num_zero_ += 1;
         } else {

From 26c500ea4c30305e0fdb35aebe2752e44cfab4a4 Mon Sep 17 00:00:00 2001
From: Alexander Condello <arcondello@gmail.com>
Date: Thu, 20 Aug 2026 15:15:24 -0700
Subject: [PATCH 08/10] Update operator[] instead of operator() in simplex.cpp

---
 dwave/optimization/src/simplex.cpp | 74 +++++++++++++-----------------
 1 file changed, 33 insertions(+), 41 deletions(-)

diff --git a/dwave/optimization/src/simplex.cpp b/dwave/optimization/src/simplex.cpp
index 29c8060d0..14e2fc34c 100644
--- a/dwave/optimization/src/simplex.cpp
+++ b/dwave/optimization/src/simplex.cpp
@@ -35,20 +35,12 @@ class Matrix {
     Matrix(std::span<const double> data, ssize_t n, ssize_t m) :
         Matrix(std::vector<double>(data.begin(), data.end()), n, m) {}
 
-    double& operator()(ssize_t i, ssize_t j) {
-        if (i < 0) i += n_;
-        if (j < 0) j += m_;
-        assert(i >= 0 && i < n_);
-        assert(j >= 0 && j < m_);
-        return buffer_[i * m_ + j];
-    }
-
-    const double& operator()(ssize_t i, ssize_t j) const {
-        if (i < 0) i += n_;
-        if (j < 0) j += m_;
-        assert(i >= 0 && i < n_);
-        assert(j >= 0 && j < m_);
-        return buffer_[i * m_ + j];
+    auto&& operator[](this auto&& self, ssize_t i, ssize_t j) {
+        if (i < 0) i += self.n_;
+        if (j < 0) j += self.m_;
+        assert(i >= 0 && i < self.n_);
+        assert(j >= 0 && j < self.m_);
+        return self.buffer_[i * self.m_ + j];
     }
 
     ssize_t n() const { return n_; }
@@ -67,7 +59,7 @@ ssize_t _pivot_col(Matrix& T, double tolerance, bool bland) {
     double min_val = std::numeric_limits<double>::infinity();
     ssize_t col = -1;
     for (ssize_t j = 0; j < T.m() - 1; j++) {
-        double val = T(-1, j);
+        double val = T[-1, j];
         if (val < -tolerance) {
             // Bland's rule: return the first index with negative value
             if (bland) return j;
@@ -95,9 +87,9 @@ ssize_t _pivot_row(
     double min_q = std::numeric_limits<double>::infinity();
     std::vector<ssize_t> min_rows;
     for (ssize_t i = 0; i < T.n() - k; i++) {
-        if (T(i, pivcol) <= tolerance) continue;
+        if (T[i, pivcol] <= tolerance) continue;
 
-        double q = T(i, -1) / T(i, pivcol);
+        double q = T[i, -1] / T[i, pivcol];
         if (q < min_q) {
             min_rows.clear();
             min_rows.push_back(i);
@@ -132,18 +124,18 @@ void _apply_pivot(
 ) {
     basis[pivrow] = pivcol;
 
-    double pivval = T(pivrow, pivcol);
+    double pivval = T[pivrow, pivcol];
 
     for (ssize_t j = 0; j < T.m(); j++) {
-        T(pivrow, j) /= pivval;
+        T[pivrow, j] /= pivval;
     }
 
     for (ssize_t irow = 0; irow < T.n(); irow++) {
         if (irow == pivrow) continue;
 
-        double val = T(irow, pivcol);
+        double val = T[irow, pivcol];
         for (ssize_t j = 0; j < T.m(); j++) {
-            T(irow, j) -= T(pivrow, j) * val;
+            T[irow, j] -= T[pivrow, j] * val;
         }
     }
 }
@@ -167,7 +159,7 @@ SolveResult _solve_simplex(
             if (basis[pivrow] <= T.m() - 2) continue;
 
             for (ssize_t col = 0; col < T.m() - 1; col++) {
-                if (std::abs(T(pivrow, col)) > tolerance) {
+                if (std::abs(T[pivrow, col]) > tolerance) {
                     _apply_pivot(T, basis, pivrow, col, tolerance);
                     status.num_iterations++;
                     break;
@@ -225,31 +217,31 @@ Matrix construct_T(
     // Copy in A to T[:n, :m] and b.T to T[:, -1]
     for (ssize_t i = 0; i < A.n(); i++) {
         double sign = b[i] < 0 ? -1.0 : 1.0;
-        T(i, -1) = sign * b[i];
+        T[i, -1] = sign * b[i];
         for (ssize_t j = 0; j < A.m(); j++) {
-            T(i, j) = sign * A(i, j);
+            T[i, j] = sign * A[i, j];
         }
     }
 
     // T[:, m:m+n] = I
     for (ssize_t i = 0; i < A.n(); i++) {
-        T(i, A.m() + i) = 1;
+        T[i, A.m() + i] = 1;
     }
 
     // Row objective
     for (ssize_t i = 0; i < A.m(); i++) {
-        T(A.n(), i) = c[i];
+        T[A.n(), i] = c[i];
     }
-    T(A.n(), -1) = c0;
+    T[A.n(), -1] = c0;
 
     // Row pseudo objective
     for (ssize_t j = 0; j < A.m(); j++) {
         for (ssize_t i = 0; i < A.n(); i++) {
-            T(A.n() + 1, j) -= T(i, j);
+            T[A.n() + 1, j] -= T[i, j];
         }
     }
     for (ssize_t i = 0; i < A.n(); i++) {
-        T(-1, -1) -= T(i, -1);
+        T[-1, -1] -= T[i, -1];
     }
 
     return T;
@@ -277,18 +269,18 @@ SolveResult _linprog_simplex(
     SolveResult status =
         _solve_simplex(T, A.n(), basis, max_iterations, tolerance, phase, bland, 0);
 
-    if (std::abs(T(-1, -1)) < tolerance) {
+    if (std::abs(T[-1, -1]) < tolerance) {
         Matrix newT(A.n() + 1, A.m() + 1);
 
         // newT[:, :m] = T[:-1, :m]
         for (ssize_t i = 0; i < newT.n(); i++) {
             for (ssize_t j = 0; j < A.m(); j++) {
-                newT(i, j) = T(i, j);
+                newT[i, j] = T[i, j];
             }
         }
         // newT[:, -1] = T[:, -1]
         for (ssize_t i = 0; i < newT.n(); i++) {
-            newT(i, -1) = T(i, -1);
+            newT[i, -1] = T[i, -1];
         }
 
         std::swap(T, newT);
@@ -311,7 +303,7 @@ SolveResult _linprog_simplex(
     std::vector<double> solution(A.m(), 0.0);
     for (ssize_t i = 0; i < A.n(); i++) {
         if (basis[i] < A.m()) {
-            solution[basis[i]] = T(i, -1);
+            solution[basis[i]] = T[i, -1];
         }
     }
     result.set_partial_solution(std::move(solution));
@@ -415,7 +407,7 @@ LP translate_LP_to_simple(
         if (!lb_is_unbounded(b_lb[i])) {
             // Copy the flipped constraint
             for (ssize_t j = 0; j < A.m(); j++) {
-                A_(A_row, j) = -A(i, j);
+                A_[A_row, j] = -A[i, j];
             }
             // Copy the flipped bound
             b[A_row] = -b_lb[i];
@@ -425,7 +417,7 @@ LP translate_LP_to_simple(
         if (!ub_is_unbounded(b_ub[i])) {
             // Copy the constraint
             for (ssize_t j = 0; j < A.m(); j++) {
-                A_(A_row, j) = A(i, j);
+                A_[A_row, j] = A[i, j];
             }
             // Copy the bound
             b[A_row] = b_ub[i];
@@ -437,7 +429,7 @@ LP translate_LP_to_simple(
     for (ssize_t i = 0; i < A_eq.n(); i++) {
         ssize_t A_row = i + A_constraint_count + upper_bounded_var_count;
         for (ssize_t j = 0; j < A_eq.m(); j++) {
-            A_(A_row, j) = A_eq(i, j);
+            A_[A_row, j] = A_eq[i, j];
         }
         b[A_row] = b_eq[i];
     }
@@ -449,7 +441,7 @@ LP translate_LP_to_simple(
         if (lb_is_unbounded(lb_[j]) && ub_is_unbounded(ub_[j])) {
             // Free variable, substitute xi = xi+ - xi-
             for (ssize_t i = 0; i < A_.n(); i++) {
-                A_(i, free_variable_index) = -A_(i, j);
+                A_[i, free_variable_index] = -A_[i, j];
             }
             c_[free_variable_index] = -c[j];
             free_variable_index++;
@@ -461,12 +453,12 @@ LP translate_LP_to_simple(
             c_[j] = -c_[j];
 
             for (ssize_t i = 0; i < A_.n(); i++) {
-                A_(i, j) *= -1;
+                A_[i, j] *= -1;
             }
         }
 
         if (!ub_is_unbounded(ub_[j])) {
-            A_(A_row, j) = 1;
+            A_[A_row, j] = 1;
             b[A_row] = ub_[j];
             A_row++;
         }
@@ -477,7 +469,7 @@ LP translate_LP_to_simple(
 
     // Add slack variables for inequalities
     for (ssize_t i = 0; i < slack_var_count; i++) {
-        A_(i, free_variable_index + i) = 1;
+        A_[i, free_variable_index + i] = 1;
     }
 
     // Substitute the lower bounds
@@ -487,7 +479,7 @@ LP translate_LP_to_simple(
             c0 += lb_[j] * c[j];
 
             for (ssize_t i = 0; i < A_.n(); i++) {
-                b[i] -= A_(i, j) * lb_[j];
+                b[i] -= A_[i, j] * lb_[j];
             }
         }
     }

From b5c619e336a7b780ab540c675c0f3b51e02f9322 Mon Sep 17 00:00:00 2001
From: Alexander Condello <arcondello@gmail.com>
Date: Thu, 20 Aug 2026 15:30:44 -0700
Subject: [PATCH 09/10] Use std::views::zip for clarity

---
 dwave/optimization/src/nodes/binaryop.cpp | 6 ++----
 1 file changed, 2 insertions(+), 4 deletions(-)

diff --git a/dwave/optimization/src/nodes/binaryop.cpp b/dwave/optimization/src/nodes/binaryop.cpp
index 5f3279b37..5384781a1 100644
--- a/dwave/optimization/src/nodes/binaryop.cpp
+++ b/dwave/optimization/src/nodes/binaryop.cpp
@@ -268,10 +268,8 @@ void BinaryOpNode<BinaryOp>::initialize_state(State& state) const {
         // This is the easy case - all we need to do is iterate over both as flat arrays
         values.reserve(lhs_ptr->size(state));
 
-        auto it = lhs_ptr->begin(state);
-        for (const double val : rhs_ptr->view(state)) {
-            values.emplace_back(op(*it, val));  // order is important
-            ++it;
+        for (const auto& [x, y] : std::views::zip(lhs_ptr->view(state), rhs_ptr->view(state))) {
+            values.emplace_back(op(x, y));
         }
 
     } else if (lhs_ptr->size() == 1) {

From ea9442e8438ab9e23e304f78e573a08650a11846 Mon Sep 17 00:00:00 2001
From: Alexander Condello <arcondello@gmail.com>
Date: Thu, 20 Aug 2026 15:40:42 -0700
Subject: [PATCH 10/10] Use std::ranges::to for adapting ranges

---
 .../dwave-optimization/nodes/manipulation.hpp    | 16 ++++++++++------
 .../include/dwave-optimization/nodes/numbers.hpp | 12 ++++++++----
 dwave/optimization/src/nodes/_state.hpp          | 12 ++++++------
 dwave/optimization/src/nodes/flow.cpp            |  7 ++++---
 4 files changed, 28 insertions(+), 19 deletions(-)

diff --git a/dwave/optimization/include/dwave-optimization/nodes/manipulation.hpp b/dwave/optimization/include/dwave-optimization/nodes/manipulation.hpp
index 21cfc549e..1b4fab40e 100644
--- a/dwave/optimization/include/dwave-optimization/nodes/manipulation.hpp
+++ b/dwave/optimization/include/dwave-optimization/nodes/manipulation.hpp
@@ -270,9 +270,9 @@ class ReshapeNode : public ArrayOutputMixin<EqualityMixin<ArrayNode, ReshapeNode
     ///
     /// @param array_ptr The array to be reshaped.
     /// @param shape The new shape. Must have the same size as the original shape.
-    template <std::ranges::range Range>
-    ReshapeNode(ArrayNode* node_ptr, Range&& shape) :
-        ReshapeNode(node_ptr, std::vector<ssize_t>(shape.begin(), shape.end())) {}
+    template <std::ranges::range R>
+    ReshapeNode(ArrayNode* node_ptr, R&& shape) :
+        ReshapeNode(node_ptr, std::ranges::to<std::vector<ssize_t>>(std::forward<R>(shape))) {}
 
     /// @copydoc Array::buff()
     double const* buff(const State& state) const override;
@@ -348,9 +348,13 @@ class ResizeNode : public ArrayOutputMixin<EqualityMixin<ArrayNode, ResizeNode>>
     /// @param array_ptr The array to be resized.
     /// @param shape The new shape. Must not be dynamic.
     /// @param fill_value The value to use for missing values.
-    template <std::ranges::range Range>
-    ResizeNode(ArrayNode* node_ptr, Range&& shape, double fill_value = 0) :
-        ResizeNode(node_ptr, std::vector<ssize_t>(shape.begin(), shape.end()), fill_value) {}
+    template <std::ranges::range R>
+    ResizeNode(ArrayNode* node_ptr, R&& shape, double fill_value = 0) :
+        ResizeNode(
+            node_ptr,
+            std::ranges::to<std::vector<ssize_t>>(std::forward<R>(shape)),
+            fill_value
+        ) {}
 
     /// @copydoc Array::buff()
     double const* buff(const State& state) const override;
diff --git a/dwave/optimization/include/dwave-optimization/nodes/numbers.hpp b/dwave/optimization/include/dwave-optimization/nodes/numbers.hpp
index ea4c18f25..21155753d 100644
--- a/dwave/optimization/include/dwave-optimization/nodes/numbers.hpp
+++ b/dwave/optimization/include/dwave-optimization/nodes/numbers.hpp
@@ -112,8 +112,10 @@ class NumberNode : public ArrayOutputMixin<ArrayNode>, public DecisionNode {
 
     // Initialize a state from an existing container, making a copy.
     template <std::ranges::range R>
-    void initialize_state(State& state, const R& values) const {
-        return initialize_state(state, std::vector<double>(values.begin(), values.end()));
+    void initialize_state(State& state, R&& values) const {
+        return initialize_state(
+            state, std::ranges::to<std::vector<double>>(std::forward<R>(values))
+        );
     }
 
     /// @copydoc Node::propagate()
@@ -410,8 +412,10 @@ class BinaryNode : public IntegerNode {
 
     /// Initialize a state from an existing container, making a copy.
     template <std::ranges::range R>
-    void initialize_state(State& state, const R& values) const {
-        return initialize_state(state, std::vector<double>(values.begin(), values.end()));
+    void initialize_state(State& state, R&& values) const {
+        return initialize_state(
+            state, std::ranges::to<std::vector<double>>(std::forward<R>(values))
+        );
     }
 
     /// ************************** BinaryNode methods **************************
diff --git a/dwave/optimization/src/nodes/_state.hpp b/dwave/optimization/src/nodes/_state.hpp
index 0ba4e597c..ff1105c0c 100644
--- a/dwave/optimization/src/nodes/_state.hpp
+++ b/dwave/optimization/src/nodes/_state.hpp
@@ -35,9 +35,9 @@ class ArrayStateData {
     explicit ArrayStateData(std::vector<double>&& values) noexcept :
         buffer(std::move(values)), size_(buffer.size()), previous_size_(buffer.size()) {}
 
-    template <std::ranges::range Range>
-    explicit ArrayStateData(Range&& values) noexcept :
-        ArrayStateData(std::vector<double>(values.begin(), values.end())) {}
+    template <std::ranges::range R>
+    explicit ArrayStateData(R&& values) noexcept :
+        ArrayStateData(std::ranges::to<std::vector<double>>(std::forward<R>(values))) {}
 
     // Assign new values to the state starting from an offset, tracking the changes from the
     // previous state to the new. If the original buffer extends past the new range of values,
@@ -272,9 +272,9 @@ class ArrayNodeStateData : public ArrayStateData, public NodeStateData {
     explicit ArrayNodeStateData(std::vector<double>&& values) noexcept :
         ArrayStateData(std::move(values)), NodeStateData() {}
 
-    template <std::ranges::random_access_range Range>
-    explicit ArrayNodeStateData(Range&& values) noexcept :
-        ArrayNodeStateData(std::vector<double>(values.begin(), values.begin() + values.size())) {}
+    template <std::ranges::range R>
+    explicit ArrayNodeStateData(R&& values) noexcept :
+        ArrayNodeStateData(std::ranges::to<std::vector<double>>(std::forward<R>(values))) {}
 
     std::unique_ptr<NodeStateData> copy() const override {
         return std::make_unique<ArrayNodeStateData>(*this);
diff --git a/dwave/optimization/src/nodes/flow.cpp b/dwave/optimization/src/nodes/flow.cpp
index 0776410b6..9095c3eeb 100644
--- a/dwave/optimization/src/nodes/flow.cpp
+++ b/dwave/optimization/src/nodes/flow.cpp
@@ -178,12 +178,13 @@ SizeInfo ExtractNode::sizeinfo() const { return this->sizeinfo_; }
 
 struct WhereNodeData : ArrayNodeStateData {
     // Initialize the state with the values given
-    explicit WhereNodeData(const std::ranges::view auto& values) noexcept :
-        ArrayNodeStateData(std::vector<double>(values.begin(), values.begin() + values.size())) {}
-
     explicit WhereNodeData(std::vector<double>&& values) noexcept :
         ArrayNodeStateData(std::move(values)) {}
 
+    template<std::ranges::range R>
+    explicit WhereNodeData(R&& values) noexcept :
+        ArrayNodeStateData(std::ranges::to<std::vector<double>>(std::forward<R>(values))) {}
+
     // Update the buffer according to the given diffs
     void apply_diffs(
         std::ranges::view auto&& condition,