solver_linear.cpp

浏览该文件的文档.

#include "internal/expression_transform.hpp"
#include "solver.hpp"
#include "internal/symbolic_ast.hpp"
#include "poly_utils.hpp"
#include "internal/expression_analysis.hpp"
#include "internal/exact_matrix.hpp"
#include "assumption_context.hpp"
#include <iostream>
#include <vector>
#include <string>
#include <map>
#include <algorithm>
#include <numeric>
#include <cmath>
#include <set>
#include <queue>
#include <unordered_map>
#include <optional>
#include <limits>
#include "internal/solver_support.hpp"
 
namespace LMCAS {
using namespace solver_detail;
 
 
std::shared_ptr<SymbolicExpr> solver_detail::to_ptr(const SymbolicExpr& expr) {
    return LMCAS::detail::make_expression_ptr(expr);
}
 
static bool get_integer_value(const std::shared_ptr<const SymbolicNode>& node, BigInt& value) {
    auto num = std::dynamic_pointer_cast<const NumberNode>(node);
    if (!num) {
        return false;
    }
    if (std::holds_alternative<BigInt>(num->value())) {
        value = std::get<BigInt>(num->value());
        return true;
    }
    if (std::holds_alternative<Rational>(num->value())) {
        const auto& r = std::get<Rational>(num->value());
        if (!r.is_integer()) {
            return false;
        }
        value = r.to_bigint();
        return true;
    }
    return false;
}
 
static std::shared_ptr<const NumberNode> add_number_nodes(const std::shared_ptr<const NumberNode>& a, const std::shared_ptr<const NumberNode>& b) {
    if (std::holds_alternative<Rational>(a->value()) || std::holds_alternative<Rational>(b->value())) {
        Rational r1 = std::holds_alternative<Rational>(a->value()) ? std::get<Rational>(a->value()) : Rational(std::get<BigInt>(a->value()));
        Rational r2 = std::holds_alternative<Rational>(b->value()) ? std::get<Rational>(b->value()) : Rational(std::get<BigInt>(b->value()));
        return LMCAS::detail::make_node<NumberNode>(r1 + r2);
    }
 
    return LMCAS::detail::make_node<NumberNode>(std::get<BigInt>(a->value()) + std::get<BigInt>(b->value()));
}
 
static std::shared_ptr<const NumberNode> multiply_number_nodes(const std::shared_ptr<const NumberNode>& a, const std::shared_ptr<const NumberNode>& b) {
    if (std::holds_alternative<Rational>(a->value()) || std::holds_alternative<Rational>(b->value())) {
        Rational r1 = std::holds_alternative<Rational>(a->value()) ? std::get<Rational>(a->value()) : Rational(std::get<BigInt>(a->value()));
        Rational r2 = std::holds_alternative<Rational>(b->value()) ? std::get<Rational>(b->value()) : Rational(std::get<BigInt>(b->value()));
        return LMCAS::detail::make_node<NumberNode>(r1 * r2);
    }
 
    return LMCAS::detail::make_node<NumberNode>(std::get<BigInt>(a->value()) * std::get<BigInt>(b->value()));
}
 
struct NodeLess {
    bool operator()(const std::shared_ptr<const SymbolicNode>& a, const std::shared_ptr<const SymbolicNode>& b) const {
        if (!a || !b) {
            return a < b;
        }
        return a->compare(*b) < 0;
    }
};
 
std::shared_ptr<SymbolicExpr> solver_detail::multiply_no_expand(
    const std::shared_ptr<const SymbolicNode>& term,
    const std::vector<std::shared_ptr<const SymbolicNode>>& den_factors
) {
    std::vector<std::shared_ptr<const SymbolicNode>> factors;
    if (auto mul = std::dynamic_pointer_cast<const MultiplyNode>(term)) {
        factors.insert(factors.end(), mul->operands().begin(), mul->operands().end());
    } else if (term) {
        factors.push_back(term);
    }
    factors.insert(factors.end(), den_factors.begin(), den_factors.end());
 
    std::shared_ptr<const NumberNode> const_acc =
        LMCAS::detail::make_node<NumberNode>(BigInt(1));
    std::map<std::shared_ptr<const SymbolicNode>, std::shared_ptr<const NumberNode>, NodeLess> bases;
 
    for (const auto& op : factors) {
        if (!op) {
            continue;
        }
        if (auto num = std::dynamic_pointer_cast<const NumberNode>(op)) {
            const_acc = multiply_number_nodes(const_acc, num);
            continue;
        }
 
        std::shared_ptr<const SymbolicNode> base = op;
        std::shared_ptr<const NumberNode> exp = LMCAS::detail::make_node<NumberNode>(BigInt(1));
        if (auto pow = std::dynamic_pointer_cast<const PowerNode>(op)) {
            base = pow->base();
            if (auto e_num = std::dynamic_pointer_cast<const NumberNode>(pow->exponent())) {
                exp = e_num;
            }
        }
 
        auto it = bases.find(base);
        if (it == bases.end()) {
            bases[base] = exp;
        } else {
            bases[base] = add_number_nodes(it->second, exp);
        }
    }
 
    std::vector<std::shared_ptr<const SymbolicNode>> final_ops;
    if (!const_acc->is_one()) final_ops.push_back(const_acc);
 
    for (const auto& [base, exp] : bases) {
        if (exp->is_zero()) {
            continue;
        }
        if (exp->is_one()) final_ops.push_back(base);
        else final_ops.push_back(LMCAS::detail::make_node<PowerNode>(base, exp));
    }
 
    if (final_ops.empty()) {
        return SymbolicExpr::number(1);
    }
    if (final_ops.size() == 1) {
        return LMCAS::detail::make_expression_ptr(final_ops[0]);
    }
    return LMCAS::detail::make_expression_ptr(LMCAS::detail::make_node<MultiplyNode>(final_ops));
}
 
bool solver_detail::is_polynomial_node(const std::shared_ptr<const SymbolicNode>& node) {
    if (!node) {
        return false;
    }
    if (std::dynamic_pointer_cast<const NumberNode>(node)) {
        return true;
    }
    if (std::dynamic_pointer_cast<const VariableNode>(node)) {
        return true;
    }
 
    if (auto add = std::dynamic_pointer_cast<const AddNode>(node)) {
        for (const auto& op : add->operands()) {
            if (!is_polynomial_node(op)) {
                return false;
            }
        }
        return true;
    }
 
    if (auto mul = std::dynamic_pointer_cast<const MultiplyNode>(node)) {
        for (const auto& op : mul->operands()) {
            if (!is_polynomial_node(op)) {
                return false;
            }
        }
        return true;
    }
 
    if (auto pow = std::dynamic_pointer_cast<const PowerNode>(node)) {
        BigInt exp(0);
        if (!get_integer_value(pow->exponent(), exp)) {
            return false;
        }
        if (exp < 0) {
            return false;
        }
        return is_polynomial_node(pow->base());
    }
 
    return false;
}
 
static bool collect_operand_denominators(
    const std::vector<std::shared_ptr<const SymbolicNode>>& operands,
    std::vector<std::shared_ptr<const SymbolicNode>>& den_factors,
    std::vector<std::shared_ptr<SymbolicExpr>>& den_constraints) {
    for (const auto& operand : operands) {
        if (!collect_denominator_factors(operand, den_factors, den_constraints)) {
            return false;
        }
    }
    return true;
}
 
static bool collect_power_denominators(
    const PowerNode& power,
    std::vector<std::shared_ptr<const SymbolicNode>>& den_factors,
    std::vector<std::shared_ptr<SymbolicExpr>>& den_constraints) {
    BigInt exp(0);
    if (!get_integer_value(power.exponent(), exp)) {
        return false;
    }
 
    if (exp < 0) {
        if (!is_polynomial_node(power.base())) {
            return false;
        }
        BigInt k = -exp;
        if (k == 1) {
            den_factors.push_back(power.base());
        } else {
            den_factors.push_back(SymbolicFactory::create_power(power.base(), SymbolicFactory::create_number(k)));
        }
        den_constraints.push_back(LMCAS::detail::make_expression_ptr(power.base()));
        return true;
    }
 
    if (!collect_denominator_factors(power.base(), den_factors, den_constraints)) {
        return false;
    }
    return true;
}
 
bool solver_detail::collect_denominator_factors(
    const std::shared_ptr<const SymbolicNode>& node,
    std::vector<std::shared_ptr<const SymbolicNode>>& den_factors,
    std::vector<std::shared_ptr<SymbolicExpr>>& den_constraints
) {
    if (!node) {
        return false;
    }
 
    if (std::dynamic_pointer_cast<const NumberNode>(node) || std::dynamic_pointer_cast<const VariableNode>(node)) {
        return true;
    }
 
    if (auto add = std::dynamic_pointer_cast<const AddNode>(node)) {
        return collect_operand_denominators(add->operands(), den_factors, den_constraints);
    }
 
    if (auto mul = std::dynamic_pointer_cast<const MultiplyNode>(node)) {
        return collect_operand_denominators(mul->operands(), den_factors, den_constraints);
    }
 
    if (auto pow = std::dynamic_pointer_cast<const PowerNode>(node)) {
        return collect_power_denominators(*pow, den_factors, den_constraints);
    }
 
    if (auto func = std::dynamic_pointer_cast<const FunctionNode>(node)) {
        return collect_operand_denominators(func->arguments(), den_factors, den_constraints);
    }
 
    return false;
}
 
static std::pair<SymbolicExpr, SymbolicExpr> isolate_linear_coeff(const SymbolicExpr& expr, const std::string& var, ComputationContext& context) {
 
    detail::RewriteBudget transform_budget(context, (context).limits().max_recursion_depth,
                                            (context).limits().max_ast_nodes, "transform");
 
    auto expr_ptr = to_ptr(expr);
    auto A_ptr = detail::make_expression_ptr(detail::differentiate_node(detail::node(expr_ptr), var, transform_budget));
 
    std::vector<std::shared_ptr<const SymbolicNode>> mops;
    mops.push_back(LMCAS::detail::node(A_ptr));
    mops.push_back(SymbolicFactory::create_variable(var));
    auto term_Ax = SymbolicFactory::create_multiply(mops);
 
    std::vector<std::shared_ptr<const SymbolicNode>> nops;
    nops.push_back(SymbolicFactory::create_number(BigInt(-1)));
    nops.push_back(term_Ax);
    auto neg_term = SymbolicFactory::create_multiply(nops);
 
    std::vector<std::shared_ptr<const SymbolicNode>> aops;
    aops.push_back(LMCAS::detail::node(expr));
    aops.push_back(neg_term);
    auto B_node = SymbolicFactory::create_add(aops);
 
    auto B_expr = LMCAS::detail::expression_from_node(B_node);
    auto B_simp = detail::make_expression_ptr(detail::normalize_node(detail::node(to_ptr(B_expr)), transform_budget));
 
    auto A_expr = LMCAS::detail::expression_from_node(LMCAS::detail::node(A_ptr));
    auto A_simp = detail::make_expression_ptr(detail::normalize_node(detail::node(to_ptr(A_expr)), transform_budget));
 
    return {LMCAS::detail::expression_from_node(LMCAS::detail::node(A_simp)),
            LMCAS::detail::expression_from_node(LMCAS::detail::node(B_simp))};
}
 
std::map<std::string, SymbolicExpr> Solver::solve_linear_system(
    const std::vector<SymbolicExpr>& equations_in,
    const std::vector<std::string>& variables)
{
    ComputationContext context;
    return solver_detail::project_unconditional_linear_system(equations_in, variables, context);
}
 
std::map<std::string, SymbolicExpr> solver_detail::project_unconditional_linear_system(
    const std::vector<SymbolicExpr>& equations,
    const std::vector<std::string>& variables, ComputationContext& context) {
    try {
        return solve_linear_system_impl(equations, variables, context);
    } catch (const CasError& error) {
        // 容器接口以空映射表示求解未决；受检算法保留内部错误分类。
        if (error.code == CasErrc::Inconclusive ||
            error.code == CasErrc::UnsupportedExpression) {
            return {};
        }
        throw;
    }
}
 
std::map<std::string, SymbolicExpr> solver_detail::solve_linear_system_impl(
    const std::vector<SymbolicExpr>& equations_in,
    const std::vector<std::string>& variables, ComputationContext& context) {
    detail::RewriteBudget transform_budget(context, context.limits().max_recursion_depth,
                                            context.limits().max_ast_nodes, "solve_linear_system");
 
    size_t num_vars = variables.size();
    size_t num_eqs = equations_in.size();
 
    std::vector<std::vector<std::shared_ptr<SymbolicExpr>>> matrix(num_eqs, std::vector<std::shared_ptr<SymbolicExpr>>(num_vars + 1));
 
    for (size_t i = 0; i < num_eqs; ++i) {
        SymbolicExpr current_eq = equations_in[i];
 
        SymbolicExpr constant_part = current_eq;
 
        for (size_t j = 0; j < num_vars; ++j) {
            auto [coeff, remainder] = isolate_linear_coeff(constant_part, variables[j], context);
            matrix[i][j] = LMCAS::detail::make_expression_ptr(coeff);
 
            constant_part = remainder;
        }
 
        std::vector<std::shared_ptr<const SymbolicNode>> ops;
        ops.push_back(SymbolicFactory::create_number(BigInt(-1)));
        ops.push_back(LMCAS::detail::node(constant_part));
        auto neg_const = SymbolicFactory::create_multiply(ops);
        matrix[i][num_vars] = LMCAS::detail::make_expression_ptr(neg_const);
    }
 
    detail::ExactMatrixData augmented{
        num_eqs, num_vars + 1, {}};
    augmented.entries.reserve(num_eqs * (num_vars + 1));
    for (const auto& row : matrix) {
        augmented.entries.insert(
            augmented.entries.end(), row.begin(), row.end());
    }
 
    auto solved = detail::solve_linear_exact(
        std::move(augmented), num_vars, context,
        "solve_linear_system");
    if (!solved) throw solved.error();
    if (std::holds_alternative<detail::InconsistentLinearSolution>(solved.value())) {
        return {};
    }
 
    std::vector<std::shared_ptr<SymbolicExpr>> values;
    if (auto* unique =
            std::get_if<detail::UniqueLinearSolution>(&solved.value())) {
        values = std::move(unique->values);
    } else {
        auto& parametric =
            std::get<detail::ParametricLinearSolution>(solved.value());
        values = parametric.particular;
        for (std::size_t parameter = 0;
             parameter < parametric.free_columns.size(); ++parameter) {
            auto free_variable = SymbolicExpr::variable(
                variables[parametric.free_columns[parameter]]);
            for (std::size_t column = 0; column < num_vars; ++column) {
                auto term = SymbolicExpr::multiply(
                    parametric.nullspace_basis[parameter][column],
                    free_variable);
                values[column] = detail::make_expression_ptr(detail::normalize_node(detail::node(SymbolicExpr::add(
                    values[column], term)), transform_budget));
            }
        }
    }
 
    std::map<std::string, SymbolicExpr> solution;
    for (std::size_t column = 0; column < num_vars; ++column) {
        solution.insert_or_assign(
            variables[column], *values[column]);
    }
    return solution;
}
 
 
} // namespace LMCAS