LMCAS
Lamina Computer Algebra System
载入中...
搜索中...
未找到
transcendental_solver_substitution.cpp
浏览该文件的文档.
1#include <stdexcept>
5#include "poly_utils.hpp"
8#include <vector>
9
10namespace LMCAS {
11namespace {
12
13static bool function_candidate(const FunctionNode& function, const std::string& var) {
14 if (function.arguments().size() != 1 ||
15 !expression_depends_on_variable(function.arguments()[0], var)) {
16 return false;
17 }
18 using Type = FunctionNode::FuncType;
19 return function.type() == Type::Exp || function.type() == Type::Sin ||
20 function.type() == Type::Cos || function.type() == Type::Tan;
21}
22
23static bool power_candidate(const PowerNode& power, const std::string& var) {
24 if (!expression_depends_on_variable(power.base(), var) ||
26 return false;
27 }
28 auto exponent = std::dynamic_pointer_cast<const NumberNode>(power.exponent());
29 BigInt value;
31}
32
33
35 const std::shared_ptr<const SymbolicNode>& node,
36 const std::string& var,
37 std::vector<std::shared_ptr<SymbolicExpr>>& candidates) {
38
39 if (!node) {
40 return;
41 }
42
43 if (auto func = std::dynamic_pointer_cast<const FunctionNode>(node)) {
44 if (function_candidate(*func, var)) {
45 candidates.push_back(LMCAS::detail::make_expression_ptr(node));
46 }
47
48 for (auto& arg : func->arguments()) {
49 collect_transcendental_subexprs(arg, var, candidates);
50 }
51 return;
52 }
53
54 if (auto pow = std::dynamic_pointer_cast<const PowerNode>(node)) {
55 if (power_candidate(*pow, var)) {
56 candidates.push_back(LMCAS::detail::make_expression_ptr(pow->base()));
57 }
58
59 collect_transcendental_subexprs(pow->base(), var, candidates);
60 collect_transcendental_subexprs(pow->exponent(), var, candidates);
61 return;
62 }
63
64 if (auto add = std::dynamic_pointer_cast<const AddNode>(node)) {
65 for (auto& op : add->operands()) {
66 collect_transcendental_subexprs(op, var, candidates);
67 }
68 return;
69 }
70
71 if (auto mul = std::dynamic_pointer_cast<const MultiplyNode>(node)) {
72 for (auto& op : mul->operands()) {
73 collect_transcendental_subexprs(op, var, candidates);
74 }
75 return;
76 }
77}
78
79static void deduplicate_candidates(std::vector<std::shared_ptr<SymbolicExpr>>& candidates) {
80 std::vector<std::shared_ptr<SymbolicExpr>> unique;
81 std::vector<std::string> seen;
82 for (auto& c : candidates) {
83 std::string s = c->to_string();
84 bool found = false;
85 for (auto& existing : seen) {
86 if (existing == s) { found = true; break; }
87 }
88 if (!found) {
89 seen.push_back(s);
90 unique.push_back(c);
91 }
92 }
93 candidates = std::move(unique);
94}
95
96static BigInt product_exp_multiplier(const MultiplyNode& product, const std::string& var) {
97 std::shared_ptr<const NumberNode> number;
98 bool has_variable = false;
99 bool has_other = false;
100 for (const auto& operand : product.operands()) {
101 if (auto numeric = std::dynamic_pointer_cast<const NumberNode>(operand)) {
102 if (!number) {
103 number = numeric;
104 } else {
105 has_other = true;
106 }
107 } else if (auto variable = std::dynamic_pointer_cast<const VariableNode>(operand)) {
108 if (!variable->is_constant() && variable->name() == var) {
109 has_variable = true;
110 } else {
111 has_other = true;
112 }
113 } else {
114 has_other = true;
115 }
116 }
117 BigInt value;
118 if (has_variable && number && !has_other &&
119 try_get_integer_value(number, value) && value > 0) {
120 return value;
121 }
122 return 0;
123}
124
125static BigInt extract_exp_multiplier(const std::shared_ptr<const SymbolicNode>& node, const std::string& var) {
126 auto func = std::dynamic_pointer_cast<const FunctionNode>(node);
127 if (!func || func->type() != FunctionNode::FuncType::Exp || func->arguments().size() != 1) {
128 return 0;
129 }
130
131 auto arg = func->arguments()[0];
132
133 if (auto v = std::dynamic_pointer_cast<const VariableNode>(arg)) {
134 if (!v->is_constant() && v->name() == var) {
135 return 1;
136 }
137 }
138
139 if (auto mul = std::dynamic_pointer_cast<const MultiplyNode>(arg)) {
140 return product_exp_multiplier(*mul, var);
141 }
142
143 return 0;
144}
145
146static std::shared_ptr<const SymbolicNode> rewrite_exp_as_u_power(
147 const std::shared_ptr<const SymbolicNode>& node,
148 const std::string& var,
149 const std::string& u_var) {
150
151 if (!node) {
152 return node;
153 }
154
155 const BigInt k = extract_exp_multiplier(node, var);
156 if (k > 0) {
157 auto u_node = LMCAS::detail::make_node<VariableNode>(u_var);
158 if (k == 1) {
159 return u_node;
160 }
162 }
163
164 if (auto add = std::dynamic_pointer_cast<const AddNode>(node)) {
165 std::vector<std::shared_ptr<const SymbolicNode>> new_ops;
166 new_ops.reserve(add->operands().size());
167 for (auto& op : add->operands()) {
168 new_ops.push_back(rewrite_exp_as_u_power(op, var, u_var));
169 }
170 return SymbolicFactory::create_add(std::move(new_ops));
171 }
172
173 if (auto mul = std::dynamic_pointer_cast<const MultiplyNode>(node)) {
174 std::vector<std::shared_ptr<const SymbolicNode>> new_ops;
175 new_ops.reserve(mul->operands().size());
176 for (auto& op : mul->operands()) {
177 new_ops.push_back(rewrite_exp_as_u_power(op, var, u_var));
178 }
179 return SymbolicFactory::create_multiply(std::move(new_ops));
180 }
181
182 if (auto pow = std::dynamic_pointer_cast<const PowerNode>(node)) {
183 auto new_base = rewrite_exp_as_u_power(pow->base(), var, u_var);
184 auto new_exp = rewrite_exp_as_u_power(pow->exponent(), var, u_var);
185 return SymbolicFactory::create_power(new_base, new_exp);
186 }
187
188 if (auto func = std::dynamic_pointer_cast<const FunctionNode>(node)) {
189 std::vector<std::shared_ptr<const SymbolicNode>> new_args;
190 new_args.reserve(func->arguments().size());
191 for (auto& arg : func->arguments()) {
192 new_args.push_back(rewrite_exp_as_u_power(arg, var, u_var));
193 }
194 return LMCAS::detail::make_node<FunctionNode>(func->type(), std::move(new_args));
195 }
196
197 return node;
198}
199
200static bool has_exp_k_var_terms(const std::shared_ptr<const SymbolicNode>& node, const std::string& var) {
201 if (!node) {
202 return false;
203 }
204
205 if (extract_exp_multiplier(node, var) > 0) {
206 return true;
207 }
208
209 if (auto add = std::dynamic_pointer_cast<const AddNode>(node)) {
210 for (auto& op : add->operands()) {
211 if (has_exp_k_var_terms(op, var)) {
212 return true;
213 }
214 }
215 }
216 if (auto mul = std::dynamic_pointer_cast<const MultiplyNode>(node)) {
217 for (auto& op : mul->operands()) {
218 if (has_exp_k_var_terms(op, var)) {
219 return true;
220 }
221 }
222 }
223 if (auto pow = std::dynamic_pointer_cast<const PowerNode>(node)) {
224 if (has_exp_k_var_terms(pow->base(), var)) {
225 return true;
226 }
227 if (has_exp_k_var_terms(pow->exponent(), var)) {
228 return true;
229 }
230 }
231 if (auto func = std::dynamic_pointer_cast<const FunctionNode>(node)) {
232 for (auto& arg : func->arguments()) {
233 if (has_exp_k_var_terms(arg, var)) {
234 return true;
235 }
236 }
237 }
238
239 return false;
240}
241
243 std::shared_ptr<const SymbolicNode> target;
244 std::string u_name;
245
246 std::shared_ptr<const SymbolicNode> replace(const std::shared_ptr<const SymbolicNode>& node) {
247 if (!node) {
248 return node;
249 }
250
251 if (node->equals(*target)) {
253 }
254
255 if (auto add = std::dynamic_pointer_cast<const AddNode>(node)) {
256 std::vector<std::shared_ptr<const SymbolicNode>> new_ops;
257 new_ops.reserve(add->operands().size());
258 for (auto& op : add->operands()) {
259 new_ops.push_back(replace(op));
260 }
261 return SymbolicFactory::create_add(std::move(new_ops));
262 }
263 if (auto mul = std::dynamic_pointer_cast<const MultiplyNode>(node)) {
264 std::vector<std::shared_ptr<const SymbolicNode>> new_ops;
265 new_ops.reserve(mul->operands().size());
266 for (auto& op : mul->operands()) {
267 new_ops.push_back(replace(op));
268 }
269 return SymbolicFactory::create_multiply(std::move(new_ops));
270 }
271 if (auto pow = std::dynamic_pointer_cast<const PowerNode>(node)) {
272 auto new_base = replace(pow->base());
273 auto new_exp = replace(pow->exponent());
274 return SymbolicFactory::create_power(new_base, new_exp);
275 }
276 if (auto func = std::dynamic_pointer_cast<const FunctionNode>(node)) {
277 std::vector<std::shared_ptr<const SymbolicNode>> new_args;
278 new_args.reserve(func->arguments().size());
279 for (auto& arg : func->arguments()) {
280 new_args.push_back(replace(arg));
281 }
282 return LMCAS::detail::make_node<FunctionNode>(func->type(), std::move(new_args));
283 }
284
285 return node;
286 }
287 };
288
290 const std::shared_ptr<SymbolicExpr>& candidate,
291 const std::shared_ptr<SymbolicExpr>& expr, const std::string& var) {
292 auto variable = std::dynamic_pointer_cast<const VariableNode>(LMCAS::detail::node(candidate));
293 if (!variable || variable->is_constant() || variable->name() != var) {
294 return false;
295 }
296 auto polynomial = symbolic_to_poly<SymbolicPolyCoeff>(expr, var);
297 if (!polynomial) {
298 if (polynomial.error().code == CasErrc::UnsupportedExpression) {
299 return false;
300 }
301 return Result<bool>::failure(polynomial.error());
302 }
303 return polynomial.value().degree() >= 2;
304}
305
307 const std::shared_ptr<SymbolicExpr>& expr,
308 const std::string& var, const std::string& u_var) {
310 return false;
311 }
312 auto polynomial = symbolic_to_poly<SymbolicPolyCoeff>(expr, u_var);
313 if (!polynomial) {
314 if (polynomial.error().code == CasErrc::UnsupportedExpression) {
315 return false;
316 }
317 return Result<bool>::failure(polynomial.error());
318 }
319 return polynomial.value().degree() >= 2;
320}
321
322}
323
325 const std::shared_ptr<SymbolicExpr>& expr,
326 const std::string& var) try {
328 detail::RewriteBudget transform_budget(context, (context).limits().max_recursion_depth,
329 (context).limits().max_ast_nodes, "transform");
330
331 if (!expr || !LMCAS::detail::node(expr)) {
332 return std::optional<SubstitutionResult>{};
333 }
334
336 return std::optional<SubstitutionResult>{};
337 }
338
339 const std::string u_var = "_u_subst";
340
341 std::vector<std::shared_ptr<SymbolicExpr>> candidates;
342 collect_transcendental_subexprs(LMCAS::detail::node(expr), var, candidates);
343 deduplicate_candidates(candidates);
344
345 for (auto& h : candidates) {
346
347 auto redundant = redundant_variable_candidate(h, expr, var);
348 if (!redundant) {
349 return Result<std::optional<SubstitutionResult>>::failure(redundant.error());
350 }
351 if (redundant.value()) {
352 continue;
353 }
354
355
356
357 SubstNodeVisitor visitor;
358 visitor.target = LMCAS::detail::node(h);
359 visitor.u_name = u_var;
360
361 auto substituted_node = visitor.replace(LMCAS::detail::node(expr));
362 auto substituted_expr = LMCAS::detail::make_expression_ptr(substituted_node);
363 substituted_expr = detail::make_expression_ptr(detail::normalize_node(detail::node(substituted_expr), transform_budget));
364
365 auto polynomial = is_substitution_polynomial(substituted_expr, var, u_var);
366 if (!polynomial) {
367 return Result<std::optional<SubstitutionResult>>::failure(polynomial.error());
368 }
369 if (polynomial.value()) {
370 return std::optional<SubstitutionResult>{{h, substituted_expr, u_var}};
371 }
372 }
373
374 if (has_exp_k_var_terms(LMCAS::detail::node(expr), var)) {
375 auto rewritten_node = rewrite_exp_as_u_power(LMCAS::detail::node(expr), var, u_var);
376 auto rewritten_expr = LMCAS::detail::make_expression_ptr(rewritten_node);
377 rewritten_expr = detail::make_expression_ptr(detail::normalize_node(detail::node(rewritten_expr), transform_budget));
378
379 auto polynomial = is_substitution_polynomial(rewritten_expr, var, u_var);
380 if (!polynomial) {
381 return Result<std::optional<SubstitutionResult>>::failure(polynomial.error());
382 }
383 if (polynomial.value()) {
384 auto exp_x = SymbolicExpr::exp(SymbolicExpr::variable(var));
385 return std::optional<SubstitutionResult>{{exp_x, rewritten_expr, u_var}};
386 }
387 }
388
389 return std::optional<SubstitutionResult>{};
390} catch (const CasError& error) {
391 return Result<std::optional<SubstitutionResult>>::failure(error);
392} catch (const std::bad_alloc&) {
393 return Result<std::optional<SubstitutionResult>>::failure(CasErrc::ResourceLimit, "allocation failed", "detect_substitution");
394} catch (const std::length_error& error) {
395 return Result<std::optional<SubstitutionResult>>::failure(CasErrc::ResourceLimit, error.what(), "detect_substitution");
396} catch (const std::exception& error) {
397 return Result<std::optional<SubstitutionResult>>::failure(CasErrc::InternalInvariant, error.what(), "detect_substitution");
398}
399
400}
函数节点,表示数学函数调用(三角函数、对数、特殊函数等)。
乘法节点,表示多个操作数的乘积。
const std::vector< std::shared_ptr< const SymbolicNode > > & operands() const noexcept
幂运算节点,表示 base^exponent。
const std::shared_ptr< const SymbolicNode > & exponent() const noexcept
const std::shared_ptr< const SymbolicNode > & base() const noexcept
static Result failure(CasErrc code, std::string message, std::string operation={})
static std::shared_ptr< const SymbolicNode > create_add(std::vector< std::shared_ptr< const SymbolicNode > > ops, detail::RewriteBudget *budget=nullptr)
创建加法节点,自动扁平化嵌套加法并消除零项.
static std::shared_ptr< const SymbolicNode > create_power(std::shared_ptr< const SymbolicNode > base, std::shared_ptr< const SymbolicNode > exponent)
创建幂运算节点,自动处理指数为 0/1 及底数为 0/1 的情况.
static std::shared_ptr< const SymbolicNode > create_multiply(std::vector< std::shared_ptr< const SymbolicNode > > ops, detail::RewriteBudget *budget=nullptr)
创建乘法节点,自动扁平化嵌套乘法并消除单位元.
static BigInt extract_exp_multiplier(const std::shared_ptr< const SymbolicNode > &node, const std::string &var)
static Result< bool > is_substitution_polynomial(const std::shared_ptr< SymbolicExpr > &expr, const std::string &var, const std::string &u_var)
static bool has_exp_k_var_terms(const std::shared_ptr< const SymbolicNode > &node, const std::string &var)
static void deduplicate_candidates(std::vector< std::shared_ptr< SymbolicExpr > > &candidates)
static Result< bool > redundant_variable_candidate(const std::shared_ptr< SymbolicExpr > &candidate, const std::shared_ptr< SymbolicExpr > &expr, const std::string &var)
static BigInt product_exp_multiplier(const MultiplyNode &product, const std::string &var)
static void collect_transcendental_subexprs(const std::shared_ptr< const SymbolicNode > &node, const std::string &var, std::vector< std::shared_ptr< SymbolicExpr > > &candidates)
static std::shared_ptr< const SymbolicNode > rewrite_exp_as_u_power(const std::shared_ptr< const SymbolicNode > &node, const std::string &var, const std::string &u_var)
static bool function_candidate(const FunctionNode &function, const std::string &var)
std::shared_ptr<::LMCAS::SymbolicExpr > make_expression_ptr(SymbolicNodePtr root)
const SymbolicNodePtr & node(const ::LMCAS::SymbolicExpr &expression) noexcept
std::shared_ptr< const Node > make_node(Args &&... args)
SymbolicNodePtr normalize_node(const SymbolicNodePtr &node, RewriteBudget &budget, const FactsQuery &facts, Domain domain)
ExprResult pow(const ExprPtr &base, const ExprPtr &exponent, ComputationContext &context)
const std::string & variable
__attribute__((visibility("default"))) VectorCalculusExprResult surface_integral_scalar_checked(const std __attribute__((visibility("default"))) VectorCalculusExprResult surface_integral_scalar_checked(const std __attribute__((visibility("default"))) std const VectorField const std::string const std::string & v
const std::vector< Polynomial< BigInt > > const BigInt ComputationContext & context
bool try_get_integer_value(const std::shared_ptr< const NumberNode > &node, BigInt &value)
ExprResult mul(const ExprPtr &lhs, const ExprPtr &rhs, ComputationContext &context)
const VectorField const VectorField & c
bool expression_depends_on_variable(const detail::SymbolicNodePtr &expression, const std::string &variable, detail::RewriteBudget *budget)
ExprResult add(const ExprPtr &lhs, const ExprPtr &rhs, ComputationContext &context)
const ExprPtr RelationOp op
const ExprPtr & exponent
Result< std::optional< SubstitutionResult > > detect_substitution(const std::shared_ptr< SymbolicExpr > &expr, const std::string &var)
const std::string const ExprPtr & value
ExprResult function(const std::string &name, std::vector< ExprPtr > arguments)
多项式转换与精确符号多项式代数。
超越方程求解:三角、指数、对数方程的符号反演。
std::shared_ptr< const SymbolicNode > replace(const std::shared_ptr< const SymbolicNode > &node)
内部 AST 定义。