From 761bf647ba605b72ac31cbcced834fdd88fd991a Mon Sep 17 00:00:00 2001 From: petlenz Date: Tue, 15 Sep 2026 23:29:37 +0200 Subject: [PATCH 1/5] Fix #357: limits report a sign only when it is provable Limit algebra: - finite + finite is positive or negative only when both signs agree; mixed signs give unknown - log(c) for finite positive c is unknown (negative below 1) - 1/0, 0^(-c) and log(0) need the zero to be approached from above; callers pass that knowledge, otherwise the result is unknown - c^0 with c -> 0 or c -> inf is indeterminate Scalar limit visitor: - a symbol other than the limit variable takes its sign from its assumptions, otherwise unknown - sin, tan, asin, atan of an argument tending to 0 tend to 0; cos and acos tend to a positive value; other finite arguments of sin, cos, tan and a positive argument of acos give unknown - asin of an infinite argument is unknown; atan and asin keep the argument's sign - sign(x) at 0 follows the side of approach Tensor-to-scalar limit visitor: - invariants independent of the limit variable take their sign from their assumptions, otherwise unknown - a scalar wrapper takes the sign of its scalar value - pow and log of the exact limit variable use the side of approach Signed-off-by: petlenz --- include/numsim_cas/core/limit_algebra.h | 9 +- .../scalar/visitors/scalar_limit_visitor.h | 3 + .../visitors/tensor_to_scalar_limit_visitor.h | 1 + src/numsim_cas/core/limit_algebra.cpp | 41 +++-- .../scalar/visitors/scalar_limit_visitor.cpp | 136 +++++++++++---- .../tensor_to_scalar_limit_visitor.cpp | 74 +++++--- tests/LimitVisitorTest.h | 161 ++++++++++++++++++ 7 files changed, 348 insertions(+), 77 deletions(-) diff --git a/include/numsim_cas/core/limit_algebra.h b/include/numsim_cas/core/limit_algebra.h index cde32413..393fa3c6 100644 --- a/include/numsim_cas/core/limit_algebra.h +++ b/include/numsim_cas/core/limit_algebra.h @@ -10,11 +10,14 @@ class limit_algebra { static limit_result combine_add(limit_result a, limit_result b); static limit_result combine_mul(limit_result a, limit_result b); static limit_result apply_neg(limit_result a); - static limit_result apply_log(limit_result a); - static limit_result apply_pow(limit_result base, limit_result exponent); + // zero_from_above: the argument that tends to zero is known to stay + // positive; without it 1/0, 0^(-c) and log(0) have no provable sign. + static limit_result apply_log(limit_result a, bool zero_from_above); + static limit_result apply_pow(limit_result base, limit_result exponent, + bool zero_from_above); static limit_result apply_sqrt(limit_result a); static limit_result apply_abs(limit_result a); - static limit_result apply_reciprocal(limit_result a); + static limit_result apply_reciprocal(limit_result a, bool zero_from_above); static limit_result apply_exp(limit_result a); }; diff --git a/include/numsim_cas/scalar/visitors/scalar_limit_visitor.h b/include/numsim_cas/scalar/visitors/scalar_limit_visitor.h index 67ed30f9..56e55cdc 100644 --- a/include/numsim_cas/scalar/visitors/scalar_limit_visitor.h +++ b/include/numsim_cas/scalar/visitors/scalar_limit_visitor.h @@ -57,6 +57,9 @@ class scalar_limit_visitor final : public scalar_visitor_const_t, void operator()(scalar_if_then_else const &) override; private: + bool zero_from_above(expr_holder_t const &expr) const; + bool zero_from_below(expr_holder_t const &expr) const; + expr_holder_t m_limit_var; limit_target m_target; limit_result m_result; diff --git a/include/numsim_cas/tensor_to_scalar/visitors/tensor_to_scalar_limit_visitor.h b/include/numsim_cas/tensor_to_scalar/visitors/tensor_to_scalar_limit_visitor.h index 4d88feef..1fb3a08a 100644 --- a/include/numsim_cas/tensor_to_scalar/visitors/tensor_to_scalar_limit_visitor.h +++ b/include/numsim_cas/tensor_to_scalar/visitors/tensor_to_scalar_limit_visitor.h @@ -56,6 +56,7 @@ class tensor_to_scalar_limit_visitor final private: bool depends_on_limit_var(t2s_holder_t const &expr) const; + bool zero_from_above(t2s_holder_t const &expr) const; dependency_mode m_mode; t2s_holder_t m_limit_var_t2s; diff --git a/src/numsim_cas/core/limit_algebra.cpp b/src/numsim_cas/core/limit_algebra.cpp index 60f86d65..213aea25 100644 --- a/src/numsim_cas/core/limit_algebra.cpp +++ b/src/numsim_cas/core/limit_algebra.cpp @@ -110,8 +110,12 @@ limit_result limit_algebra::combine_add(limit_result a, limit_result b) { if (is_finite(a.dir) && is_infinite(b.dir)) return b; - // finite + finite: finite (can't determine exact sign in general) - return {dir::finite_positive}; + // finite + finite: the sign is known only when both signs agree + if (a.dir == dir::finite_positive && b.dir == dir::finite_positive) + return {dir::finite_positive}; + if (a.dir == dir::finite_negative && b.dir == dir::finite_negative) + return {dir::finite_negative}; + return {dir::unknown}; } limit_result limit_algebra::combine_mul(limit_result a, limit_result b) { @@ -173,17 +177,19 @@ limit_result limit_algebra::apply_neg(limit_result a) { return {flip_sign(a.dir), a.rate}; } -limit_result limit_algebra::apply_log(limit_result a) { +limit_result limit_algebra::apply_log(limit_result a, bool zero_from_above) { if (a.dir == dir::indeterminate || a.dir == dir::unknown) return a; switch (a.dir) { case dir::zero: - // log(0+) = -inf (logarithmic) + // log(0+) = -inf (logarithmic); log is undefined when approached from below + if (!zero_from_above) + return {dir::unknown}; return {dir::neg_infinity, {gtype::logarithmic, 1.0}}; case dir::finite_positive: - // log(finite_positive) = finite - return {dir::finite_positive}; + // log(c) is negative for c < 1 and zero at c = 1 + return {dir::unknown}; case dir::finite_negative: // log(negative) is undefined in reals return {dir::unknown}; @@ -198,8 +204,8 @@ limit_result limit_algebra::apply_log(limit_result a) { } } -limit_result limit_algebra::apply_pow(limit_result base, - limit_result exponent) { +limit_result limit_algebra::apply_pow(limit_result base, limit_result exponent, + bool zero_from_above) { if (base.dir == dir::indeterminate || exponent.dir == dir::indeterminate) return {dir::indeterminate}; if (base.dir == dir::unknown || exponent.dir == dir::unknown) @@ -208,8 +214,12 @@ limit_result limit_algebra::apply_pow(limit_result base, // Only handle constant/finite exponents for now if (is_finite(exponent.dir) || exponent.dir == dir::zero) { if (exponent.dir == dir::zero) { - // x^0 = 1 - return {dir::finite_positive}; + // c^0 = 1, but 0^0 and inf^0 are indeterminate forms + if (base.dir == dir::finite_positive) + return {dir::finite_positive}; + if (base.dir == dir::zero || is_infinite(base.dir)) + return {dir::indeterminate}; + return {dir::unknown}; } bool exp_positive = is_positive(exponent.dir); @@ -220,7 +230,9 @@ limit_result limit_algebra::apply_pow(limit_result base, // 0^(+c) = 0 return {dir::zero}; } else { - // 0^(-c) = +inf (polynomial) + // 0^(-c) = +inf (polynomial) when the base stays positive + if (!zero_from_above) + return {dir::unknown}; return {dir::pos_infinity, {gtype::polynomial, 1.0}}; } case dir::finite_positive: @@ -301,13 +313,16 @@ limit_result limit_algebra::apply_abs(limit_result a) { } } -limit_result limit_algebra::apply_reciprocal(limit_result a) { +limit_result limit_algebra::apply_reciprocal(limit_result a, + bool zero_from_above) { if (a.dir == dir::indeterminate || a.dir == dir::unknown) return a; switch (a.dir) { case dir::zero: - // 1/0 = +inf (polynomial, degree 1) + // 1/0+ = +inf (polynomial, degree 1) + if (!zero_from_above) + return {dir::unknown}; return {dir::pos_infinity, {gtype::polynomial, 1.0}}; case dir::finite_positive: return {dir::finite_positive}; diff --git a/src/numsim_cas/scalar/visitors/scalar_limit_visitor.cpp b/src/numsim_cas/scalar/visitors/scalar_limit_visitor.cpp index f1bff02f..784cb31c 100644 --- a/src/numsim_cas/scalar/visitors/scalar_limit_visitor.cpp +++ b/src/numsim_cas/scalar/visitors/scalar_limit_visitor.cpp @@ -1,5 +1,6 @@ #include +#include #include #include #include @@ -30,6 +31,18 @@ limit_result target_to_limit(limit_target target) { } // namespace +bool scalar_limit_visitor::zero_from_above(expr_holder_t const &expr) const { + if (expr == m_limit_var) + return m_target.target == limit_target::point::zero_plus; + return is_positive(expr); +} + +bool scalar_limit_visitor::zero_from_below(expr_holder_t const &expr) const { + if (expr == m_limit_var) + return m_target.target == limit_target::point::zero_minus; + return is_negative(expr); +} + limit_result scalar_limit_visitor::apply(expr_holder_t const &expr) { if (!expr.is_valid()) return {dir::zero}; @@ -44,10 +57,15 @@ limit_result scalar_limit_visitor::apply(expr_holder_t const &expr) { // ─── Leaf nodes ─────────────────────────────────────────────────── -void scalar_limit_visitor::operator()([[maybe_unused]] scalar const &) { - // If we reach here, this symbol is NOT the limit variable - // (the limit variable case is handled in apply() before dispatching) - m_result = {dir::finite_positive}; +void scalar_limit_visitor::operator()(scalar const &v) { + // Not the limit variable (handled in apply()), so a constant of unknown sign + // unless its assumptions say otherwise. + if (v.assumptions().contains(positive{})) + m_result = {dir::finite_positive}; + else if (v.assumptions().contains(negative{})) + m_result = {dir::finite_negative}; + else + m_result = {dir::unknown}; } void scalar_limit_visitor::operator()([[maybe_unused]] scalar_zero const &) { @@ -109,43 +127,42 @@ void scalar_limit_visitor::operator()(scalar_negative const &v) { } void scalar_limit_visitor::operator()(scalar_pow const &v) { - m_result = apply_pow(apply(v.expr_lhs()), apply(v.expr_rhs())); + m_result = apply_pow(apply(v.expr_lhs()), apply(v.expr_rhs()), + zero_from_above(v.expr_lhs())); } // ─── Functions ──────────────────────────────────────────────────── +// For a finite nonzero argument sin, cos and tan can take either sign; at +// infinity sin and cos oscillate. void scalar_limit_visitor::operator()(scalar_sin const &v) { auto child = apply(v.expr()); - if (child.dir == dir::indeterminate || child.dir == dir::unknown) { + if (child.dir == dir::indeterminate) m_result = child; - } else if (child.dir == dir::pos_infinity || child.dir == dir::neg_infinity) { - // sin oscillates => indeterminate + else if (child.dir == dir::zero) + m_result = {dir::zero}; + else m_result = {dir::unknown}; - } else { - // finite input => bounded output in [-1, 1] - m_result = {dir::finite_positive}; - } } void scalar_limit_visitor::operator()(scalar_cos const &v) { auto child = apply(v.expr()); - if (child.dir == dir::indeterminate || child.dir == dir::unknown) { + if (child.dir == dir::indeterminate) m_result = child; - } else if (child.dir == dir::pos_infinity || child.dir == dir::neg_infinity) { - m_result = {dir::unknown}; - } else { + else if (child.dir == dir::zero) m_result = {dir::finite_positive}; - } + else + m_result = {dir::unknown}; } void scalar_limit_visitor::operator()(scalar_tan const &v) { auto child = apply(v.expr()); - if (child.dir == dir::indeterminate || child.dir == dir::unknown) { + if (child.dir == dir::indeterminate) m_result = child; - } else { - // tan can diverge at pi/2 + n*pi, treat as unknown in general + else if (child.dir == dir::zero) + m_result = {dir::zero}; + else m_result = {dir::unknown}; - } } void scalar_limit_visitor::operator()(scalar_exp const &v) { @@ -153,7 +170,7 @@ void scalar_limit_visitor::operator()(scalar_exp const &v) { } void scalar_limit_visitor::operator()(scalar_log const &v) { - m_result = apply_log(apply(v.expr())); + m_result = apply_log(apply(v.expr()), zero_from_above(v.expr())); } void scalar_limit_visitor::operator()(scalar_sqrt const &v) { @@ -166,43 +183,90 @@ void scalar_limit_visitor::operator()(scalar_abs const &v) { void scalar_limit_visitor::operator()(scalar_sign const &v) { auto child = apply(v.expr()); - if (child.dir == dir::indeterminate || child.dir == dir::unknown) { + switch (child.dir) { + case dir::indeterminate: m_result = child; - } else { - // sign is bounded in {-1, 0, 1} + break; + case dir::zero: + // sign jumps at 0: the limit is 1 only when the argument stays positive + if (zero_from_above(v.expr())) + m_result = {dir::finite_positive}; + else if (zero_from_below(v.expr())) + m_result = {dir::finite_negative}; + else + m_result = {dir::unknown}; + break; + case dir::finite_positive: + case dir::pos_infinity: m_result = {dir::finite_positive}; + break; + case dir::finite_negative: + case dir::neg_infinity: + m_result = {dir::finite_negative}; + break; + default: + m_result = {dir::unknown}; } } +// asin and acos are undefined outside [-1, 1], so infinite arguments are +// unknown. void scalar_limit_visitor::operator()(scalar_asin const &v) { auto child = apply(v.expr()); - if (child.dir == dir::indeterminate || child.dir == dir::unknown) { + switch (child.dir) { + case dir::indeterminate: m_result = child; - } else { - // asin is bounded [-pi/2, pi/2] for finite input + break; + case dir::zero: + m_result = {dir::zero}; + break; + case dir::finite_positive: m_result = {dir::finite_positive}; + break; + case dir::finite_negative: + m_result = {dir::finite_negative}; + break; + default: + m_result = {dir::unknown}; } } void scalar_limit_visitor::operator()(scalar_acos const &v) { auto child = apply(v.expr()); - if (child.dir == dir::indeterminate || child.dir == dir::unknown) { + switch (child.dir) { + case dir::indeterminate: m_result = child; - } else { + break; + case dir::zero: + case dir::finite_negative: + // acos maps [-1, 0] into [pi/2, pi] m_result = {dir::finite_positive}; + break; + default: + // acos(1) = 0, so a positive argument leaves the sign open + m_result = {dir::unknown}; } } void scalar_limit_visitor::operator()(scalar_atan const &v) { auto child = apply(v.expr()); - if (child.dir == dir::indeterminate || child.dir == dir::unknown) { + switch (child.dir) { + case dir::indeterminate: m_result = child; - } else if (child.dir == dir::neg_infinity) { - // atan(-inf) = -pi/2 - m_result = {dir::finite_negative}; - } else { - // atan(finite) or atan(+inf) = finite_positive + break; + case dir::zero: + m_result = {dir::zero}; + break; + case dir::finite_positive: + case dir::pos_infinity: m_result = {dir::finite_positive}; + break; + case dir::finite_negative: + case dir::neg_infinity: + m_result = {dir::finite_negative}; + break; + default: + m_result = {dir::unknown}; } } diff --git a/src/numsim_cas/tensor_to_scalar/visitors/tensor_to_scalar_limit_visitor.cpp b/src/numsim_cas/tensor_to_scalar/visitors/tensor_to_scalar_limit_visitor.cpp index e5dd0578..a1907331 100644 --- a/src/numsim_cas/tensor_to_scalar/visitors/tensor_to_scalar_limit_visitor.cpp +++ b/src/numsim_cas/tensor_to_scalar/visitors/tensor_to_scalar_limit_visitor.cpp @@ -1,6 +1,7 @@ #include #include +#include #include #include @@ -8,6 +9,20 @@ namespace numsim::cas { using dir = limit_result::direction; +namespace { + +// A value that does not depend on the limit variable: only its assumptions +// can give it a sign. +limit_result constant_sign(tensor_to_scalar_expression const &e) { + if (e.assumptions().contains(positive{})) + return {dir::finite_positive}; + if (e.assumptions().contains(negative{})) + return {dir::finite_negative}; + return {dir::unknown}; +} + +} // namespace + // ─── Constructors ───────────────────────────────────────────────── tensor_to_scalar_limit_visitor::tensor_to_scalar_limit_visitor( @@ -65,10 +80,16 @@ bool tensor_to_scalar_limit_visitor::depends_on_limit_var( return depends_on_tensor(expr, m_tensor_var); } +bool tensor_to_scalar_limit_visitor::zero_from_above( + t2s_holder_t const &expr) const { + if (m_mode == dependency_mode::exact_match && expr == m_limit_var_t2s) + return m_target.target == limit_target::point::zero_plus; + return expr.get().assumptions().contains(positive{}); +} + // ─── T2S functions ──────────────────────────────────────────────── -void tensor_to_scalar_limit_visitor::operator()( - [[maybe_unused]] tensor_trace const &v) { +void tensor_to_scalar_limit_visitor::operator()(tensor_trace const &v) { // trace depends on tensor child if (m_mode == dependency_mode::tensor_dependency) { if (contains_expression(v.expr(), m_tensor_var)) { @@ -77,22 +98,20 @@ void tensor_to_scalar_limit_visitor::operator()( return; } } - m_result = {dir::finite_positive}; + m_result = constant_sign(v); } -void tensor_to_scalar_limit_visitor::operator()( - [[maybe_unused]] tensor_dot const &v) { +void tensor_to_scalar_limit_visitor::operator()(tensor_dot const &v) { if (m_mode == dependency_mode::tensor_dependency) { if (contains_expression(v.expr(), m_tensor_var)) { m_result = {dir::unknown}; return; } } - m_result = {dir::finite_positive}; + m_result = constant_sign(v); } -void tensor_to_scalar_limit_visitor::operator()( - [[maybe_unused]] tensor_det const &v) { +void tensor_to_scalar_limit_visitor::operator()(tensor_det const &v) { if (m_mode == dependency_mode::tensor_dependency) { if (contains_expression(v.expr(), m_tensor_var)) { // det depends on tensor: behavior depends on limit target @@ -101,11 +120,10 @@ void tensor_to_scalar_limit_visitor::operator()( return; } } - m_result = {dir::finite_positive}; + m_result = constant_sign(v); } -void tensor_to_scalar_limit_visitor::operator()( - [[maybe_unused]] tensor_norm const &v) { +void tensor_to_scalar_limit_visitor::operator()(tensor_norm const &v) { if (m_mode == dependency_mode::tensor_dependency) { if (contains_expression(v.expr(), m_tensor_var)) { // norm(F) as F -> infinity => +infinity (polynomial) @@ -118,11 +136,11 @@ void tensor_to_scalar_limit_visitor::operator()( return; } } - m_result = {dir::finite_positive}; + m_result = constant_sign(v); } void tensor_to_scalar_limit_visitor::operator()( - [[maybe_unused]] tensor_to_scalar_eigenvalue const &v) { + tensor_to_scalar_eigenvalue const &v) { // An eigenvalue's sign and magnitude aren't recoverable from the AST // generically (unlike norm/det), so any tensor dependency is unknown. if (m_mode == dependency_mode::tensor_dependency && @@ -130,11 +148,11 @@ void tensor_to_scalar_limit_visitor::operator()( m_result = {dir::unknown}; return; } - m_result = {dir::finite_positive}; + m_result = constant_sign(v); } void tensor_to_scalar_limit_visitor::operator()( - [[maybe_unused]] tensor_to_scalar_divided_difference const &v) { + tensor_to_scalar_divided_difference const &v) { // A divided difference of eigenvalues — like an eigenvalue, not // recoverable from the AST; any tensor dependency is unknown. if (m_mode == dependency_mode::tensor_dependency && @@ -142,7 +160,7 @@ void tensor_to_scalar_limit_visitor::operator()( m_result = {dir::unknown}; return; } - m_result = {dir::finite_positive}; + m_result = constant_sign(v); } void tensor_to_scalar_limit_visitor::operator()( @@ -155,7 +173,7 @@ void tensor_to_scalar_limit_visitor::operator()( return; } } - m_result = {dir::finite_positive}; + m_result = constant_sign(v); } // ─── Arithmetic ─────────────────────────────────────────────────── @@ -188,11 +206,12 @@ void tensor_to_scalar_limit_visitor::operator()( } void tensor_to_scalar_limit_visitor::operator()(tensor_to_scalar_pow const &v) { - m_result = apply_pow(apply(v.expr_lhs()), apply(v.expr_rhs())); + m_result = apply_pow(apply(v.expr_lhs()), apply(v.expr_rhs()), + zero_from_above(v.expr_lhs())); } void tensor_to_scalar_limit_visitor::operator()(tensor_to_scalar_log const &v) { - m_result = apply_log(apply(v.expr())); + m_result = apply_log(apply(v.expr()), zero_from_above(v.expr())); } void tensor_to_scalar_limit_visitor::operator()(tensor_to_scalar_exp const &v) { @@ -217,12 +236,17 @@ void tensor_to_scalar_limit_visitor::operator()( } void tensor_to_scalar_limit_visitor::operator()( - tensor_to_scalar_scalar_wrapper const & /*v*/) { - // Delegate to scalar limit visitor - // Scalar expressions don't depend on tensor variables, - // so they should evaluate to finite - // Scalar sub-expressions don't depend on tensor or T2S limit variables - m_result = {dir::finite_positive}; + tensor_to_scalar_scalar_wrapper const &v) { + // A scalar sub-expression is constant with respect to the limit variable. + auto const &e = v.expr(); + if (is_same(e)) + m_result = {dir::zero}; + else if (is_positive(e)) + m_result = {dir::finite_positive}; + else if (is_negative(e)) + m_result = {dir::finite_negative}; + else + m_result = {dir::unknown}; } // if_then_else (#135 / #210): limit depends on the condition's eventual diff --git a/tests/LimitVisitorTest.h b/tests/LimitVisitorTest.h index af79925f..b1c4ab9e 100644 --- a/tests/LimitVisitorTest.h +++ b/tests/LimitVisitorTest.h @@ -287,6 +287,143 @@ TEST(ScalarLimit, NegationFlips) { EXPECT_EQ(result.dir, dir::neg_infinity); } +// ─── Signs are reported only when provable ───────────────────────── + +struct limit_algebra_probe : limit_algebra { + using limit_algebra::apply_log; + using limit_algebra::apply_pow; + using limit_algebra::apply_reciprocal; + using limit_algebra::combine_add; + using limit_algebra::combine_mul; +}; + +TEST(LimitAlgebra, AddSigns) { + using A = limit_algebra_probe; + struct row { + dir a, b, expected; + }; + for (auto [a, b, expected] : { + row{dir::finite_positive, dir::finite_positive, + dir::finite_positive}, + row{dir::finite_negative, dir::finite_negative, + dir::finite_negative}, + row{dir::finite_positive, dir::finite_negative, dir::unknown}, + row{dir::finite_negative, dir::finite_positive, dir::unknown}, + row{dir::zero, dir::finite_negative, dir::finite_negative}, + row{dir::pos_infinity, dir::finite_negative, dir::pos_infinity}, + row{dir::pos_infinity, dir::neg_infinity, dir::indeterminate}, + }) { + EXPECT_EQ(A::combine_add({a}, {b}).dir, expected) + << static_cast(a) << " + " << static_cast(b); + } +} + +TEST(LimitAlgebra, MulSigns) { + using A = limit_algebra_probe; + struct row { + dir a, b, expected; + }; + for (auto [a, b, expected] : { + row{dir::finite_positive, dir::finite_negative, + dir::finite_negative}, + row{dir::finite_negative, dir::finite_negative, + dir::finite_positive}, + row{dir::finite_negative, dir::pos_infinity, dir::neg_infinity}, + row{dir::zero, dir::pos_infinity, dir::indeterminate}, + row{dir::unknown, dir::pos_infinity, dir::unknown}, + }) { + EXPECT_EQ(A::combine_mul({a}, {b}).dir, expected) + << static_cast(a) << " * " << static_cast(b); + } +} + +TEST(LimitAlgebra, ZeroNeedsApproachSide) { + using A = limit_algebra_probe; + limit_result zero{dir::zero}, neg{dir::finite_negative}; + EXPECT_EQ(A::apply_reciprocal(zero, false).dir, dir::unknown); + EXPECT_EQ(A::apply_reciprocal(zero, true).dir, dir::pos_infinity); + EXPECT_EQ(A::apply_pow(zero, neg, false).dir, dir::unknown); + EXPECT_EQ(A::apply_pow(zero, neg, true).dir, dir::pos_infinity); + EXPECT_EQ(A::apply_log(zero, false).dir, dir::unknown); + EXPECT_EQ(A::apply_log(zero, true).dir, dir::neg_infinity); +} + +TEST(LimitAlgebra, LogAndPowOfFiniteValues) { + using A = limit_algebra_probe; + // log(c) < 0 for c < 1 + EXPECT_EQ(A::apply_log({dir::finite_positive}, false).dir, dir::unknown); + // 0^0 and inf^0 are indeterminate forms + EXPECT_EQ(A::apply_pow({dir::zero}, {dir::zero}, true).dir, + dir::indeterminate); + EXPECT_EQ(A::apply_pow({dir::pos_infinity}, {dir::zero}, false).dir, + dir::indeterminate); + EXPECT_EQ(A::apply_pow({dir::finite_positive}, {dir::zero}, false).dir, + dir::finite_positive); +} + +TEST(ScalarLimit, OddFunctionsAtZero) { + auto x = make_expression("x"); + scalar_limit_visitor v(x, {pt::zero_plus}); + EXPECT_EQ(v.apply(sin(x)).dir, dir::zero); + EXPECT_EQ(v.apply(tan(x)).dir, dir::zero); + EXPECT_EQ(v.apply(asin(x)).dir, dir::zero); + EXPECT_EQ(v.apply(atan(x)).dir, dir::zero); + EXPECT_EQ(v.apply(cos(x)).dir, dir::finite_positive); + EXPECT_EQ(v.apply(acos(x)).dir, dir::finite_positive); +} + +TEST(ScalarLimit, ReciprocalOfSinAtZeroIsNotFinite) { + auto x = make_expression("x"); + scalar_limit_visitor v(x, {pt::zero_plus}); + auto r = v.apply(pow(sin(x), make_scalar_constant(-1))).dir; + EXPECT_NE(r, dir::finite_positive); + EXPECT_NE(r, dir::finite_negative); + EXPECT_NE(r, dir::zero); +} + +TEST(ScalarLimit, ZeroFromBelowHasNoLogOrReciprocal) { + auto x = make_expression("x"); + scalar_limit_visitor v(x, {pt::zero_minus}); + EXPECT_EQ(v.apply(log(x)).dir, dir::unknown); + EXPECT_EQ(v.apply(pow(x, make_scalar_constant(-1))).dir, dir::unknown); +} + +TEST(ScalarLimit, SignAtZeroFollowsApproachSide) { + auto x = make_expression("x"); + scalar_limit_visitor right(x, {pt::zero_plus}); + scalar_limit_visitor left(x, {pt::zero_minus}); + EXPECT_EQ(right.apply(sign(x)).dir, dir::finite_positive); + EXPECT_EQ(left.apply(sign(x)).dir, dir::finite_negative); +} + +TEST(ScalarLimit, AssumedSignOfOtherSymbols) { + auto x = make_expression("x"); + auto [a, b, c, y] = make_scalar_variable("a", "b", "c", "y"); + a.assumption(negative{}); + c.assumption(negative{}); + b.assumption(positive{}); + scalar_limit_visitor v(x, {pt::pos_infinity}); + EXPECT_EQ(v.apply(a).dir, dir::finite_negative); + EXPECT_EQ(v.apply(b).dir, dir::finite_positive); + EXPECT_EQ(v.apply(y).dir, dir::unknown); + EXPECT_EQ(v.apply(a + c).dir, dir::finite_negative); + EXPECT_EQ(v.apply(a + b).dir, dir::unknown); + EXPECT_EQ(v.apply(a * x).dir, dir::neg_infinity); +} + +TEST(ScalarLimit, InverseTrigSigns) { + auto x = make_expression("x"); + auto [a] = make_scalar_variable("a"); + a.assumption(negative{}); + scalar_limit_visitor v(x, {pt::pos_infinity}); + EXPECT_EQ(v.apply(atan(a)).dir, dir::finite_negative); + EXPECT_EQ(v.apply(asin(a)).dir, dir::finite_negative); + EXPECT_EQ(v.apply(acos(a)).dir, dir::finite_positive); + // asin is undefined outside [-1, 1] + EXPECT_EQ(v.apply(asin(x)).dir, dir::unknown); + EXPECT_EQ(v.apply(cos(a)).dir, dir::unknown); +} + // ═══════════════════════════════════════════════════════════════════ // T2S limit visitor tests (exact match mode) // ═══════════════════════════════════════════════════════════════════ @@ -396,6 +533,30 @@ TEST(T2sLimit, TensorDepScalarWrapperFinite) { EXPECT_EQ(result.dir, dir::finite_positive); } +TEST(T2sLimit, ConstantSignsComeFromValuesAndAssumptions) { + auto F = make_expression("F", 3, 2); + auto G = make_expression("G", 3, 2); + auto P = make_expression("P", 3, 2); + P.assumption(positive_definite{}); + tensor_to_scalar_limit_visitor v(F, {pt::pos_infinity}); + auto neg = make_expression( + make_expression(-2.0)); + EXPECT_EQ(v.apply(neg).dir, dir::finite_negative); + EXPECT_EQ(v.apply(trace(G)).dir, dir::unknown); + EXPECT_EQ(v.apply(det(P)).dir, dir::finite_positive); +} + +TEST(T2sLimit, ExactMatchReciprocalNeedsApproachSide) { + auto F = make_expression("F", 3, 2); + auto J = det(F); + auto inv_J = pow(J, make_scalar_constant(-1)); + tensor_to_scalar_limit_visitor right(J, {pt::zero_plus}); + tensor_to_scalar_limit_visitor left(J, {pt::zero_minus}); + EXPECT_EQ(right.apply(inv_J).dir, dir::pos_infinity); + EXPECT_EQ(left.apply(inv_J).dir, dir::unknown); + EXPECT_EQ(left.apply(log(J)).dir, dir::unknown); +} + // ═══════════════════════════════════════════════════════════════════ // Growth rate tracking // ═══════════════════════════════════════════════════════════════════ From b6f32cdf6c9568cf97cb2be347409d70e3b61eea Mon Sep 17 00:00:00 2001 From: petlenz Date: Tue, 15 Sep 2026 23:32:28 +0200 Subject: [PATCH 2/5] Correct two limit tests that asserted an unprovable positive sign ScalarLimit.OtherVariableFinite: y has no assumptions, so as x -> 0+ it is a constant of either sign; finite_positive is not provable (y = -1 is a counterexample). Expected result is now unknown. T2sLimit.TensorDepIndependentExpr: det(G) is independent of F but a general 3x3 determinant takes either sign (det(-I) = -1). Expected result is now unknown. Signed-off-by: petlenz --- tests/LimitVisitorTest.h | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/tests/LimitVisitorTest.h b/tests/LimitVisitorTest.h index b1c4ab9e..9cb6f427 100644 --- a/tests/LimitVisitorTest.h +++ b/tests/LimitVisitorTest.h @@ -139,7 +139,8 @@ TEST(ScalarLimit, OtherVariableFinite) { auto y = make_expression("y"); scalar_limit_visitor v(x, {pt::zero_plus}); auto result = v.apply(y); - EXPECT_EQ(result.dir, dir::finite_positive); + // y is constant in x but may have either sign + EXPECT_EQ(result.dir, dir::unknown); } // ─── log(x) as x -> 0+ ──────────────────────────────────────────── @@ -509,10 +510,10 @@ TEST(T2sLimit, TensorDepNormToPosInfinity) { TEST(T2sLimit, TensorDepIndependentExpr) { auto F = make_expression("F", 3, 2); auto G = make_expression("G", 3, 2); - auto expr = det(G); // independent of F + auto expr = det(G); // independent of F, but det(G) may have either sign tensor_to_scalar_limit_visitor v(F, {pt::pos_infinity}); auto result = v.apply(expr); - EXPECT_EQ(result.dir, dir::finite_positive); + EXPECT_EQ(result.dir, dir::unknown); } TEST(T2sLimit, TensorDepDetIsUnknown) { From 36b9569a243210e5c3c0551e50581a9b98acf1b5 Mon Sep 17 00:00:00 2001 From: petlenz Date: Wed, 16 Sep 2026 22:02:01 +0200 Subject: [PATCH 3/5] Collapse apply_log's identical unknown branches into the default Signed-off-by: petlenz --- src/numsim_cas/core/limit_algebra.cpp | 11 ++--------- 1 file changed, 2 insertions(+), 9 deletions(-) diff --git a/src/numsim_cas/core/limit_algebra.cpp b/src/numsim_cas/core/limit_algebra.cpp index 213aea25..ddaad0b2 100644 --- a/src/numsim_cas/core/limit_algebra.cpp +++ b/src/numsim_cas/core/limit_algebra.cpp @@ -187,19 +187,12 @@ limit_result limit_algebra::apply_log(limit_result a, bool zero_from_above) { if (!zero_from_above) return {dir::unknown}; return {dir::neg_infinity, {gtype::logarithmic, 1.0}}; - case dir::finite_positive: - // log(c) is negative for c < 1 and zero at c = 1 - return {dir::unknown}; - case dir::finite_negative: - // log(negative) is undefined in reals - return {dir::unknown}; case dir::pos_infinity: // log(+inf) = +inf (logarithmic -- slower than any polynomial) return {dir::pos_infinity, {gtype::logarithmic, 1.0}}; - case dir::neg_infinity: - // log(-inf) undefined in reals - [[fallthrough]]; default: + // log(c) has no provable sign for finite positive c (negative below 1, + // zero at 1); for a negative value or -inf it is undefined in the reals return {dir::unknown}; } } From d1da09beedde9879e2d4cf74263e5eddabefebab Mon Sep 17 00:00:00 2001 From: petlenz Date: Wed, 16 Sep 2026 22:25:10 +0200 Subject: [PATCH 4/5] Review fixes on #357: domain guards, nonnegative arguments, c^0 - asin and acos claim a sign only for an argument tending to zero or one provably within [-1, 1]; asin(x + 2) is NaN, not finite positive - a zero limit counts as approached from above when the expression is nonnegative by construction (abs, exp, sqrt, even integer power) or assumed nonnegative, which restores 1/x^2, 1/|x| and log(|x|) - an even integer exponent makes a power nonnegative from either side - c^0 is 1 for any finite nonzero c, negative bases included Signed-off-by: petlenz --- .../scalar/visitors/scalar_limit_visitor.h | 1 + src/numsim_cas/core/limit_algebra.cpp | 4 +- .../scalar/visitors/scalar_limit_visitor.cpp | 109 ++++++++++++------ tests/LimitVisitorTest.h | 49 +++++++- 4 files changed, 123 insertions(+), 40 deletions(-) diff --git a/include/numsim_cas/scalar/visitors/scalar_limit_visitor.h b/include/numsim_cas/scalar/visitors/scalar_limit_visitor.h index 56e55cdc..2bd6895a 100644 --- a/include/numsim_cas/scalar/visitors/scalar_limit_visitor.h +++ b/include/numsim_cas/scalar/visitors/scalar_limit_visitor.h @@ -59,6 +59,7 @@ class scalar_limit_visitor final : public scalar_visitor_const_t, private: bool zero_from_above(expr_holder_t const &expr) const; bool zero_from_below(expr_holder_t const &expr) const; + void inverse_trig(expr_holder_t const &arg, bool is_asin); expr_holder_t m_limit_var; limit_target m_target; diff --git a/src/numsim_cas/core/limit_algebra.cpp b/src/numsim_cas/core/limit_algebra.cpp index ddaad0b2..542b1c55 100644 --- a/src/numsim_cas/core/limit_algebra.cpp +++ b/src/numsim_cas/core/limit_algebra.cpp @@ -207,8 +207,8 @@ limit_result limit_algebra::apply_pow(limit_result base, limit_result exponent, // Only handle constant/finite exponents for now if (is_finite(exponent.dir) || exponent.dir == dir::zero) { if (exponent.dir == dir::zero) { - // c^0 = 1, but 0^0 and inf^0 are indeterminate forms - if (base.dir == dir::finite_positive) + // c^0 = 1 for any finite nonzero c; 0^0 and inf^0 are indeterminate + if (base.dir == dir::finite_positive || base.dir == dir::finite_negative) return {dir::finite_positive}; if (base.dir == dir::zero || is_infinite(base.dir)) return {dir::indeterminate}; diff --git a/src/numsim_cas/scalar/visitors/scalar_limit_visitor.cpp b/src/numsim_cas/scalar/visitors/scalar_limit_visitor.cpp index 784cb31c..c3213542 100644 --- a/src/numsim_cas/scalar/visitors/scalar_limit_visitor.cpp +++ b/src/numsim_cas/scalar/visitors/scalar_limit_visitor.cpp @@ -1,6 +1,7 @@ #include #include +#include #include #include #include @@ -15,6 +16,43 @@ scalar_limit_visitor::scalar_limit_visitor(expr_holder_t const &limit_var, namespace { +// Nonnegative by construction, so a limit of zero is approached from above. +bool is_structurally_nonnegative( + expression_holder const &e) { + if (is_same(e) || is_same(e) || + is_same(e)) + return true; + if (is_same(e)) { + auto exponent = try_int_constant(e.get().expr_rhs()); + return exponent && *exponent % 2 == 0; + } + return false; +} + +// asin and acos are real only on [-1, 1]; outside it the value is NaN. +bool is_within_unit_interval(expression_holder const &e) { + if (is_same(e) || is_same(e) || + is_same(e)) + return true; + auto value = domain_traits::try_numeric(e); + if (!value) + return false; + auto magnitude = std::visit( + [](auto const &x) -> double { + using V = std::decay_t; + if constexpr (std::is_same_v>) { + return std::abs(x.real()); + } else if constexpr (std::is_same_v) { + return std::abs(static_cast(x.num) / + static_cast(x.den)); + } else { + return std::abs(static_cast(x)); + } + }, + value->raw()); + return magnitude <= 1.0; +} + limit_result target_to_limit(limit_target target) { using pt = limit_target::point; switch (target.target) { @@ -34,7 +72,8 @@ limit_result target_to_limit(limit_target target) { bool scalar_limit_visitor::zero_from_above(expr_holder_t const &expr) const { if (expr == m_limit_var) return m_target.target == limit_target::point::zero_plus; - return is_positive(expr); + return is_positive(expr) || is_nonnegative(expr) || + is_structurally_nonnegative(expr); } bool scalar_limit_visitor::zero_from_below(expr_holder_t const &expr) const { @@ -127,8 +166,11 @@ void scalar_limit_visitor::operator()(scalar_negative const &v) { } void scalar_limit_visitor::operator()(scalar_pow const &v) { - m_result = apply_pow(apply(v.expr_lhs()), apply(v.expr_rhs()), - zero_from_above(v.expr_lhs())); + // An even integer exponent makes the power nonnegative from either side. + auto exponent = try_int_constant(v.expr_rhs()); + bool from_above = + zero_from_above(v.expr_lhs()) || (exponent && *exponent % 2 == 0); + m_result = apply_pow(apply(v.expr_lhs()), apply(v.expr_rhs()), from_above); } // ─── Functions ──────────────────────────────────────────────────── @@ -209,43 +251,42 @@ void scalar_limit_visitor::operator()(scalar_sign const &v) { } } -// asin and acos are undefined outside [-1, 1], so infinite arguments are -// unknown. -void scalar_limit_visitor::operator()(scalar_asin const &v) { - auto child = apply(v.expr()); - switch (child.dir) { - case dir::indeterminate: +// asin and acos are real only on [-1, 1], so a sign needs an argument that is +// provably in range: asin(0) = 0, acos(0) = pi/2. +void scalar_limit_visitor::inverse_trig(expr_holder_t const &arg, + bool is_asin) { + auto child = apply(arg); + if (child.dir == dir::indeterminate) { m_result = child; - break; - case dir::zero: - m_result = {dir::zero}; - break; - case dir::finite_positive: - m_result = {dir::finite_positive}; - break; - case dir::finite_negative: - m_result = {dir::finite_negative}; - break; - default: + return; + } + if (child.dir == dir::zero) { + m_result = + is_asin ? limit_result{dir::zero} : limit_result{dir::finite_positive}; + return; + } + if (!is_within_unit_interval(arg)) { m_result = {dir::unknown}; + return; + } + if (child.dir == dir::finite_negative) { + // asin maps [-1, 0) below zero; acos maps it into (pi/2, pi] + m_result = {is_asin ? dir::finite_negative : dir::finite_positive}; + return; } + // asin(c) > 0 for c in (0, 1]; acos(1) = 0 leaves acos open + if (is_asin && child.dir == dir::finite_positive) + m_result = {dir::finite_positive}; + else + m_result = {dir::unknown}; +} + +void scalar_limit_visitor::operator()(scalar_asin const &v) { + inverse_trig(v.expr(), true); } void scalar_limit_visitor::operator()(scalar_acos const &v) { - auto child = apply(v.expr()); - switch (child.dir) { - case dir::indeterminate: - m_result = child; - break; - case dir::zero: - case dir::finite_negative: - // acos maps [-1, 0] into [pi/2, pi] - m_result = {dir::finite_positive}; - break; - default: - // acos(1) = 0, so a positive argument leaves the sign open - m_result = {dir::unknown}; - } + inverse_trig(v.expr(), false); } void scalar_limit_visitor::operator()(scalar_atan const &v) { diff --git a/tests/LimitVisitorTest.h b/tests/LimitVisitorTest.h index 9cb6f427..d363ccfa 100644 --- a/tests/LimitVisitorTest.h +++ b/tests/LimitVisitorTest.h @@ -360,6 +360,9 @@ TEST(LimitAlgebra, LogAndPowOfFiniteValues) { dir::indeterminate); EXPECT_EQ(A::apply_pow({dir::finite_positive}, {dir::zero}, false).dir, dir::finite_positive); + // (-2)^0 = 1 + EXPECT_EQ(A::apply_pow({dir::finite_negative}, {dir::zero}, false).dir, + dir::finite_positive); } TEST(ScalarLimit, OddFunctionsAtZero) { @@ -416,15 +419,53 @@ TEST(ScalarLimit, InverseTrigSigns) { auto x = make_expression("x"); auto [a] = make_scalar_variable("a"); a.assumption(negative{}); + auto half = make_scalar_constant(-0.5); scalar_limit_visitor v(x, {pt::pos_infinity}); EXPECT_EQ(v.apply(atan(a)).dir, dir::finite_negative); - EXPECT_EQ(v.apply(asin(a)).dir, dir::finite_negative); - EXPECT_EQ(v.apply(acos(a)).dir, dir::finite_positive); - // asin is undefined outside [-1, 1] - EXPECT_EQ(v.apply(asin(x)).dir, dir::unknown); + EXPECT_EQ(v.apply(asin(half)).dir, dir::finite_negative); + EXPECT_EQ(v.apply(acos(half)).dir, dir::finite_positive); + EXPECT_EQ(v.apply(asin(sign(a))).dir, dir::finite_negative); EXPECT_EQ(v.apply(cos(a)).dir, dir::unknown); } +// asin and acos are NaN outside [-1, 1], so an argument of unknown magnitude +// has no provable sign. +TEST(ScalarLimit, InverseTrigOutsideDomainIsUnknown) { + auto x = make_expression("x"); + auto [a] = make_scalar_variable("a"); + a.assumption(negative{}); + scalar_limit_visitor v(x, {pt::zero_plus}); + EXPECT_EQ(v.apply(asin(x + make_scalar_constant(2))).dir, dir::unknown); + EXPECT_EQ(v.apply(acos(x - make_scalar_constant(2))).dir, dir::unknown); + EXPECT_EQ(v.apply(asin(a)).dir, dir::unknown); + EXPECT_EQ(v.apply(acos(a)).dir, dir::unknown); + scalar_limit_visitor at_infinity(x, {pt::pos_infinity}); + EXPECT_EQ(at_infinity.apply(asin(x)).dir, dir::unknown); + // in range: asin(0) = 0 and acos(0) = pi/2 need no magnitude bound + EXPECT_EQ(v.apply(asin(x)).dir, dir::zero); + EXPECT_EQ(v.apply(acos(x)).dir, dir::finite_positive); +} + +// A structurally nonnegative argument reaches zero from above, whichever side +// the limit variable approaches from. +TEST(ScalarLimit, NonnegativeArgumentsKeepTheirZeroSide) { + auto x = make_expression("x"); + scalar_limit_visitor v(x, {pt::zero_minus}); + EXPECT_EQ( + v.apply(pow(pow(x, make_scalar_constant(2)), make_scalar_constant(-1))) + .dir, + dir::pos_infinity); + EXPECT_EQ(v.apply(pow(abs(x), make_scalar_constant(-1))).dir, + dir::pos_infinity); + EXPECT_EQ(v.apply(log(abs(x))).dir, dir::neg_infinity); + EXPECT_EQ(v.apply(log(sqrt(abs(x)))).dir, dir::neg_infinity); + // an odd power keeps the sign of x, so the side stays unknown + EXPECT_EQ( + v.apply(pow(pow(x, make_scalar_constant(3)), make_scalar_constant(-1))) + .dir, + dir::unknown); +} + // ═══════════════════════════════════════════════════════════════════ // T2S limit visitor tests (exact match mode) // ═══════════════════════════════════════════════════════════════════ From c54430e026a061238813be9f88d5fe29de9d1808 Mon Sep 17 00:00:00 2001 From: petlenz Date: Wed, 16 Sep 2026 23:02:22 +0200 Subject: [PATCH 5/5] Review fixes on #357: roots and powers that are undefined near zero - sqrt of an argument tending to zero is unknown unless the approach is provably from above; sqrt is nonnegative only where it is defined, so it reaches zero from above exactly when its operand does - a non-integer power of a base that may reach zero from below is NaN - c^t for c < 0 and t -> 0 has no real limit, so the zero-exponent branch claims a value only for a positive base Signed-off-by: petlenz --- include/numsim_cas/core/limit_algebra.h | 2 +- src/numsim_cas/core/limit_algebra.cpp | 10 ++-- .../scalar/visitors/scalar_limit_visitor.cpp | 19 +++++--- .../tensor_to_scalar_limit_visitor.cpp | 2 +- tests/LimitVisitorTest.h | 47 ++++++++++++++++++- 5 files changed, 67 insertions(+), 13 deletions(-) diff --git a/include/numsim_cas/core/limit_algebra.h b/include/numsim_cas/core/limit_algebra.h index 393fa3c6..0703034a 100644 --- a/include/numsim_cas/core/limit_algebra.h +++ b/include/numsim_cas/core/limit_algebra.h @@ -15,7 +15,7 @@ class limit_algebra { static limit_result apply_log(limit_result a, bool zero_from_above); static limit_result apply_pow(limit_result base, limit_result exponent, bool zero_from_above); - static limit_result apply_sqrt(limit_result a); + static limit_result apply_sqrt(limit_result a, bool zero_from_above); static limit_result apply_abs(limit_result a); static limit_result apply_reciprocal(limit_result a, bool zero_from_above); static limit_result apply_exp(limit_result a); diff --git a/src/numsim_cas/core/limit_algebra.cpp b/src/numsim_cas/core/limit_algebra.cpp index 542b1c55..fe090d58 100644 --- a/src/numsim_cas/core/limit_algebra.cpp +++ b/src/numsim_cas/core/limit_algebra.cpp @@ -207,8 +207,9 @@ limit_result limit_algebra::apply_pow(limit_result base, limit_result exponent, // Only handle constant/finite exponents for now if (is_finite(exponent.dir) || exponent.dir == dir::zero) { if (exponent.dir == dir::zero) { - // c^0 = 1 for any finite nonzero c; 0^0 and inf^0 are indeterminate - if (base.dir == dir::finite_positive || base.dir == dir::finite_negative) + // c^0 = 1 for finite positive c; for c < 0 the real power c^t is + // undefined off the rationals, and 0^0 / inf^0 are indeterminate + if (base.dir == dir::finite_positive) return {dir::finite_positive}; if (base.dir == dir::zero || is_infinite(base.dir)) return {dir::indeterminate}; @@ -261,12 +262,15 @@ limit_result limit_algebra::apply_pow(limit_result base, limit_result exponent, return {dir::unknown}; } -limit_result limit_algebra::apply_sqrt(limit_result a) { +limit_result limit_algebra::apply_sqrt(limit_result a, bool zero_from_above) { if (a.dir == dir::indeterminate || a.dir == dir::unknown) return a; switch (a.dir) { case dir::zero: + // sqrt is NaN along an approach to zero from below + if (!zero_from_above) + return {dir::unknown}; return {dir::zero}; case dir::finite_positive: return {dir::finite_positive}; diff --git a/src/numsim_cas/scalar/visitors/scalar_limit_visitor.cpp b/src/numsim_cas/scalar/visitors/scalar_limit_visitor.cpp index c3213542..e146180c 100644 --- a/src/numsim_cas/scalar/visitors/scalar_limit_visitor.cpp +++ b/src/numsim_cas/scalar/visitors/scalar_limit_visitor.cpp @@ -19,8 +19,7 @@ namespace { // Nonnegative by construction, so a limit of zero is approached from above. bool is_structurally_nonnegative( expression_holder const &e) { - if (is_same(e) || is_same(e) || - is_same(e)) + if (is_same(e) || is_same(e)) return true; if (is_same(e)) { auto exponent = try_int_constant(e.get().expr_rhs()); @@ -72,6 +71,10 @@ limit_result target_to_limit(limit_target target) { bool scalar_limit_visitor::zero_from_above(expr_holder_t const &expr) const { if (expr == m_limit_var) return m_target.target == limit_target::point::zero_plus; + // sqrt is nonnegative exactly where it is defined, so it reaches zero from + // above only if its operand does. + if (is_same(expr)) + return zero_from_above(expr.get().expr()); return is_positive(expr) || is_nonnegative(expr) || is_structurally_nonnegative(expr); } @@ -168,9 +171,13 @@ void scalar_limit_visitor::operator()(scalar_negative const &v) { void scalar_limit_visitor::operator()(scalar_pow const &v) { // An even integer exponent makes the power nonnegative from either side. auto exponent = try_int_constant(v.expr_rhs()); - bool from_above = - zero_from_above(v.expr_lhs()) || (exponent && *exponent % 2 == 0); - m_result = apply_pow(apply(v.expr_lhs()), apply(v.expr_rhs()), from_above); + bool base_from_above = zero_from_above(v.expr_lhs()); + bool from_above = base_from_above || (exponent && *exponent % 2 == 0); + auto base = apply(v.expr_lhs()); + m_result = apply_pow(base, apply(v.expr_rhs()), from_above); + // a non-integer power of a base that may reach zero from below is NaN + if (base.dir == dir::zero && !base_from_above && !exponent) + m_result = {dir::unknown}; } // ─── Functions ──────────────────────────────────────────────────── @@ -216,7 +223,7 @@ void scalar_limit_visitor::operator()(scalar_log const &v) { } void scalar_limit_visitor::operator()(scalar_sqrt const &v) { - m_result = apply_sqrt(apply(v.expr())); + m_result = apply_sqrt(apply(v.expr()), zero_from_above(v.expr())); } void scalar_limit_visitor::operator()(scalar_abs const &v) { diff --git a/src/numsim_cas/tensor_to_scalar/visitors/tensor_to_scalar_limit_visitor.cpp b/src/numsim_cas/tensor_to_scalar/visitors/tensor_to_scalar_limit_visitor.cpp index a1907331..fa7842fd 100644 --- a/src/numsim_cas/tensor_to_scalar/visitors/tensor_to_scalar_limit_visitor.cpp +++ b/src/numsim_cas/tensor_to_scalar/visitors/tensor_to_scalar_limit_visitor.cpp @@ -220,7 +220,7 @@ void tensor_to_scalar_limit_visitor::operator()(tensor_to_scalar_exp const &v) { void tensor_to_scalar_limit_visitor::operator()( tensor_to_scalar_sqrt const &v) { - m_result = apply_sqrt(apply(v.expr())); + m_result = apply_sqrt(apply(v.expr()), zero_from_above(v.expr())); } // ─── Constants ──────────────────────────────────────────────────── diff --git a/tests/LimitVisitorTest.h b/tests/LimitVisitorTest.h index d363ccfa..862ce9c2 100644 --- a/tests/LimitVisitorTest.h +++ b/tests/LimitVisitorTest.h @@ -294,6 +294,7 @@ struct limit_algebra_probe : limit_algebra { using limit_algebra::apply_log; using limit_algebra::apply_pow; using limit_algebra::apply_reciprocal; + using limit_algebra::apply_sqrt; using limit_algebra::combine_add; using limit_algebra::combine_mul; }; @@ -360,9 +361,11 @@ TEST(LimitAlgebra, LogAndPowOfFiniteValues) { dir::indeterminate); EXPECT_EQ(A::apply_pow({dir::finite_positive}, {dir::zero}, false).dir, dir::finite_positive); - // (-2)^0 = 1 + // c^t with c < 0 has no real limit as t -> 0 EXPECT_EQ(A::apply_pow({dir::finite_negative}, {dir::zero}, false).dir, - dir::finite_positive); + dir::unknown); + EXPECT_EQ(A::apply_sqrt({dir::zero}, false).dir, dir::unknown); + EXPECT_EQ(A::apply_sqrt({dir::zero}, true).dir, dir::zero); } TEST(ScalarLimit, OddFunctionsAtZero) { @@ -446,6 +449,46 @@ TEST(ScalarLimit, InverseTrigOutsideDomainIsUnknown) { EXPECT_EQ(v.apply(acos(x)).dir, dir::finite_positive); } +// sqrt and non-integer powers are NaN along an approach to zero from below. +TEST(ScalarLimit, RootsNeedANonnegativeApproach) { + auto x = make_expression("x"); + scalar_limit_visitor from_below(x, {pt::zero_minus}); + scalar_limit_visitor from_above(x, {pt::zero_plus}); + EXPECT_EQ(from_below.apply(sqrt(x)).dir, dir::unknown); + EXPECT_EQ(from_above.apply(sqrt(x)).dir, dir::zero); + EXPECT_EQ(from_below.apply(sqrt(sin(x))).dir, dir::unknown); + EXPECT_EQ(from_below.apply(pow(x, make_scalar_constant(0.5))).dir, + dir::unknown); + EXPECT_EQ(from_above.apply(pow(x, make_scalar_constant(0.5))).dir, dir::zero); + // an integer power stays defined from either side + EXPECT_EQ(from_below.apply(pow(x, make_scalar_constant(3))).dir, dir::zero); + // sqrt of a nonnegative argument keeps its side + EXPECT_EQ(from_below.apply(sqrt(abs(x))).dir, dir::zero); + EXPECT_EQ(from_above.apply(log(sqrt(x))).dir, dir::neg_infinity); +} + +// sqrt() is nonnegative only where it is defined, so it cannot lend its sign +// to an argument that reaches zero from below. +TEST(ScalarLimit, SignOfRootFollowsTheOperand) { + auto x = make_expression("x"); + scalar_limit_visitor from_below(x, {pt::zero_minus}); + scalar_limit_visitor from_above(x, {pt::zero_plus}); + EXPECT_EQ(from_below.apply(sign(sqrt(x))).dir, dir::unknown); + EXPECT_EQ(from_above.apply(sign(sqrt(x))).dir, dir::finite_positive); + EXPECT_EQ(from_below.apply(sign(sqrt(abs(x)))).dir, dir::finite_positive); +} + +// c^t with c < 0 and t -> 0 has no real limit. +TEST(ScalarLimit, ZeroPowerOfANegativeBaseIsUnknown) { + auto x = make_expression("x"); + auto [n, p] = make_scalar_variable("n", "p"); + n.assumption(negative{}); + p.assumption(positive{}); + scalar_limit_visitor v(x, {pt::zero_plus}); + EXPECT_EQ(v.apply(pow(n, x)).dir, dir::unknown); + EXPECT_EQ(v.apply(pow(p, x)).dir, dir::finite_positive); +} + // A structurally nonnegative argument reaches zero from above, whichever side // the limit variable approaches from. TEST(ScalarLimit, NonnegativeArgumentsKeepTheirZeroSide) {