sim_lib_numbers_cas_diff/implementation/
integrate.rs

1//! Symbolic integration of `CasExpr` trees: the recursive `integrate` rules
2//! covering the builtin operators and power-of-variable cases.
3
4use sim_kernel::{Cx, Error, Result, Symbol, Value};
5use sim_lib_numbers_cas::{CasExpr, simplify_expr};
6use sim_lib_numbers_core::domains;
7
8use super::diff::math;
9
10/// The `integrate-sym` symbol: the symbolic integration entry point.
11pub fn integrate_sym_symbol() -> Symbol {
12    Symbol::new("integrate-sym")
13}
14
15/// Integrate a [`CasExpr`] with respect to `var`, returning a simplified
16/// antiderivative tree (without constant of integration).
17///
18/// Covers the builtin arithmetic operators and power-of-variable cases; inputs
19/// outside the supported set produce an error.
20pub fn integrate_cas(cx: &mut Cx, expr: &CasExpr, var: &Symbol) -> Result<CasExpr> {
21    let integral = match expr {
22        CasExpr::Num(value) => op(
23            math("mul"),
24            vec![CasExpr::Num(value.clone()), CasExpr::Var(var.clone())],
25        ),
26        CasExpr::Var(symbol) if symbol == var => integrate_power_of_var(cx, 1, var)?,
27        CasExpr::Var(symbol) => op(
28            math("mul"),
29            vec![CasExpr::Var(symbol.clone()), CasExpr::Var(var.clone())],
30        ),
31        CasExpr::Op(operator, args) if *operator == math("add") => {
32            op(math("add"), integrate_all(cx, args, var)?)
33        }
34        CasExpr::Op(operator, args) if *operator == math("sub") => {
35            op(math("sub"), integrate_all(cx, args, var)?)
36        }
37        CasExpr::Op(operator, args) if *operator == math("mul") => integrate_mul(cx, args, var)?,
38        CasExpr::Op(operator, args) if *operator == math("pow") => integrate_pow(cx, args, var)?,
39        _ => {
40            return Err(Error::Eval(format!(
41                "{} only supports constants, sums, scalar products, and powers of the integration variable",
42                integrate_sym_symbol()
43            )));
44        }
45    };
46    simplify_expr(cx, integral)
47}
48
49fn integrate_all(cx: &mut Cx, args: &[CasExpr], var: &Symbol) -> Result<Vec<CasExpr>> {
50    args.iter().map(|arg| integrate_cas(cx, arg, var)).collect()
51}
52
53fn integrate_mul(cx: &mut Cx, args: &[CasExpr], var: &Symbol) -> Result<CasExpr> {
54    let mut constant = Vec::new();
55    let mut variable = Vec::new();
56    for arg in args {
57        if depends_on(arg, var) {
58            variable.push(arg.clone());
59        } else {
60            constant.push(arg.clone());
61        }
62    }
63    if variable.len() != 1 {
64        return Err(Error::Eval(format!(
65            "{} only handles products with exactly one variable-dependent factor",
66            integrate_sym_symbol()
67        )));
68    }
69    let mut out = constant;
70    out.push(integrate_cas(cx, &variable[0], var)?);
71    Ok(op(math("mul"), out))
72}
73
74fn integrate_pow(cx: &mut Cx, args: &[CasExpr], var: &Symbol) -> Result<CasExpr> {
75    let [base, exponent] = two_args(args)?;
76    if !matches!(base, CasExpr::Var(symbol) if symbol == var) {
77        return Err(Error::Eval(format!(
78            "{} only supports powers of the integration variable",
79            integrate_sym_symbol()
80        )));
81    }
82    let exponent = literal_i64(cx, exponent)?.ok_or_else(|| {
83        Error::Eval(format!(
84            "{} only supports integer exponents for symbolic powers",
85            integrate_sym_symbol()
86        ))
87    })?;
88    integrate_power_of_var(cx, exponent, var)
89}
90
91fn integrate_power_of_var(cx: &mut Cx, exponent: i64, var: &Symbol) -> Result<CasExpr> {
92    if exponent == -1 {
93        return Ok(op(Symbol::new("ln"), vec![CasExpr::Var(var.clone())]));
94    }
95    let next = exponent.checked_add(1).ok_or_else(|| {
96        Error::Eval(format!(
97            "{} exponent {exponent} overflows when raised for integration",
98            integrate_sym_symbol()
99        ))
100    })?;
101    Ok(op(
102        math("mul"),
103        vec![
104            CasExpr::Num(rational_constant(cx, 1, next)?),
105            op(
106                math("pow"),
107                vec![
108                    CasExpr::Var(var.clone()),
109                    CasExpr::Num(integer_constant(cx, next)?),
110                ],
111            ),
112        ],
113    ))
114}
115
116fn depends_on(expr: &CasExpr, var: &Symbol) -> bool {
117    match expr {
118        CasExpr::Num(_) => false,
119        CasExpr::Var(symbol) => symbol == var,
120        CasExpr::Op(_, args) => args.iter().any(|arg| depends_on(arg, var)),
121    }
122}
123
124fn literal_i64(cx: &mut Cx, expr: &CasExpr) -> Result<Option<i64>> {
125    let CasExpr::Num(value) = expr else {
126        return Ok(None);
127    };
128    let display = value.object().display(cx)?;
129    Ok(display.parse::<i64>().ok())
130}
131
132fn integer_constant(cx: &mut Cx, value: i64) -> Result<Value> {
133    if cx
134        .registry()
135        .number_domain_by_symbol(&domains::i64())
136        .is_some()
137    {
138        return cx
139            .factory()
140            .number_literal(domains::i64(), value.to_string());
141    }
142    cx.factory()
143        .number_literal(domains::f64(), format!("{value}.0"))
144}
145
146fn rational_constant(cx: &mut Cx, num: i64, den: i64) -> Result<Value> {
147    if cx
148        .registry()
149        .number_domain_by_symbol(&domains::rational())
150        .is_some()
151    {
152        return cx
153            .factory()
154            .number_literal(domains::rational(), format!("{num}/{den}"));
155    }
156    let value = num as f64 / den as f64;
157    cx.factory()
158        .number_literal(domains::f64(), value.to_string())
159}
160
161fn two_args(args: &[CasExpr]) -> Result<[&CasExpr; 2]> {
162    let [left, right] = args else {
163        return Err(Error::Eval(format!(
164            "{} expects exactly two operands",
165            math("pow")
166        )));
167    };
168    Ok([left, right])
169}
170
171use super::diff::op;
sim_lib_numbers_cas_diff/implementation/integrate.rs

sim_lib_numbers_cas_diff/implementation/
integrate.rs
