Fix #355: bound parser recursion depth — nested input raises parse_error

.github/workflows/clang-tidy-check.yml

@@ -30,4 +30,4 @@ jobs:
       - name: Run clang-tidy
         run: |
           find src -name '*.cpp' -print0 \
-          | xargs -0 clang-tidy-18 -p build
+          | xargs -0 clang-tidy-18 -p build --warnings-as-errors='*'

CMakeLists.txt

@@ -194,9 +194,10 @@ endif()
 # Sanitizers
 if(NUMSIM_CAS_SANITIZERS)
     target_compile_options(${PROJECT_NAME} PUBLIC
-        -fsanitize=address,undefined -fno-omit-frame-pointer)
+        -fsanitize=address,undefined -fno-sanitize-recover=undefined
+        -fno-omit-frame-pointer)
     target_link_options(${PROJECT_NAME} PUBLIC
-        -fsanitize=address,undefined)
+        -fsanitize=address,undefined -fno-sanitize-recover=undefined)
 endif()
 
 # Optional convenience for Windows
@@ -274,9 +275,10 @@ if(NUMSIM_CAS_BUILD_PARSER)
 
   if(NUMSIM_CAS_SANITIZERS)
     target_compile_options(NumSim_CAS_Parser PUBLIC
-      -fsanitize=address,undefined -fno-omit-frame-pointer)
+      -fsanitize=address,undefined -fno-sanitize-recover=undefined
+      -fno-omit-frame-pointer)
     target_link_options(NumSim_CAS_Parser PUBLIC
-      -fsanitize=address,undefined)
+      -fsanitize=address,undefined -fno-sanitize-recover=undefined)
   endif()
 
   if(WIN32)

include/numsim_cas/core/hash_functions.h

@@ -1,7 +1,9 @@
 #ifndef HASH_FUNCTIONS_H
 #define HASH_FUNCTIONS_H
 
+#include <bit>
 #include <complex>
+#include <cstdint>
 #include <string>
 #include <vector>
 
@@ -15,6 +17,18 @@ inline void hash_combine(std::size_t &seed, const T &value) {
           static_cast<std::size_t>(0x9e3779b9) + (seed << 6) + (seed >> 2);
 }
 
+// Doubles hash by bit pattern: the generic static_cast is UB for negative
+// values and collapses all fractions in (0,1) onto 0 (#361).
+inline void hash_combine(std::size_t &seed, double value) {
+  if (value == 0.0)
+    value = 0.0; // normalize -0.0
+  hash_combine(seed, std::bit_cast<std::uint64_t>(value));
+}
+
+inline void hash_combine(std::size_t &seed, float value) {
+  hash_combine(seed, static_cast<double>(value));
+}
+
 inline void hash_combine(std::size_t &seed, const std::string &value) {
   for (const auto &c : value) {
     hash_combine(seed, c);

include/numsim_cas/core/scalar_number.h

@@ -96,7 +96,7 @@ class scalar_number {
                                           scalar_number const &exp);
 
 private:
-  explicit scalar_number(variant_t vv) : v_(std::move(vv)) {
+  explicit scalar_number(variant_t vv) : v_(vv) {
     if (auto *r = std::get_if<rational_t>(&v_)) {
       if (r->den == 1)
         v_ = r->num;
@@ -106,6 +106,27 @@ class scalar_number {
   variant_t v_;
 };
 
+// Value-normalizing hash: numerically equal constants must hash equal
+// across variant alternatives (int64 2 vs double 2.0 vs rational 2/1),
+// while doubles otherwise hash by bit pattern (#361).
+inline void hash_combine(std::size_t &seed, scalar_number const &value) {
+  std::visit(
+      [&](auto const &x) {
+        using T = std::decay_t<decltype(x)>;
+        if constexpr (std::is_same_v<T, double>) {
+          // guard BEFORE casting: the int64 cast is UB for NaN/inf and
+          // |x| >= 2^63 (review on #361)
+          if (x >= -9.2e18 && x <= 9.2e18 &&
+              x == static_cast<double>(static_cast<std::int64_t>(x))) {
+            hash_combine(seed, static_cast<std::int64_t>(x));
+            return;
+          }
+        }
+        hash_combine(seed, x);
+      },
+      value.raw());
+}
+
 } // namespace numsim::cas
 
 #endif // SCALAR_NUMBER_H

include/numsim_cas/parser/parse_error.h

@@ -37,7 +37,7 @@ class parse_error : public cas_error {
 public:
   /// Construct with optional source context. Pass an empty `source` and
   /// `byte_offset = 0` for errors raised outside the parser.
-  parse_error(std::string message, std::size_t byte_offset,
+  parse_error(std::string const &message, std::size_t byte_offset,
               std::string_view source);
 
   /// Byte offset into the source where the error was detected, clamped

include/numsim_cas/scalar/scalar_constant.h

@@ -61,8 +61,7 @@ class scalar_constant final : public scalar_node_base_t<scalar_constant> {
   void update_hash_value() const noexcept override {
     this->m_hash_value = 0;
     hash_combine(this->m_hash_value, this->id());
-    std::visit([&](auto const &x) { hash_combine(this->m_hash_value, x); },
-               m_value.raw());
+    hash_combine(this->m_hash_value, m_value);
   }
 
 private:

include/numsim_cas/scalar/visitors/scalar_rebuild_visitor.h

@@ -14,7 +14,7 @@ class scalar_rebuild_visitor : public scalar_visitor_const_t {
 public:
   using expr_holder_t = expression_holder<scalar_expression>;
 
-  virtual ~scalar_rebuild_visitor() = default;
+  ~scalar_rebuild_visitor() override = default;
 
   virtual expr_holder_t apply(expr_holder_t const &expr) {
     if (expr.is_valid()) {

include/numsim_cas/tensor/data/tensor_data_to_scalar_wrapper.h

@@ -4,6 +4,7 @@
 #include "spectral_decomposition_cache.h"
 #include "tensor_data.h"
 #include <numsim_cas/core/cas_error.h>
+#include <numsim_cas/tensor/sequence.h>
 
 #include <algorithm>
 #include <array>
@@ -58,18 +59,47 @@ class tensor_data_dcontract_wrapper final
     : public tensor_data_eval_up_unary<tensor_data_dcontract_wrapper<ValueType>,
                                        ValueType> {
 public:
+  // #353 — the contraction sequences are part of the node semantics:
+  // {1,2}/{2,1} is A_ij*B_ji (A : B^T), not A : B.
   tensor_data_dcontract_wrapper(tensor_data_base<ValueType> const &lhs,
-                                tensor_data_base<ValueType> const &rhs)
-      : m_lhs(lhs), m_rhs(rhs) {}
+                                tensor_data_base<ValueType> const &rhs,
+                                sequence const &lhs_indices,
+                                sequence const &rhs_indices)
+      : m_lhs(lhs), m_rhs(rhs), m_lhs_indices(lhs_indices),
+        m_rhs_indices(rhs_indices) {}
 
   template <std::size_t Dim, std::size_t Rank> ValueType evaluate_imp() {
+    // The dispatch picks Dim/Rank from the LHS; a mixed-rank/dim node
+    // (constructible through the weak dot_product precondition, #360)
+    // would type-pun the RHS cast into silent garbage (review on #353).
+    if (m_lhs.rank() != m_rhs.rank() || m_lhs.dim() != m_rhs.dim()) {
+      throw evaluation_error(
+          "tensor_data_dcontract_wrapper: operand rank/dim mismatch");
+    }
+    if (m_lhs_indices.size() != Rank || m_rhs_indices.size() != Rank) {
+      throw evaluation_error(
+          "tensor_data_dcontract_wrapper: sequence size != operand rank");
+    }
     if constexpr (Rank == 2) {
       using Tensor = tensor_data<ValueType, Dim, Rank>;
       auto const &l = static_cast<const Tensor &>(m_lhs).data();
       auto const &r = static_cast<const Tensor &>(m_rhs).data();
-      return static_cast<ValueType>(tmech::dcontract(l, r));
+      const bool straight{m_lhs_indices == m_rhs_indices};
+      if (straight) {
+        // {1,2}/{1,2} = A_ij B_ij; {2,1}/{2,1} = A_ji B_ji = same sum
+        return static_cast<ValueType>(tmech::dcontract(l, r));
+      }
+      // {1,2}/{2,1} (either orientation) = A_ij B_ji
+      return static_cast<ValueType>(tmech::dcontract(l, tmech::trans(r)));
+    } else if constexpr (Rank == 1) {
+      using Tensor = tensor_data<ValueType, Dim, Rank>;
+      auto const &l = static_cast<const Tensor &>(m_lhs).data();
+      auto const &r = static_cast<const Tensor &>(m_rhs).data();
+      return static_cast<ValueType>(tmech::dot(l, r));
     } else {
-      throw evaluation_error("tensor_data_dcontract_wrapper: requires rank 2");
+      throw evaluation_error(
+          "tensor_data_dcontract_wrapper: rank > 2 contraction not "
+          "implemented (#353 follow-up in #383)");
     }
   }
 
@@ -86,6 +116,8 @@ class tensor_data_dcontract_wrapper final
 private:
   tensor_data_base<ValueType> const &m_lhs;
   tensor_data_base<ValueType> const &m_rhs;
+  sequence m_lhs_indices;
+  sequence m_rhs_indices;
 };
 
 // ─── Eigenvalue wrapper: dispatches runtime (dim,rank), computes the

include/numsim_cas/tensor/identity_tensor.h

@@ -106,10 +106,9 @@ class identity_tensor final : public tensor_node_base_t<identity_tensor> {
       a.insert(positive_semidefinite{}); // PD ⇒ PSD, mirrors
                                          // assume_positive_definite()
       // PD ⇒ symmetric, via the same helper assume_positive_definite()
-      // uses. Mechanically a no-op today because space_for_rank above
-      // already wrote a qualifying space tag at both supported ranks
-      // (rank-2 Symmetric ⇒ ProjKind::Sym early-return; rank-4
-      // MinorMajor ⇒ holds_alternative early-return). The call exists
+      // uses. Mechanically a no-op today: rank-2 already carries
+      // Symmetric, and the helper's rank()!=2 gate skips rank 4 (whose
+      // tag is Major since #351). The call exists
       // so that any future drift in detail::set_symmetric_unless_more_specific
       // (e.g. it grows a sibling-field write that we'd otherwise miss)
       // tracks automatically. Pairs with the helper in tensor_assume.h.
@@ -154,16 +153,16 @@ class identity_tensor final : public tensor_node_base_t<identity_tensor> {
 
 private:
   // Closed-form structural classification by rank. Rank-2 is Kronecker δ_ij
-  // (Symmetric). Rank-4 minor identity is δ_ik·δ_jl (MinorMajor — fully
-  // symmetric at rank-4). Higher ranks return nullopt; the variant has no
-  // general "all-pairs-minor" alternative and higher-rank identity is
-  // rarely queried. Open decision tracked in
-  // docs/sympy-assumption-redesign.md.
+  // (Symmetric). Rank-4 identity δ_ik·δ_jl has MAJOR symmetry only:
+  // swapping (ij)<->(kl) gives δ_ki·δ_lj = δ_ik·δ_jl, but swapping i<->j
+  // gives δ_jk·δ_il ≠ δ_ik·δ_jl (the minor-symmetric identity is P_sym,
+  // a different tensor). The former MinorMajor tag routed inv() through
+  // the symmetric Voigt path, evaluating inv(-I4) wrongly (#351).
   static std::optional<tensor_space> space_for_rank(std::size_t rank) noexcept {
     if (rank == 2)
       return tensor_space{Symmetric{}, AnyTraceTag{}};
     if (rank == 4)
-      return tensor_space{MinorMajor{}, AnyTraceTag{}};
+      return tensor_space{Major{}, AnyTraceTag{}};
     return std::nullopt;
   }
 };

include/numsim_cas/tensor/tensor_std.h

@@ -39,6 +39,27 @@ namespace numsim::cas {
 
 template <tensor_expr_holder ExprLHS, scalar_expr_holder ExprRHS>
 [[nodiscard]] auto pow(ExprLHS &&expr_lhs, ExprRHS &&expr_rhs) {
+  // Matrix powers are rank-2 only; the evaluator's repeated contraction
+  // and the diff kernels are meaningless for other ranks (#350).
+  if (expr_lhs.get().rank() != 2) {
+    throw invalid_expression_error(
+        "pow: tensor operand must be rank 2 (got rank " +
+        std::to_string(expr_lhs.get().rank()) + ")");
+  }
+  // Non-integer constant exponents have no matrix-power meaning here;
+  // isotropic tensor functions (#227) cover fractional powers (#350).
+  {
+    // strip any depth of negation before the literal check, matching
+    // try_int_constant's own recursion (round-2 review on #350)
+    expression_holder<scalar_expression> probe{expr_rhs};
+    while (is_same<scalar_negative>(probe)) {
+      probe = probe.template get<scalar_negative>().expr();
+    }
+    if (is_same<scalar_constant>(probe) && !try_int_constant(expr_rhs)) {
+      throw invalid_expression_error(
+          "pow: tensor exponent must be an integer constant");
+    }
+  }
   // pow(0, n) → 0
   if (is_same<tensor_zero>(expr_lhs))
     return make_expression<tensor_zero>(expr_lhs.get().dim(),

include/numsim_cas/tensor/visitors/tensor_evaluator.h

@@ -1,6 +1,8 @@
 #ifndef TENSOR_EVALUATOR_H
 #define TENSOR_EVALUATOR_H
 
+#include <cmath>
+
 #include <algorithm>
 #include <cstring>
 #include <map>
@@ -255,8 +257,20 @@ class tensor_evaluator final : public tensor_visitor_const_t {
   // ─── Tensor functions (tmech wrappers) ─────────────────────
 
   void operator()(tensor_pow const &visitable) override {
+    if (visitable.rank() != 2) {
+      // the repeated single-index contraction is rank-2-only; a rebuilt
+      // tree can bypass the factory gate (review on #350: substitution
+      // recreated a rank-4 pow and corrupted the heap)
+      throw evaluation_error("tensor_pow: rank-2 operand required");
+    }
     auto base_data = apply(visitable.expr_lhs());
     const auto exp_val = m_scalar_eval.apply(visitable.expr_rhs());
+    if (exp_val != std::round(exp_val) || exp_val < -1e6 || exp_val > 1e6) {
+      // silently truncating pow(A, 0.5) to the identity was #350; the
+      // range check keeps the int cast below defined for huge values
+      throw evaluation_error(
+          "tensor_pow: exponent must evaluate to a moderate integer");
+    }
     const auto n = static_cast<int>(exp_val);
     const auto d = visitable.dim();
     const auto r = visitable.rank();
@@ -284,6 +298,14 @@ class tensor_evaluator final : public tensor_visitor_const_t {
       ip.evaluate(d, r, r);
       accumulated = std::move(temp);
     }
+    if (n < 0) {
+      // A^-n = inv(A^n); without this, pow(X,-1) returned X itself (#350)
+      auto inverted = make_tensor_data<ValueType>(d, r);
+      tensor_data_unary_wrapper<tmech_ops::inv, ValueType> iv(*inverted,
+                                                              *accumulated);
+      iv.evaluate(d, r);
+      accumulated = std::move(inverted);
+    }
     m_result = std::move(accumulated);
   }
 
@@ -385,6 +407,14 @@ class tensor_evaluator final : public tensor_visitor_const_t {
   template <typename Op>
   void eval_projector_unary(inner_product_wrapper const &visitable) {
     auto rhs_data = apply(visitable.expr_rhs());
+    // round-3 review: the wrapper's dim() is the projector's; a
+    // dim-changing substitution produced an operand whose buffer this
+    // shortcut then over-read (ASan heap-buffer-overflow)
+    if (rhs_data->dim() != visitable.dim() ||
+        rhs_data->rank() != visitable.rank()) {
+      throw evaluation_error(
+          "projector contraction: operand dim/rank mismatch");
+    }
     m_result = make_tensor_data<ValueType>(visitable.dim(), visitable.rank());
     tensor_data_unary_wrapper<Op, ValueType> op(*m_result, *rhs_data);
     op.evaluate(visitable.dim(), visitable.rank());

include/numsim_cas/tensor/visitors/tensor_rebuild_visitor.h

@@ -4,6 +4,7 @@
 #include <numsim_cas/basic_functions.h>
 #include <numsim_cas/core/operators.h>
 #include <numsim_cas/eigen_decomposition.h>
+#include <numsim_cas/tensor/projector_algebra.h>
 #include <numsim_cas/tensor/tensor_definitions.h>
 #include <numsim_cas/tensor/tensor_functions.h>
 #include <numsim_cas/tensor/tensor_operators.h>
@@ -25,10 +26,14 @@ class tensor_rebuild_visitor : public tensor_visitor_const_t {
       expr.get<tensor_visitable_t>().accept(*this);
       // #93 — reconstructions below build fresh nodes (variadic ctors)
       // that drop post-construction space(). Restore from source, but not
-      // over a self-computed space (e.g. tensor_add's child-join).
+      // over a self-computed space (e.g. tensor_add's child-join), and only
+      // for structurally unchanged rebuilds: a subclass that swapped
+      // children (substitution) must not inherit the source's space (#352).
       if (m_result.is_valid() && !m_result.get().space()) {
-        if (auto const &sp = expr.get().space())
-          m_result.data()->set_space(*sp);
+        if (auto const &sp = expr.get().space()) {
+          if (m_result == expr || same_projector_contraction(expr, m_result))
+            m_result.data()->set_space(*sp);
+        }
       }
       return std::move(m_result);
     }
@@ -41,6 +46,22 @@ class tensor_rebuild_visitor : public tensor_visitor_const_t {
 
   virtual t2s_holder_t apply_t2s(t2s_holder_t const &expr) { return expr; }
 
+  // skew(X)/sym(X)/...: the space comes from the projector, not the
+  // argument, so it survives child substitution (#93/#352).
+  static bool same_projector_contraction(tensor_holder_t const &a,
+                                         tensor_holder_t const &b) {
+    auto ia = as_projector_contraction(a);
+    auto ib = as_projector_contraction(b);
+    if (!(ia && ib && *ia->proj == *ib->proj)) {
+      return false;
+    }
+    // round-2 review: the wrapper's dim() comes from the projector LHS, so
+    // node-level shape checks are inert - compare the actual arguments
+    // (a dim-changing substitution reached a heap overflow at evaluation)
+    return ia->argument.get().rank() == ib->argument.get().rank() &&
+           ia->argument.get().dim() == ib->argument.get().dim();
+  }
+
   // Leaf nodes: return as-is
   void operator()(tensor const &) override { m_result = m_current; }
   void operator()(tensor_zero const &) override { m_result = m_current; }

include/numsim_cas/tensor/wrappers/tensor_inv.h

@@ -12,8 +12,8 @@ class tensor_inv final : public unary_op<tensor_node_base_t<tensor_inv>> {
   using base = unary_op<tensor_node_base_t<tensor_inv>>;
 
   template <typename Expr>
-  explicit tensor_inv(
-      Expr &&_expr) // NOLINT(bugprone-forwarding-reference-overload)
+  // NOLINTNEXTLINE(bugprone-forwarding-reference-overload)
+  explicit tensor_inv(Expr &&_expr)
       : base(std::forward<Expr>(_expr), _expr.get().dim(), _expr.get().rank()) {
     // ── Rank gate (#292) ──────────────────────────────────────────
     // Mirror the inv() factory's rank gate (tensor_functions.h:467)

include/numsim_cas/tensor_to_scalar/visitors/tensor_to_scalar_evaluator.h

@@ -160,7 +160,8 @@ class tensor_to_scalar_evaluator final
     auto rhs_data = m_tensor_eval.apply(v.expr_rhs());
     const auto dim = lhs_data->dim();
     const auto rank = lhs_data->rank();
-    tensor_data_dcontract_wrapper<ValueType> op(*lhs_data, *rhs_data);
+    tensor_data_dcontract_wrapper<ValueType> op(
+        *lhs_data, *rhs_data, v.indices_lhs(), v.indices_rhs());
     m_result = op.evaluate(dim, rank);
   }
 

include/numsim_cas/tensor_to_scalar/visitors/tensor_to_scalar_rebuild_visitor.h

@@ -19,7 +19,7 @@ class tensor_to_scalar_rebuild_visitor
   using scalar_holder_t = expression_holder<scalar_expression>;
   using tensor_holder_t = expression_holder<tensor_expression>;
 
-  virtual ~tensor_to_scalar_rebuild_visitor() = default;
+  ~tensor_to_scalar_rebuild_visitor() override = default;
 
   virtual t2s_holder_t apply(t2s_holder_t const &expr) {
     if (expr.is_valid()) {