series_engine.cpp

浏览该文件的文档.

#include <stdexcept>
#include "internal/expression_transform.hpp"
/**
 * @file series_engine.cpp
 */
 
#include "series_engine.hpp"
#include "symbolic.hpp"
#include "internal/symbolic_ast.hpp"
#include "integration.hpp"
#include "internal/expression_analysis.hpp"
#include "internal/normalization_utils.hpp"
#include "internal/series_support.hpp"
#include <cstdint>
#include <memory>
#include <string>
#include <cstdio>
#include <typeinfo>
#include <algorithm>
#include <limits>
#include <optional>
 
 
namespace LMCAS {
 
namespace {
 
 
Result<void> validate_power_series_order(
    int order,
    ComputationContext& context,
    const std::string& operation)
{
    auto step = context.consume_steps(1, operation);
    if (!step) {
        return step;
    }
    if (order <= 0) {
        return Result<void>::failure(
            CasErrc::InvalidArgument,
            "power series truncation order must be positive",
            operation);
    }
    return Result<void>::success();
}
 
 
 
 
} // namespace
namespace detail::series_support {
 
Result<void> validate_series_variable(const std::string& var,
                                      ComputationContext& context,
                                      const std::string& operation)
{
    auto step = context.consume_steps(1, operation);
    if (!step) {
        return step;
    }
    if (var.empty()) {
        return Result<void>::failure(CasErrc::InvalidArgument,
                                     "series variable name cannot be empty",
                                     operation);
    }
    return Result<void>::success();
}
 
std::optional<int> integer_value_from_node(const std::shared_ptr<const SymbolicNode>& node)
{
    BigInt integer;
    if (!try_get_integer_value(std::dynamic_pointer_cast<const NumberNode>(node), integer)) {
        return std::nullopt;
    }
    auto value = integer.try_to_int64();
    if (!value || *value < std::numeric_limits<int>::min() ||
        *value > std::numeric_limits<int>::max()) {
        return std::nullopt;
    }
    return static_cast<int>(*value);
}
 
static std::optional<int> laurent_power_of_power(
    const PowerNode& power, const std::string& var) {
    auto exponent = integer_value_from_node(power.exponent());
    if (!exponent) {
        return std::nullopt;
    }
    auto base = std::dynamic_pointer_cast<const VariableNode>(power.base());
    if (base && !base->is_constant() && base->name() == var) {
        return *exponent;
    }
    auto nested = supported_laurent_integer_power(power.base(), var);
    if (!nested) {
        return std::nullopt;
    }
    const BigInt product = BigInt(*nested) * BigInt(*exponent);
    if (product < BigInt(std::numeric_limits<int>::min()) ||
        product > BigInt(std::numeric_limits<int>::max())) {
        return std::nullopt;
    }
    return static_cast<int>(*product.try_to_int64());
}
 
std::optional<int> supported_laurent_integer_power(
    const std::shared_ptr<const SymbolicNode>& node,
    const std::string& var)
{
    if (!node) {
        return std::nullopt;
    }
    if (std::dynamic_pointer_cast<const NumberNode>(node)) {
        return 0;
    }
    if (auto variable = std::dynamic_pointer_cast<const VariableNode>(node)) {
        return !variable->is_constant() && variable->name() == var
            ? std::optional<int>(1) : std::nullopt;
    }
    if (auto power = std::dynamic_pointer_cast<const PowerNode>(node)) {
        return laurent_power_of_power(*power, var);
    }
    if (auto multiply = std::dynamic_pointer_cast<const MultiplyNode>(node)) {
        BigInt total_power(0);
        for (const auto& operand : multiply->operands()) {
            if (std::dynamic_pointer_cast<const NumberNode>(operand)) {
                continue;
            }
            auto power = supported_laurent_integer_power(operand, var);
            if (!power) {
                return std::nullopt;
            }
            total_power += BigInt(*power);
        }
        if (total_power < BigInt(std::numeric_limits<int>::min()) ||
            total_power > BigInt(std::numeric_limits<int>::max())) {
            return std::nullopt;
        }
        return static_cast<int>(*total_power.try_to_int64());
    }
    return std::nullopt;
 
 
}
 
Result<void> validate_power_series_coefficients(
    const std::vector<std::shared_ptr<SymbolicExpr>>& coeffs,
    const std::string& operation,
    const std::string& name)
{
    for (size_t i = 0; i < coeffs.size(); ++i) {
        if (!coeffs[i] || !LMCAS::detail::node(coeffs[i])) {
            return Result<void>::failure(
                CasErrc::InvalidArgument,
                name + " contains a null coefficient at index " + std::to_string(i),
                operation);
        }
    }
    return Result<void>::success();
}
} // namespace detail::series_support
 
using detail::series_support::validate_power_series_coefficients;
 
std::vector<std::shared_ptr<SymbolicExpr>> power_series_add(
    const std::vector<std::shared_ptr<SymbolicExpr>>& a,
    const std::vector<std::shared_ptr<SymbolicExpr>>& b) {
    ComputationContext context;
    detail::RewriteBudget transform_budget(context, (context).limits().max_recursion_depth,
                                            (context).limits().max_ast_nodes, "transform");
    size_t len = std::max(a.size(), b.size());
    std::vector<std::shared_ptr<SymbolicExpr>> result(len);
    for (size_t i = 0; i < len; ++i) {
        auto ai = (i < a.size() && a[i]) ? a[i] : SymbolicExpr::number(0);
        auto bi = (i < b.size() && b[i]) ? b[i] : SymbolicExpr::number(0);
        result[i] = detail::make_expression_ptr(detail::normalize_node(detail::node(SymbolicExpr::add(ai, bi)), transform_budget));
    }
    return result;
}
 
PowerSeriesResult power_series_multiply_checked(
    const std::vector<std::shared_ptr<SymbolicExpr>>& a,
    const std::vector<std::shared_ptr<SymbolicExpr>>& b, int order,
    ComputationContext& context) try {
    detail::RewriteBudget transform_budget(context, (context).limits().max_recursion_depth,
                                            (context).limits().max_ast_nodes, "transform");
    const std::string operation = "power_series_multiply";
    auto order_check = validate_power_series_order(order, context, operation);
    if (!order_check) {
        return PowerSeriesResult::failure(order_check.error());
    }
    auto a_check = validate_power_series_coefficients(a, operation, "left series");
    if (!a_check) {
        return PowerSeriesResult::failure(a_check.error());
    }
    auto b_check = validate_power_series_coefficients(b, operation, "right series");
    if (!b_check) {
        return PowerSeriesResult::failure(b_check.error());
    }
 
    size_t n = static_cast<size_t>(order);
    std::vector<std::shared_ptr<SymbolicExpr>> result(n, SymbolicExpr::number(0));
    for (size_t k = 0; k < n; ++k) {
        auto step = context.consume_steps(1, operation);
        if (!step) {
            return PowerSeriesResult::failure(step.error());
        }
        std::vector<std::shared_ptr<const SymbolicNode>> terms;
        for (size_t j = 0; j <= k; ++j) {
            if (j >= a.size() || (k-j) >= b.size()) {
                continue;
            }
            auto aj = a[j] ? a[j] : SymbolicExpr::number(0);
            auto bkj = b[k-j] ? b[k-j] : SymbolicExpr::number(0);
            terms.push_back(LMCAS::detail::node(SymbolicExpr::multiply(aj, bkj)));
        }
        if (terms.empty()) {
            result[k] = SymbolicExpr::number(0);
        } else if (terms.size() == 1) {
            result[k] = detail::make_expression_ptr(detail::normalize_node(terms[0], transform_budget));
        } else {
            result[k] = detail::make_expression_ptr(detail::normalize_node(LMCAS::detail::make_node<AddNode>(terms), transform_budget));
        }
    }
    return PowerSeriesResult::success(std::move(result));
} catch (const CasError& error) {
    return PowerSeriesResult::failure(error);
} catch (const std::bad_alloc&) {
    return PowerSeriesResult::failure(CasErrc::ResourceLimit, "allocation failed", "power_series_multiply");
} catch (const std::length_error& error) {
    return PowerSeriesResult::failure(CasErrc::ResourceLimit, error.what(), "power_series_multiply");
} catch (const std::exception& error) {
    return PowerSeriesResult::failure(CasErrc::InternalInvariant, error.what(), "power_series_multiply");
}
 
PowerSeriesResult power_series_multiply_checked(
    const std::vector<std::shared_ptr<SymbolicExpr>>& a,
    const std::vector<std::shared_ptr<SymbolicExpr>>& b, int order) {
    ComputationContext context;
    return power_series_multiply_checked(a, b, order, context);
}
 
 
static Result<void> composition_inner_powers(
    const std::vector<std::shared_ptr<SymbolicExpr>>& g,
    std::vector<std::vector<std::shared_ptr<SymbolicExpr>>>& gp,
    int order, ComputationContext& context, const std::string& operation) {
    const size_t n = gp.size();
    gp[0].resize(n, SymbolicExpr::number(0));
    gp[0][0] = SymbolicExpr::number(1);
    if (n > 1) {
        gp[1].resize(n, SymbolicExpr::number(0));
        for (size_t i = 0; i < std::min(g.size(), n); ++i) {
            gp[1][i] = g[i] ? g[i] : SymbolicExpr::number(0);
        }
    }
    for (size_t k = 2; k < n; ++k) {
        auto step = context.consume_steps(1, operation);
        if (!step) {
            return step;
        }
        auto multiplied = power_series_multiply_checked(gp[k-1], gp[1], order, context);
        if (!multiplied) {
            return Result<void>::failure(multiplied.error());
        }
        gp[k] = std::move(multiplied.value());
    }
    return Result<void>::success();
}
 
PowerSeriesResult power_series_compose_checked(
    const std::vector<std::shared_ptr<SymbolicExpr>>& f,
    const std::vector<std::shared_ptr<SymbolicExpr>>& g, int order,
    ComputationContext& context) try {
    detail::RewriteBudget transform_budget(context, (context).limits().max_recursion_depth,
                                            (context).limits().max_ast_nodes, "transform");
    const std::string operation = "power_series_compose";
    auto order_check = validate_power_series_order(order, context, operation);
    if (!order_check) {
        return PowerSeriesResult::failure(order_check.error());
    }
    auto f_check = validate_power_series_coefficients(f, operation, "outer series");
    if (!f_check) {
        return PowerSeriesResult::failure(f_check.error());
    }
    auto g_check = validate_power_series_coefficients(g, operation, "inner series");
    if (!g_check) {
        return PowerSeriesResult::failure(g_check.error());
    }
 
    size_t n = static_cast<size_t>(order);
    if (!g.empty() && !g[0]->is_zero()) {
        return PowerSeriesResult::failure(
            CasErrc::DomainError,
            "power series composition around zero requires g(0) = 0",
            operation);
    }
    std::vector<std::shared_ptr<SymbolicExpr>> result(n, SymbolicExpr::number(0));
    std::vector<std::vector<std::shared_ptr<SymbolicExpr>>> gp(n);
    auto powers = composition_inner_powers(g, gp, order, context, operation);
    if (!powers) {
        return PowerSeriesResult::failure(powers.error());
    }
    for (size_t k = 0; k < std::min(f.size(), n); ++k) {
        auto step = context.consume_steps(1, operation);
        if (!step) {
            return PowerSeriesResult::failure(step.error());
        }
        auto fk = f[k] ? f[k] : SymbolicExpr::number(0);
        if (fk->is_zero()) {
            continue;
        }
        for (size_t i = 0; i < n && i < gp[k].size(); ++i)
            result[i] = detail::make_expression_ptr(detail::normalize_node(detail::node(SymbolicExpr::add(result[i], SymbolicExpr::multiply(fk, gp[k][i]))), transform_budget));
    }
    return PowerSeriesResult::success(std::move(result));
} catch (const CasError& error) {
    return PowerSeriesResult::failure(error);
} catch (const std::bad_alloc&) {
    return PowerSeriesResult::failure(CasErrc::ResourceLimit, "allocation failed", "power_series_compose");
} catch (const std::length_error& error) {
    return PowerSeriesResult::failure(CasErrc::ResourceLimit, error.what(), "power_series_compose");
} catch (const std::exception& error) {
    return PowerSeriesResult::failure(CasErrc::InternalInvariant, error.what(), "power_series_compose");
}
 
PowerSeriesResult power_series_compose_checked(
    const std::vector<std::shared_ptr<SymbolicExpr>>& f,
    const std::vector<std::shared_ptr<SymbolicExpr>>& g, int order) {
    ComputationContext context;
    return power_series_compose_checked(f, g, order, context);
}
 
 
 
std::shared_ptr<SymbolicExpr> symbolic_sum(
    const std::shared_ptr<SymbolicExpr>& body, const std::string& index,
    const std::shared_ptr<SymbolicExpr>& lower, const std::shared_ptr<SymbolicExpr>& upper) {
    ComputationContext context;
    detail::RewriteBudget transform_budget(context, (context).limits().max_recursion_depth,
                                            (context).limits().max_ast_nodes, "transform");
    if (!body || !lower || !upper) {
        return nullptr;
    }
    BigInt lv, uv;
    if (try_get_integer_value(
            std::dynamic_pointer_cast<const NumberNode>(LMCAS::detail::node(lower)), lv) &&
        try_get_integer_value(
            std::dynamic_pointer_cast<const NumberNode>(LMCAS::detail::node(upper)), uv)) {
        if (uv < lv) {
            return SymbolicExpr::number(0);
        }
        const BigInt span = uv - lv;
        if (span < BigInt(100)) {
            const auto iterations = *span.try_to_uint64() + 1;
            std::vector<std::shared_ptr<const SymbolicNode>> terms;
            BigInt k = lv;
            const BigInt increment(1);
            for (std::uint64_t iteration = 0; iteration < iterations; ++iteration, k += increment) {
                auto val = SymbolicExpr::number(k);
                auto term = detail::make_expression_ptr(detail::substitute_node(detail::node(body), index, detail::node(val), transform_budget));
                if (term) { term = detail::make_expression_ptr(detail::normalize_node(detail::node(term), transform_budget)); terms.push_back(LMCAS::detail::node(term)); }
            }
            if (terms.empty()) {
                return SymbolicExpr::number(0);
            }
            if (terms.size() == 1) {
                return detail::make_expression_ptr(detail::normalize_node(terms[0], transform_budget));
            }
            return detail::make_expression_ptr(detail::normalize_node(LMCAS::detail::make_node<AddNode>(terms), transform_budget));
        }
    }
    auto k_var = SymbolicExpr::variable(index);
    if (body->to_string() == k_var->to_string()) {
        auto u = upper; auto l = lower;
        auto u1 = SymbolicExpr::add(u, SymbolicExpr::number(1));
        auto su = SymbolicExpr::divide(SymbolicExpr::multiply(u, u1), SymbolicExpr::number(2));
        auto lm1 = SymbolicExpr::add(l, SymbolicExpr::number(-1));
        auto sl = SymbolicExpr::divide(SymbolicExpr::multiply(lm1, l), SymbolicExpr::number(2));
        return detail::make_expression_ptr(detail::normalize_node(detail::node(SymbolicExpr::add(su, SymbolicExpr::multiply(SymbolicExpr::number(-1), sl))), transform_budget));
    }
    return LMCAS::detail::make_expression_ptr(LMCAS::detail::make_node<SummationNode>(LMCAS::detail::node(body), index, LMCAS::detail::node(lower), LMCAS::detail::node(upper)));
}
 
std::shared_ptr<SymbolicExpr> symbolic_product(
    const std::shared_ptr<SymbolicExpr>& body, const std::string& index,
    const std::shared_ptr<SymbolicExpr>& lower, const std::shared_ptr<SymbolicExpr>& upper) {
    ComputationContext context;
    detail::RewriteBudget transform_budget(context, (context).limits().max_recursion_depth,
                                            (context).limits().max_ast_nodes, "transform");
    if (!body || !lower || !upper) {
        return nullptr;
    }
    BigInt lv, uv;
    if (try_get_integer_value(
            std::dynamic_pointer_cast<const NumberNode>(LMCAS::detail::node(lower)), lv) &&
        try_get_integer_value(
            std::dynamic_pointer_cast<const NumberNode>(LMCAS::detail::node(upper)), uv)) {
        if (uv < lv) {
            return SymbolicExpr::number(1);
        }
        const BigInt span = uv - lv;
        if (span < BigInt(50)) {
            const auto iterations = *span.try_to_uint64() + 1;
            auto result = SymbolicExpr::number(1);
            BigInt k = lv;
            const BigInt increment(1);
            for (std::uint64_t iteration = 0; iteration < iterations; ++iteration, k += increment) {
                auto val = SymbolicExpr::number(k);
                auto term = detail::make_expression_ptr(detail::substitute_node(detail::node(body), index, detail::node(val), transform_budget));
                if (term) { term = detail::make_expression_ptr(detail::normalize_node(detail::node(term), transform_budget)); result = SymbolicExpr::multiply(result, term); }
            }
            return detail::make_expression_ptr(detail::normalize_node(detail::node(result), transform_budget));
        }
    }
    return LMCAS::detail::make_expression_ptr(LMCAS::detail::make_node<ProductNode>(LMCAS::detail::node(body), index, LMCAS::detail::node(lower), LMCAS::detail::node(upper)));
}
 
 
 
} // namespace LMCAS