// SPDX-License-Identifier: MPL-2.0
#pragma once

// Runtime-rank CPU contract. Mapping/ownership separation follows mdspan,
// doi:10.1109/P3HPC49587.2019.00011; dependent domains follow TAMM III-IV,
// doi:10.1063/5.0142433. Neither source supplies this runtime's lifetime policy.
#include <algorithm>
#include <array>
#include <atomic>
#include <complex>
#include <cstddef>
#include <cstdint>
#include <limits>
#include <memory>
#include <mutex>
#include <stdexcept>
#include <type_traits>
#include <utility>
#include <vector>

namespace vibeqc::tensor {

inline std::size_t add(std::size_t a, std::size_t b) {
    if (b > std::numeric_limits<std::size_t>::max() - a)
        throw std::overflow_error("tensor size addition overflow");
    return a + b;
}
inline std::size_t mul(std::size_t a, std::size_t b) {
    if (a && b > std::numeric_limits<std::size_t>::max() / a)
        throw std::overflow_error("tensor size multiplication overflow");
    return a * b;
}
inline std::ptrdiff_t signed_size(std::size_t n) {
    if (n > static_cast<std::size_t>(PTRDIFF_MAX))
        throw std::overflow_error("tensor signed address overflow");
    return static_cast<std::ptrdiff_t>(n);
}
inline std::ptrdiff_t signed_add(std::ptrdiff_t a, std::ptrdiff_t b) {
    if ((b > 0 && a > PTRDIFF_MAX - b) || (b < 0 && a < PTRDIFF_MIN - b))
        throw std::overflow_error("tensor offset addition overflow");
    return a + b;
}
inline std::ptrdiff_t stride_span(std::size_t n, std::ptrdiff_t stride) {
    if (!n || !stride) return 0;
    // Do not negate PTRDIFF_MIN.
    auto magnitude = stride < 0 ? std::size_t(-(stride + 1)) + 1 : std::size_t(stride);
    auto span = mul(n, magnitude);
    if (stride < 0 && span == std::size_t(PTRDIFF_MAX) + 1) return PTRDIFF_MIN;
    auto result = signed_size(span);
    return stride < 0 ? -result : result;
}

template<class T> std::shared_ptr<T> checked_handle(std::shared_ptr<T> handle) {
    if (!handle) throw std::invalid_argument("null tensor handle");
    return handle;
}

// A tracked byte ceiling, not an allocator/RSS guarantee. Reservations include
// fixed logical control objects and requested payloads. OS, allocator and
// OpenMP runtime overhead are outside the boundary (see the CPU ADR).
class Budget {
    struct State {
        std::size_t cap, live = 0, peak = 0;
        std::mutex mutex;
        explicit State(std::size_t n) : cap(n) {
            live = sizeof(State) + sizeof(Budget);
            if (n < live) throw std::invalid_argument("tensor budget cannot hold its ledger");
            peak = live;
        }
    };
    std::shared_ptr<State> state_;
public:
    class Reservation {
        std::shared_ptr<State> state_;
        std::size_t bytes_ = 0;
        friend class Budget;
        Reservation(std::shared_ptr<State> s, std::size_t n) : state_(std::move(s)), bytes_(n) {}
    public:
        Reservation() = default;
        Reservation(const Reservation&) = delete;
        Reservation& operator=(const Reservation&) = delete;
        Reservation(Reservation&& x) noexcept : state_(std::move(x.state_)), bytes_(x.bytes_) {}
        Reservation& operator=(Reservation&& x) noexcept {
            if (this != &x) { release(); state_ = std::move(x.state_); bytes_ = x.bytes_; }
            return *this;
        }
        ~Reservation() { release(); }
        void release() noexcept {
            if (state_) {
                { std::lock_guard<std::mutex> lock(state_->mutex); state_->live -= bytes_; }
                state_.reset();
            }
        }
    };
    explicit Budget(std::size_t cap) : state_(std::make_shared<State>(cap)) {
        if (!cap) throw std::invalid_argument("tensor budget must be explicit and positive");
    }
    Reservation reserve(std::size_t bytes) const {
        std::lock_guard<std::mutex> lock(state_->mutex);
        if (bytes > state_->cap - state_->live)
            throw std::runtime_error("tensor tracked memory budget exceeded");
        state_->live += bytes;
        state_->peak = std::max(state_->peak, state_->live);
        return Reservation(state_, bytes);
    }
    bool same(const Budget& b) const { return state_ == b.state_; }
    static constexpr std::size_t baseline_bytes() { return sizeof(State) + sizeof(Budget); }
    std::size_t live() const { std::lock_guard<std::mutex> lock(state_->mutex); return state_->live; }
    std::size_t limit() const { return state_->cap; }
    std::size_t available() const { std::lock_guard<std::mutex> lock(state_->mutex); return state_->cap-state_->live; }
    std::size_t peak() const { std::lock_guard<std::mutex> lock(state_->mutex); return state_->peak; }
};

// Fixed-size arrays avoid vector capacity growth and make metadata admission
// precede allocation. The containing object's fixed bytes are charged by its owner.
template<class T> class Array {
    Budget::Reservation reservation_;
    std::unique_ptr<T[]> values_;
    std::size_t size_ = 0;
public:
    Array(const Budget& b, std::size_t n)
        : reservation_(b.reserve(mul(n, sizeof(T)))), values_(n ? new T[n] : nullptr), size_(n) {}
    Array(Array&&) = default;
    Array& operator=(Array&&) = default;
    T& operator[](std::size_t i) { return values_[i]; }
    const T& operator[](std::size_t i) const { return values_[i]; }
    T* data() { return values_.get(); }
    const T* data() const { return values_.get(); }
    std::size_t size() const { return size_; }
};

class Space {
    std::uint64_t id_ = 0, revision_ = 0;
    inline static std::atomic<std::uint64_t> next_{1};
public:
    static Space create() {
        // Saturation prevents identity reuse, even after an overflow attempt.
        auto n = next_.load();
        do {
            if (n == UINT64_MAX) throw std::overflow_error("tensor identity exhausted");
        } while (!next_.compare_exchange_weak(n, n + 1));
        Space s; s.id_ = n; return s;
    }
    Space revised() const {
        if (!id_ || revision_ == UINT64_MAX) throw std::overflow_error("tensor revision exhausted");
        Space s = *this; ++s.revision_; return s;
    }
    bool valid() const { return id_ != 0; }
    bool operator==(const Space& s) const { return id_ == s.id_ && revision_ == s.revision_; }
    bool operator!=(const Space& s) const { return !(*this == s); }
};

struct Domain {
    Space space, parent;
    std::size_t parent_coordinate = 0, extent = 0, begin = 0, step = 1;
    // Selected maps are immutable and retained by every axis/plan that uses them.
    struct Selection {
        Budget budget;
        Budget::Reservation control;
        const Space identity = Space::create();
        Array<std::size_t> coordinates;
        Selection(const Budget& b, const std::vector<std::size_t>& c)
            : budget(b), control(b.reserve(sizeof(Selection))), coordinates(b, c.size()) {
            for (std::size_t i = 0; i < c.size(); ++i) coordinates[i] = c[i];
        }
    };
    std::shared_ptr<const Selection> selection;
    static Domain range(Space s, std::size_t n, std::size_t first = 0, std::size_t increment = 1,
                        Space p = {}, std::size_t key = 0) {
        if (!s.valid() || (!increment && n > 1)) throw std::invalid_argument("invalid tensor domain");
        if (n) add(first, mul(n - 1, increment));
        return {s, p, key, n, first, increment, {}};
    }
    static Domain selected(const Budget& b, Space s, const std::vector<std::size_t>& coordinates,
                           Space p = {}, std::size_t key = 0) {
        Domain d = range(s, coordinates.size(), 0, 1, p, key);
        d.selection = std::make_shared<Selection>(b, coordinates);
        return d;
    }
    void validate(const Budget& b) const {
        if (!space.valid() || (!step && extent > 1)) throw std::invalid_argument("invalid tensor domain");
        if (selection) {
            if (selection->coordinates.size() != extent || !b.same(selection->budget))
                throw std::invalid_argument("invalid tensor selection/ledger");
        } else if (extent) add(begin, mul(extent - 1, step));
    }
    std::size_t coordinate(std::size_t i) const {
        return selection ? selection->coordinates[i] : begin + i * step;
    }
    bool operator==(const Domain& d) const {
        if (space != d.space || parent != d.parent || parent_coordinate != d.parent_coordinate || extent != d.extent)
            return false;
        if (!selection && !d.selection) return !extent || (begin == d.begin && (extent == 1 || step == d.step));
        for (std::size_t i = 0; i < extent; ++i) if (coordinate(i) != d.coordinate(i)) return false;
        return true;
    }
};

struct Axis { Domain domain; std::ptrdiff_t stride = 0; }; // stride in elements
// Fixed-size schemas at registered kernel boundaries use the same rank/axis interface.
template<std::size_t N> struct Schema {
    std::array<Axis, N> axes;
    std::size_t rank() const { return N; }
    const Axis& axis(std::size_t i) const { return axes[i]; }
};
enum class Access { read, write, read_write };
enum class DType { fp64, complex_fp64 };
template<class T> constexpr DType dtype() {
    static_assert(std::is_same_v<T, double> || std::is_same_v<T, std::complex<double>>, "FP64 types only");
    return std::is_same_v<T, double> ? DType::fp64 : DType::complex_fp64;
}

class BorrowScope {
public:
    struct State {
        std::atomic<bool> active{true};
        std::atomic<std::size_t> claimed{0};
        std::size_t imported = 0;
    };
private:
    Budget budget_;
    Budget::Reservation control_;
    std::shared_ptr<State> state_;
public:
    explicit BorrowScope(const Budget& b, std::size_t imported_live_bytes = 0)
        : budget_(b), control_(b.reserve(add(sizeof(BorrowScope) + sizeof(State), imported_live_bytes))),
          state_(std::make_shared<State>()) { state_->imported = imported_live_bytes; }
    // Import an authoritative method inventory once. Live borrowed extents
    // consume credits within it; this is not another allocation or a cap waiver.
    void claim_imported(std::size_t bytes) {
        auto old = state_->claimed.load();
        do {
            if (bytes > state_->imported - old) throw std::length_error("tensor imported inventory exhausted");
        } while (!state_->claimed.compare_exchange_weak(old, old + bytes));
    }
    bool same_budget(const Budget& b) const { return budget_.same(b); }
    BorrowScope(const BorrowScope&) = delete;
    BorrowScope& operator=(const BorrowScope&) = delete;
    ~BorrowScope() { state_->active = false; }
    std::weak_ptr<State> token() const { return state_; }
};

struct Allocation {
    Budget budget;
    Budget::Reservation reservation;
    void* base;
    std::size_t bytes;
    DType type;
    bool writable;
    std::shared_ptr<void> owner;
    std::weak_ptr<BorrowScope::State> scope;
    std::size_t imported_bytes = 0;
    Space identity = Space::create();
    std::atomic<std::uint64_t> version{0};
    std::size_t readers = 0;
    bool writer = false;
    std::mutex mutex;
    Allocation(const Budget& b, void* p, std::size_t n, DType t, bool w,
               std::shared_ptr<void> owned, std::weak_ptr<BorrowScope::State> borrowed,
               Budget::Reservation r)
        : budget(b), reservation(std::move(r)), base(p), bytes(n), type(t), writable(w),
          owner(std::move(owned)), scope(std::move(borrowed)) {}
    ~Allocation() {
        if (imported_bytes) if (auto s = scope.lock()) s->claimed.fetch_sub(imported_bytes);
    }
    void check_lifetime() const {
        if (owner) return;
        auto s = scope.lock();
        if (!s || !s->active) throw std::runtime_error("expired tensor borrow");
    }
};

template<class T> std::shared_ptr<Allocation> own(const Budget& b, std::size_t n) {
    auto bytes = mul(n, sizeof(T)); signed_size(bytes);
    auto reservation = b.reserve(add(sizeof(Allocation), bytes));
    // Allocate even for an empty owner so it has a retainable lifetime token.
    std::shared_ptr<void> storage(new T[n], [](void* p) { delete[] static_cast<T*>(p); });
    return std::make_shared<Allocation>(b, storage.get(), bytes, dtype<T>(), true, storage,
                                       std::weak_ptr<BorrowScope::State>{}, std::move(reservation));
}
template<class T> std::shared_ptr<Allocation> borrow(const Budget& b, BorrowScope& scope,
                                                    T* base, std::size_t elements, bool imported = false) {
    using Value = std::remove_const_t<T>;
    if (!scope.same_budget(b)) throw std::invalid_argument("tensor borrow ledger mismatch");
    auto bytes = mul(elements, sizeof(Value)); signed_size(bytes);
    if (bytes && !base) throw std::invalid_argument("missing tensor buffer is not exact zero");
    add(reinterpret_cast<std::uintptr_t>(base), bytes);
    if (bytes && reinterpret_cast<std::uintptr_t>(base) % alignof(Value))
        throw std::invalid_argument("misaligned tensor allocation");
    auto reservation = b.reserve(add(sizeof(Allocation), imported ? 0 : bytes));
    auto allocation = std::make_shared<Allocation>(b, const_cast<Value*>(base), bytes, dtype<Value>(),
        !std::is_const_v<T>, std::shared_ptr<void>{}, scope.token(), std::move(reservation));
    if (imported) { scope.claim_imported(bytes); allocation->imported_bytes = bytes; }
    return allocation;
}

// Leases are synchronous. They retain owned allocations and reject conflicts;
// an external importer promises no mutation/free through untracked aliases.
class Lease {
    std::shared_ptr<Allocation> allocation_;
    bool write_ = false;
public:
    Lease(std::shared_ptr<Allocation> a, bool write) : allocation_(checked_handle(std::move(a))), write_(write) {
        std::lock_guard<std::mutex> lock(allocation_->mutex);
        allocation_->check_lifetime();
        if (allocation_->writer || (write && allocation_->readers) || (write && !allocation_->writable))
            throw std::runtime_error("conflicting tensor lease");
        if (write) {
            if (allocation_->version == UINT64_MAX) throw std::overflow_error("tensor value revision exhausted");
            allocation_->writer = true; ++allocation_->version;
        } else ++allocation_->readers;
    }
    Lease(const Lease&) = delete;
    Lease& operator=(const Lease&) = delete;
    ~Lease() {
        std::lock_guard<std::mutex> lock(allocation_->mutex);
        if (write_) allocation_->writer = false; else --allocation_->readers;
    }
};

class View {
    struct Data {
        Budget::Reservation control;
        std::shared_ptr<Allocation> allocation;
        Array<Axis> axes;
        std::ptrdiff_t offset, low = 0, high = 0;
        std::size_t elements = 1, width;
        Access access;
        bool conjugate;
        Data(std::shared_ptr<Allocation> a, const std::vector<Axis>& ax, std::ptrdiff_t off,
             Access mode, bool conj)
            : control(a->budget.reserve(add(sizeof(Data), mul(ax.capacity(), sizeof(Axis))))), allocation(std::move(a)),
              axes(allocation->budget, ax.size()), offset(off),
              width(allocation->type == DType::fp64 ? sizeof(double) : sizeof(std::complex<double>)),
              access(mode), conjugate(conj) {
            allocation->check_lifetime();
            if (mode != Access::read && (!allocation->writable || conj))
                throw std::invalid_argument("invalid tensor write access/conjugation");
            bool empty = false;
            for (const auto& x : ax) {
                x.domain.validate(allocation->budget);
                if (x.domain.selection && !allocation->budget.same(x.domain.selection->budget))
                    throw std::invalid_argument("tensor domain uses a different ledger");
                empty |= x.domain.extent == 0;
            }
            if (empty) elements = 0;
            for (std::size_t i = 0; i < ax.size(); ++i) {
                axes[i] = ax[i];
                if (!empty) elements = mul(elements, ax[i].domain.extent);
            }
            if (off < 0) throw std::invalid_argument("negative tensor origin");
            if (empty) {
                if (std::size_t(off) > allocation->bytes / width) throw std::out_of_range("empty tensor origin");
                return; // No pointer formation, even for a null base.
            }
            low = high = off;
            for (const auto& x : ax) {
                auto span = stride_span(x.domain.extent - 1, x.stride);
                if (span < 0) low = signed_add(low, span); else high = signed_add(high, span);
            }
            if (low < 0 || mul(add(std::size_t(high), 1), width) > allocation->bytes)
                throw std::out_of_range("tensor view exceeds allocation");
            if (!allocation->base || reinterpret_cast<std::uintptr_t>(allocation->base) % alignof(double))
                throw std::invalid_argument("invalid tensor buffer alignment");
            // A conservative injectivity proof: ascending absolute strides
            // must each span all faster dimensions. No sorting allocation.
            if (mode != Access::read) {
                std::size_t covered = 1, previous = 0, remaining = 0;
                for (const auto& x : ax) if (x.domain.extent > 1) ++remaining;
                while (remaining--) {
                    std::size_t smallest = SIZE_MAX, selected = ax.size();
                    for (std::size_t i = 0; i < ax.size(); ++i) if (ax[i].domain.extent > 1) {
                        auto s = ax[i].stride;
                        auto magnitude = s < 0 ? std::size_t(-(s + 1)) + 1 : std::size_t(s);
                        if (magnitude > previous && magnitude < smallest) { smallest = magnitude; selected = i; }
                    }
                    if (selected == ax.size() || smallest < covered)
                        throw std::invalid_argument("overlapping tensor output");
                    covered = add(covered, mul(ax[selected].domain.extent - 1, smallest));
                    previous = smallest;
                }
            }
        }
    };
    std::shared_ptr<const Data> data_;
public:
    View(std::shared_ptr<Allocation> allocation, const std::vector<Axis>& axes,
         std::ptrdiff_t offset = 0, Access access = Access::read, bool conjugate = false)
        : data_(std::make_shared<Data>(checked_handle(std::move(allocation)), axes, offset, access, conjugate)) {}
    std::size_t rank() const { return data_->axes.size(); }
    std::size_t size() const { return data_->elements; }
    const Axis& axis(std::size_t i) const { return data_->axes[i]; }
    std::ptrdiff_t offset() const { return data_->offset; }
    Access access() const { return data_->access; }
    bool conjugate() const { return data_->conjugate; }
    const std::shared_ptr<Allocation>& allocation() const { return data_->allocation; }
    bool overlaps(const View& v) const {
        if (!size() || !v.size()) return false;
        const auto bounds = [](const Data& d) {
            auto base = reinterpret_cast<std::uintptr_t>(d.allocation->base);
            return std::make_pair(add(base, mul(std::size_t(d.low), d.width)),
                                 add(base, mul(add(std::size_t(d.high), 1), d.width)));
        };
        auto a = bounds(*data_), b = bounds(*v.data_);
        return a.first < b.second && b.first < a.second;
    }
};

// Binary indexed contraction/reduction, including repeated input labels
// (diagonals), batches and scalars. No implicit broadcasting or conjugation.
// No algebraic reassociation/CSE: TensorRight doi:10.1145/3704865 does not
// establish IEEE floating-point or layout-sensitive rewrite equivalence.
// The same checked index program describes affine operations and provider
// schedules. Provider axes have zero address strides: no affine allocation is
// inferred from a logical provider schema.
class IndexProgram {
public:
    struct Index { std::size_t label = 0; Domain domain; std::ptrdiff_t a = 0, b = 0, c = 0; bool output = false; };
private:
    Budget::Reservation control_;
    Array<Index> indices_;
    std::size_t count_ = 0, outputs_ = 1, reductions_ = 1;
public:
    template<class A, class B, class C>
    IndexProgram(const Budget& budget, const A& a, const std::vector<std::size_t>& al,
                 const B& b, const std::vector<std::size_t>& bl,
                 const C& c, const std::vector<std::size_t>& cl)
        : control_(budget.reserve(add(sizeof(IndexProgram), mul(add(add(al.capacity(), bl.capacity()), cl.capacity()), sizeof(std::size_t))))), indices_(budget, add(al.size(), bl.size())) {
        if (al.size() != a.rank() || bl.size() != b.rank() || cl.size() != c.rank())
            throw std::invalid_argument("tensor label rank mismatch");
        auto insert = [&](const auto& v, const std::vector<std::size_t>& labels, int operand) {
            for (std::size_t k = 0; k < labels.size(); ++k) {
                v.axis(k).domain.validate(budget);
                std::size_t i = 0;
                while (i < count_ && indices_[i].label != labels[k]) ++i;
                if (i == count_) {
                    if (operand == 2) throw std::invalid_argument("unknown tensor output label");
                    indices_[i].label = labels[k]; indices_[i].domain = v.axis(k).domain; ++count_;
                } else if (!(indices_[i].domain == v.axis(k).domain)) {
                    throw std::invalid_argument("tensor index domain mismatch");
                }
                auto& x = indices_[i];
                if (operand == 0) x.a = signed_add(x.a, v.axis(k).stride);
                if (operand == 1) x.b = signed_add(x.b, v.axis(k).stride);
                if (operand == 2) {
                    if (x.output) throw std::invalid_argument("duplicate tensor output label");
                    x.output = true; x.c = v.axis(k).stride;
                }
            }
        };
        insert(a, al, 0); insert(b, bl, 1); insert(c, cl, 2);
        bool empty_out = false, empty_red = false;
        for (std::size_t i = 0; i < count_; ++i)
            if (!indices_[i].domain.extent) (indices_[i].output ? empty_out : empty_red) = true;
        outputs_ = empty_out ? 0 : 1; reductions_ = empty_red ? 0 : 1;
        for (std::size_t i = 0; i < count_; ++i) {
            const auto& x = indices_[i];
            if (x.output && !empty_out) outputs_ = mul(outputs_, x.domain.extent);
            if (!x.output && !empty_red) reductions_ = mul(reductions_, x.domain.extent);
        }
        mul(outputs_, reductions_);
    }
    std::size_t count() const { return count_; }
    const Index& index(std::size_t i) const { return indices_[i]; }
    std::size_t output_elements() const { return outputs_; }
    std::size_t reduction_elements() const { return reductions_; }
};

class Contraction {
    Budget::Reservation control_;
    View a_, b_, c_;
    IndexProgram program_;
    std::uint64_t av_, bv_;
public:
    Contraction(View a, const std::vector<std::size_t>& al, View b, const std::vector<std::size_t>& bl,
                View c, const std::vector<std::size_t>& cl)
        : control_(a.allocation()->budget.reserve(sizeof(Contraction))), a_(std::move(a)), b_(std::move(b)), c_(std::move(c)),
          program_(a_.allocation()->budget, a_, al, b_, bl, c_, cl),
          av_(a_.allocation()->version), bv_(b_.allocation()->version) {
        if (!a_.allocation()->budget.same(b_.allocation()->budget) || !a_.allocation()->budget.same(c_.allocation()->budget))
            throw std::invalid_argument("tensor operation spans different ledgers");
        if (a_.access() == Access::write || b_.access() == Access::write || c_.access() == Access::read)
            throw std::invalid_argument("tensor IR access mismatch");
        if (a_.allocation()->type != b_.allocation()->type || a_.allocation()->type != c_.allocation()->type)
            throw std::invalid_argument("tensor IR dtype mismatch");
        if (a_.overlaps(c_) || b_.overlaps(c_) || a_.allocation() == c_.allocation() || b_.allocation() == c_.allocation())
            throw std::invalid_argument("tensor input/output alias requires an explicit temporary");
    }
    std::size_t output_elements() const { return program_.output_elements(); }
    std::size_t reduction_elements() const { return program_.reduction_elements(); }
    const char* kernel() const { return "cpu.indexed.binary.v1"; }
    const IndexProgram& index_program() const { return program_; }
    template<class T> void execute(T alpha = T(1), T beta = T(0)) const {
        if (dtype<T>() != a_.allocation()->type) throw std::invalid_argument("tensor execution dtype mismatch");
        if (beta != T(0) && c_.access() != Access::read_write)
            throw std::invalid_argument("beta requires read-write tensor output");
        Lease la(a_.allocation(), false), lb(b_.allocation(), false);
        if (a_.allocation()->version != av_ || b_.allocation()->version != bv_)
            throw std::runtime_error("stale tensor value binding");
        Lease lc(c_.allocation(), true);
        const auto* ap = static_cast<const T*>(a_.allocation()->base);
        const auto* bp = static_cast<const T*>(b_.allocation()->base);
        auto* cp = static_cast<T*>(c_.allocation()->base);
        auto conjugated = [](T value, bool conjugate) -> T {
            if constexpr (std::is_same_v<T, std::complex<double>>) return conjugate ? std::conj(value) : value;
            else return value;
        };
        for (std::size_t out = 0; out < program_.output_elements(); ++out) {
            auto ao = a_.offset(), bo = b_.offset(), co = c_.offset();
            auto position = out;
            for (std::size_t i = 0; i < program_.count(); ++i) if (program_.index(i).output) {
                const auto& x = program_.index(i); auto k = position % x.domain.extent; position /= x.domain.extent;
                ao = signed_add(ao, stride_span(k, x.a)); bo = signed_add(bo, stride_span(k, x.b));
                co = signed_add(co, stride_span(k, x.c));
            }
            T sum = 0;
            for (std::size_t red = 0; alpha != T(0) && red < program_.reduction_elements(); ++red) {
                auto ai = ao, bi = bo; position = red;
                for (std::size_t i = 0; i < program_.count(); ++i) if (!program_.index(i).output) {
                    const auto& x = program_.index(i); auto k = position % x.domain.extent; position /= x.domain.extent;
                    ai = signed_add(ai, stride_span(k, x.a)); bi = signed_add(bi, stride_span(k, x.b));
                }
                sum += conjugated(ap[ai], a_.conjugate()) * conjugated(bp[bi], b_.conjugate());
            }
            // beta=0 must not read an uninitialized or NaN destination.
            cp[co] = beta == T(0) ? alpha * sum : alpha * sum + beta * cp[co];
        }
    }
};

} // namespace vibeqc::tensor