alkahest_cas/calculus/
limits.rs

1//! Symbolic limits towards finite points or ±∞ via local expansions (`Series`),
2//! L'Hôpital iterations, and algebraic transforms (V2-16).
3
4use crate::calculus::series::{local_expansion, LocalExpansion};
5use crate::diff::{diff, DiffError};
6use crate::kernel::pool::POS_INFINITY_SYMBOL;
7use crate::kernel::{subs, ExprData, ExprId, ExprPool};
8use crate::poly::{poly_normal, RationalFunction};
9use crate::simplify::{simplify, simplify_expanded};
10use crate::SeriesError;
11use std::collections::HashMap;
12use std::fmt;
13
14/// Approach direction toward `point` (real-axis ordering).
15#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
16pub enum LimitDirection {
17    /// Ordinary two-sided limit.
18    Bidirectional,
19    /// Limits with `var > point` (approach from the right on the usual number line picture).
20    Plus,
21    /// Limits with `var < point`.
22    Minus,
23}
24
25#[derive(Debug)]
26pub enum LimitError {
27    /// Sub-problem rejected by [`mod@crate::calculus::series`].
28    Series(SeriesError),
29    /// Derivative unavailable for L'Hôpital.
30    Diff(DiffError),
31    /// Odd-order pole requires a one-sided direction.
32    NeedsOneSided,
33    /// Maximal L'Hôpital / recursion depth exceeded.
34    DepthExceeded,
35    /// No implemented rule applies (non-comparable growth, oscillation, …).
36    Unsupported,
37}
38
39impl fmt::Display for LimitError {
40    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
41        match self {
42            LimitError::Series(e) => write!(f, "{e}"),
43            LimitError::Diff(e) => write!(f, "{e}"),
44            LimitError::NeedsOneSided => {
45                write!(
46                    f,
47                    "two-sided limit undefined at this pole; pass direction Plus or Minus"
48                )
49            }
50            LimitError::DepthExceeded => write!(f, "limit refinement depth exceeded"),
51            LimitError::Unsupported => write!(f, "limit could not be computed with current rules"),
52        }
53    }
54}
55
56impl std::error::Error for LimitError {
57    fn source(&self) -> Option<&(dyn std::error::Error + 'static)> {
58        match self {
59            LimitError::Series(e) => Some(e),
60            LimitError::Diff(e) => Some(e),
61            _ => None,
62        }
63    }
64}
65
66impl crate::errors::AlkahestError for LimitError {
67    fn code(&self) -> &'static str {
68        match self {
69            LimitError::Series(_) => "E-LIMIT-001",
70            LimitError::Diff(_) => "E-LIMIT-002",
71            LimitError::NeedsOneSided => "E-LIMIT-003",
72            LimitError::DepthExceeded => "E-LIMIT-004",
73            LimitError::Unsupported => "E-LIMIT-005",
74        }
75    }
76
77    fn remediation(&self) -> Option<&'static str> {
78        Some(match self {
79            LimitError::Series(_) => {
80                "increase truncation order indirectly by simplifying the expression, or rewrite using standard limits"
81            }
82            LimitError::Diff(_) => {
83                "ensure primitives have differentiation rules, or simplify before taking the limit"
84            }
85            LimitError::NeedsOneSided => "use LimitDirection::Plus or Minus matching the desired one-sided approach",
86            LimitError::DepthExceeded => {
87                "try manual algebra (quotient form, cancellations) or split into simpler sub-expressions"
88            }
89            LimitError::Unsupported => {
90                "unsupported indeterminate — Gruntz-style comparability beyond this prototype is future work"
91            }
92        })
93    }
94}
95
96impl From<SeriesError> for LimitError {
97    fn from(e: SeriesError) -> Self {
98        LimitError::Series(e)
99    }
100}
101
102impl From<DiffError> for LimitError {
103    fn from(e: DiffError) -> Self {
104        LimitError::Diff(e)
105    }
106}
107
108// ---------------------------------------------------------------------------
109
110/// `limit(expr, var, point, dir)` — see [`LimitDirection`].
111///
112/// `point` may be finite or [`ExprPool::pos_infinity`]. Limits at `-∞` use
113/// `pool.mul(pool.integer(-1), pool.pos_infinity())`.
114pub fn limit(
115    expr: ExprId,
116    var: ExprId,
117    point: ExprId,
118    direction: LimitDirection,
119    pool: &ExprPool,
120) -> Result<ExprId, LimitError> {
121    let r = limit_inner(expr, var, point, direction, pool, 0)?;
122    Ok(simplify(fold_known_reals(simplify(r, pool).value, pool), pool).value)
123}
124
125/// `(g^m)^n ↦ g^{m n}` when `m,n ∈ ℤ`, so substitutions like `(1/t)^k` become `t^{-k}` Laurent heads.
126fn flatten_nested_integer_pow(expr: ExprId, pool: &ExprPool) -> ExprId {
127    match pool.get(expr) {
128        ExprData::Pow { base, exp } => {
129            let base = flatten_nested_integer_pow(base, pool);
130            let exp_fl = flatten_nested_integer_pow(exp, pool);
131            if let (
132                ExprData::Pow {
133                    base: b2,
134                    exp: inner_exp,
135                },
136                ExprData::Integer(outer_e),
137            ) = (pool.get(base), pool.get(exp_fl))
138            {
139                if let ExprData::Integer(inner_e) = pool.get(inner_exp) {
140                    let prod = inner_e.0.clone() * outer_e.0.clone();
141                    return pool.pow(flatten_nested_integer_pow(b2, pool), pool.integer(prod));
142                }
143            }
144            pool.pow(base, exp_fl)
145        }
146        ExprData::Mul(xs) => pool.mul(
147            xs.iter()
148                .map(|x| flatten_nested_integer_pow(*x, pool))
149                .collect(),
150        ),
151        ExprData::Add(xs) => pool.add(
152            xs.iter()
153                .map(|x| flatten_nested_integer_pow(*x, pool))
154                .collect(),
155        ),
156        ExprData::Func { name, args } => {
157            let na: Vec<ExprId> = args
158                .iter()
159                .map(|a| flatten_nested_integer_pow(*a, pool))
160                .collect();
161            pool.func(name.clone(), na)
162        }
163        _ => expr,
164    }
165}
166
167/// After ``x ↦ 1/t``, common forms are ``Mul(numer, denom^{-1})`` with ``Pow(t,-1)``
168/// sprinkled through both.  Clear those poles by multiplying by ``t^k`` until
169/// numerator and denominator describe an honest polynomial quotient in ``t``.
170fn canonical_polynomial_quotient_in_var(expr: ExprId, t: ExprId, pool: &ExprPool) -> ExprId {
171    let (n_raw, d_raw) = numerator_denominator(expr, pool);
172    if d_raw == pool.integer(1_i32) {
173        return expr;
174    }
175    for k in 0_i64..=40 {
176        let tk = pool.pow(t, pool.integer(k));
177        let n = simplify_expanded(pool.mul(vec![tk, n_raw]), pool).value;
178        let d = simplify_expanded(pool.mul(vec![tk, d_raw]), pool).value;
179        let (n, d) = match (poly_normal(n, vec![t], pool), poly_normal(d, vec![t], pool)) {
180            (Ok(nn), Ok(dd)) => (nn, dd),
181            _ => continue,
182        };
183        if let Ok(rf) = RationalFunction::from_symbolic(n, d, vec![t], pool) {
184            let nx = rf.numer.to_expr(pool);
185            let dx = rf.denom.to_expr(pool);
186            return simplify(pool.mul(vec![nx, pool.pow(dx, pool.integer(-1_i32))]), pool).value;
187        }
188    }
189    expr
190}
191
192fn limit_inner(
193    expr: ExprId,
194    var: ExprId,
195    point: ExprId,
196    direction: LimitDirection,
197    pool: &ExprPool,
198    depth: u32,
199) -> Result<ExprId, LimitError> {
200    const MAX_DEPTH: u32 = 48;
201    const SERIES_ORDER: u32 = 32;
202    if depth > MAX_DEPTH {
203        return Err(LimitError::DepthExceeded);
204    }
205
206    if !depends_on(expr, var, pool) {
207        if substitution_is_singular(expr, pool) {
208            return Err(LimitError::Unsupported);
209        }
210        return Ok(fold_known_reals(simplify(expr, pool).value, pool));
211    }
212
213    if let Some(r) = try_special_function_limits(expr, var, point, direction, pool)? {
214        return Ok(r);
215    }
216
217    if is_pos_infinity(point, pool) {
218        let t = pool.symbol("__lt_inf", crate::kernel::Domain::Real);
219        let inv_t = pool.pow(t, pool.integer(-1_i32));
220        let mut m = HashMap::new();
221        m.insert(var, inv_t);
222        let e2 = simplify(
223            canonical_polynomial_quotient_in_var(
224                flatten_nested_integer_pow(subs(expr, &m, pool), pool),
225                t,
226                pool,
227            ),
228            pool,
229        )
230        .value;
231        return limit_inner(
232            e2,
233            t,
234            pool.integer(0_i32),
235            LimitDirection::Plus,
236            pool,
237            depth + 1,
238        );
239    }
240
241    if is_neg_infinity(point, pool) {
242        let t = pool.symbol("__lt_ninf", crate::kernel::Domain::Real);
243        let rep = pool.mul(vec![
244            pool.integer(-1_i32),
245            pool.pow(t, pool.integer(-1_i32)),
246        ]);
247        let mut m = HashMap::new();
248        m.insert(var, rep);
249        let e2 = simplify(
250            canonical_polynomial_quotient_in_var(
251                flatten_nested_integer_pow(subs(expr, &m, pool), pool),
252                t,
253                pool,
254            ),
255            pool,
256        )
257        .value;
258        return limit_inner(
259            e2,
260            t,
261            pool.integer(0_i32),
262            LimitDirection::Plus,
263            pool,
264            depth + 1,
265        );
266    }
267
268    if let Some(r) = try_direct_substitution(expr, var, point, pool) {
269        return Ok(r);
270    }
271
272    if let Some(r) = try_x_log_x_at_zero(expr, var, point, direction, pool, depth)? {
273        return Ok(r);
274    }
275
276    if let Some(r) = try_lhopital(expr, var, point, direction, pool, depth)? {
277        return Ok(r);
278    }
279
280    if let Some(r) = try_expansion_limit(expr, var, point, direction, pool, SERIES_ORDER)? {
281        return Ok(r);
282    }
283
284    Err(LimitError::Unsupported)
285}
286
287fn try_x_log_x_at_zero(
288    expr: ExprId,
289    var: ExprId,
290    point: ExprId,
291    direction: LimitDirection,
292    pool: &ExprPool,
293    depth: u32,
294) -> Result<Option<ExprId>, LimitError> {
295    if direction == LimitDirection::Minus {
296        return Ok(None);
297    }
298    if !matches!(pool.get(point), ExprData::Integer(n) if n.0 == 0) {
299        return Ok(None);
300    }
301    let ExprData::Mul(args) = pool.get(expr) else {
302        return Ok(None);
303    };
304    if args.len() != 2 {
305        return Ok(None);
306    }
307    let (a, b) = (args[0], args[1]);
308    let log_of_var = |u: ExprId| {
309        matches!(
310            pool.get(u),
311            ExprData::Func { name, args: av } if name == "log" && av.len() == 1 && av[0] == var
312        )
313    };
314    let is_var = |u: ExprId| u == var;
315    let ok = (is_var(a) && log_of_var(b)) || (is_var(b) && log_of_var(a));
316    if !ok {
317        return Ok(None);
318    }
319    // L'Hôpital on log(x) / x^{-1}: (1/x) / (-1/x^2) = -x  → 0 as x→0+.
320    let f = pool.func("log", vec![var]);
321    let g = pool.pow(var, pool.integer(-1_i32));
322    let fp = diff(f, var, pool)?.value;
323    let gp = diff(g, var, pool)?.value;
324    let ratio = rational_quotient(fp, gp, pool);
325    Ok(Some(limit_inner(
326        ratio,
327        var,
328        point,
329        LimitDirection::Plus,
330        pool,
331        depth + 1,
332    )?))
333}
334
335fn try_special_function_limits(
336    expr: ExprId,
337    var: ExprId,
338    point: ExprId,
339    direction: LimitDirection,
340    pool: &ExprPool,
341) -> Result<Option<ExprId>, LimitError> {
342    let ExprData::Func { name, args } = pool.get(expr) else {
343        return Ok(None);
344    };
345    if args.len() != 1 || args[0] != var {
346        return Ok(None);
347    }
348    match name.as_str() {
349        "exp" => {
350            if is_pos_infinity(point, pool) {
351                return Ok(Some(pool.pos_infinity()));
352            }
353            if is_neg_infinity(point, pool) {
354                return Ok(Some(pool.integer(0_i32)));
355            }
356            if matches!(pool.get(point), ExprData::Integer(n) if n.0 == 0) {
357                return Ok(Some(pool.integer(1_i32)));
358            }
359        }
360        "log" => {
361            if is_pos_infinity(point, pool) {
362                return Ok(Some(pool.pos_infinity()));
363            }
364            if matches!(pool.get(point), ExprData::Integer(n) if n.0 == 0) {
365                if direction == LimitDirection::Plus {
366                    return Ok(Some(neg_infinity(pool)));
367                }
368                return Err(LimitError::NeedsOneSided);
369            }
370        }
371        _ => {}
372    }
373    Ok(None)
374}
375
376fn neg_infinity(pool: &ExprPool) -> ExprId {
377    pool.mul(vec![pool.integer(-1_i32), pool.pos_infinity()])
378}
379
380fn is_pos_infinity(e: ExprId, pool: &ExprPool) -> bool {
381    matches!(
382        pool.get(e),
383        ExprData::Symbol {
384            name,
385            domain: crate::kernel::Domain::Positive,
386            ..
387        } if name == POS_INFINITY_SYMBOL
388    ) || matches!(
389        pool.get(e),
390        ExprData::Symbol {
391            name,
392            domain: crate::kernel::Domain::Real,
393            ..
394        } if name == POS_INFINITY_SYMBOL
395    )
396}
397
398fn is_neg_infinity(e: ExprId, pool: &ExprPool) -> bool {
399    let ExprData::Mul(args) = pool.get(e) else {
400        return false;
401    };
402    if args.len() != 2 {
403        return false;
404    }
405    let (a, b) = (args[0], args[1]);
406    let m_one = pool.integer(-1_i32);
407    (a == m_one && is_pos_infinity(b, pool)) || (b == m_one && is_pos_infinity(a, pool))
408}
409
410fn depends_on(expr: ExprId, var: ExprId, pool: &ExprPool) -> bool {
411    if expr == var {
412        return true;
413    }
414    match pool.get(expr) {
415        ExprData::Add(xs) | ExprData::Mul(xs) => xs.iter().any(|a| depends_on(*a, var, pool)),
416        ExprData::Pow { base, exp } => depends_on(base, var, pool) || depends_on(exp, var, pool),
417        ExprData::Func { args, .. } => args.iter().any(|a| depends_on(*a, var, pool)),
418        ExprData::Piecewise { branches, default } => {
419            branches
420                .iter()
421                .any(|(c, v)| depends_on(*c, var, pool) || depends_on(*v, var, pool))
422                || depends_on(default, var, pool)
423        }
424        ExprData::Predicate { args, .. } => args.iter().any(|a| depends_on(*a, var, pool)),
425        ExprData::Forall { var: bv, body } | ExprData::Exists { var: bv, body } => {
426            bv != var && depends_on(body, var, pool)
427        }
428        ExprData::BigO(a) => depends_on(a, var, pool),
429        ExprData::Integer(_)
430        | ExprData::Rational(_)
431        | ExprData::Float(_)
432        | ExprData::Symbol { .. } => false,
433    }
434}
435
436fn try_direct_substitution(
437    expr: ExprId,
438    var: ExprId,
439    point: ExprId,
440    pool: &ExprPool,
441) -> Option<ExprId> {
442    if quotient_is_zero_over_zero(expr, var, point, pool) {
443        return None;
444    }
445    let mut m = HashMap::new();
446    m.insert(var, point);
447    let raw = subs(expr, &m, pool);
448    if is_zero_times_pole_indeterminate(raw, pool) {
449        return None;
450    }
451    let sub = fold_known_reals(simplify(raw, pool).value, pool);
452    if depends_on(sub, var, pool) || substitution_is_singular(sub, pool) {
453        None
454    } else {
455        Some(sub)
456    }
457}
458
459/// True when ``expr`` is a product quotient `n/d` with `n,d → 0` at substitution (classic `0/0`).
460fn quotient_is_zero_over_zero(expr: ExprId, var: ExprId, point: ExprId, pool: &ExprPool) -> bool {
461    let (n, d) = numerator_denominator(expr, pool);
462    if d == pool.integer(1_i32) {
463        return false;
464    }
465    let n0 = substitute_fully(n, var, point, pool);
466    let d0 = substitute_fully(d, var, point, pool);
467    is_zero_like(n0, pool) && is_zero_like(d0, pool)
468}
469
470/// `0 · (pole at 0)` style indeterminate — must not simplify to misleading `0`.
471fn is_zero_times_pole_indeterminate(expr: ExprId, pool: &ExprPool) -> bool {
472    let factors: Vec<ExprId> = if matches!(pool.get(expr), ExprData::Mul(_)) {
473        flatten_mul(expr, pool)
474    } else {
475        vec![expr]
476    };
477    let mut any_zero_factor = false;
478    let mut any_pole = false;
479    for f in factors {
480        if substitution_is_singular(f, pool) {
481            any_pole = true;
482        }
483        if matches!(pool.get(f), ExprData::Integer(z) if z.0 == 0) {
484            any_zero_factor = true;
485        }
486        if let ExprData::Func { name, args } = pool.get(f) {
487            if args.len() == 1
488                && matches!(name.as_str(), "sin" | "sinh" | "tan")
489                && matches!(pool.get(args[0]), ExprData::Integer(z) if z.0 == 0)
490            {
491                any_zero_factor = true;
492            }
493        }
494    }
495    any_zero_factor && any_pole
496}
497
498/// `true` after substitution if some sub-expression is ``0^{-n}`` (possibly nested via ``(0^{-1})^e``).
499fn substitution_is_singular(expr: ExprId, pool: &ExprPool) -> bool {
500    match pool.get(expr) {
501        ExprData::Pow { base, exp } => {
502            if let ExprData::Integer(nn) = pool.get(exp) {
503                if nn.0 < 0 {
504                    let b = simplify(base, pool).value;
505                    if matches!(pool.get(b), ExprData::Integer(z) if z.0 == 0) {
506                        return true;
507                    }
508                }
509            }
510            substitution_is_singular(base, pool) || substitution_is_singular(exp, pool)
511        }
512        ExprData::Add(xs) | ExprData::Mul(xs) => {
513            xs.iter().any(|a| substitution_is_singular(*a, pool))
514        }
515        ExprData::Func { args, .. } => args.iter().any(|a| substitution_is_singular(*a, pool)),
516        _ => false,
517    }
518}
519
520fn try_lhopital(
521    expr: ExprId,
522    var: ExprId,
523    point: ExprId,
524    direction: LimitDirection,
525    pool: &ExprPool,
526    depth: u32,
527) -> Result<Option<ExprId>, LimitError> {
528    let (nume, deno) = numerator_denominator(expr, pool);
529    if simplify(nume, pool).value == simplify(deno, pool).value {
530        return Ok(None);
531    }
532    let n0 = substitute_fully(nume, var, point, pool);
533    let d0 = substitute_fully(deno, var, point, pool);
534
535    if !is_zero_like(n0, pool) || !is_zero_like(d0, pool) {
536        return Ok(None);
537    }
538
539    let dn = diff(nume, var, pool)?.value;
540    let dd = diff(deno, var, pool)?.value;
541    if dn == nume && dd == deno {
542        return Ok(None);
543    }
544    let quot = rational_quotient(dn, dd, pool);
545    Ok(Some(limit_inner(
546        quot,
547        var,
548        point,
549        direction,
550        pool,
551        depth + 1,
552    )?))
553}
554
555fn substitute_fully(expr: ExprId, var: ExprId, point: ExprId, pool: &ExprPool) -> ExprId {
556    let mut m = HashMap::new();
557    m.insert(var, point);
558    let s = simplify(subs(expr, &m, pool), pool).value;
559    fold_known_reals(s, pool)
560}
561
562fn rational_quotient(n: ExprId, d: ExprId, pool: &ExprPool) -> ExprId {
563    simplify(pool.mul(vec![n, pool.pow(d, pool.integer(-1_i32))]), pool).value
564}
565
566fn is_zero_like(e: ExprId, pool: &ExprPool) -> bool {
567    let e = simplify(e, pool).value;
568    if matches!(pool.get(e), ExprData::Integer(n) if n.0 == 0) {
569        return true;
570    }
571    if let ExprData::Rational(r) = pool.get(e) {
572        if r.0 == 0 {
573            return true;
574        }
575    }
576    if let ExprData::Func { name, args } = pool.get(e) {
577        if args.len() == 1 && matches!(name.as_str(), "sin" | "tan" | "sinh") {
578            return is_zero_like(args[0], pool);
579        }
580    }
581    false
582}
583
584fn is_one_like(e: ExprId, pool: &ExprPool) -> bool {
585    let e = simplify(e, pool).value;
586    if matches!(pool.get(e), ExprData::Integer(n) if n.0 == 1) {
587        return true;
588    }
589    if let ExprData::Rational(r) = pool.get(e) {
590        return r.0 == 1;
591    }
592    false
593}
594
595/// Constant-fold `sin`, `cos`, `exp`, … after limits (`sin(0) → 0`, `cos(0) → 1`).
596fn fold_known_reals(expr: ExprId, pool: &ExprPool) -> ExprId {
597    let e = simplify(expr, pool).value;
598    match pool.get(e) {
599        ExprData::Add(xs) => {
600            let ys: Vec<ExprId> = xs.iter().map(|x| fold_known_reals(*x, pool)).collect();
601            simplify(pool.add(ys), pool).value
602        }
603        ExprData::Mul(xs) => {
604            let ys: Vec<ExprId> = xs.iter().map(|x| fold_known_reals(*x, pool)).collect();
605            simplify(pool.mul(ys), pool).value
606        }
607        ExprData::Pow { base, exp } => {
608            let b = fold_known_reals(base, pool);
609            let xp = fold_known_reals(exp, pool);
610            if is_one_like(b, pool) {
611                return pool.integer(1_i32);
612            }
613            simplify(pool.pow(b, xp), pool).value
614        }
615        ExprData::Func { name, args } if args.len() == 1 => {
616            let inner = fold_known_reals(args[0], pool);
617            if is_zero_like(inner, pool) {
618                match name.as_str() {
619                    "sin" | "tan" | "sinh" => return pool.integer(0_i32),
620                    "cos" | "cosh" => return pool.integer(1_i32),
621                    "exp" => return pool.integer(1_i32),
622                    _ => {}
623                }
624            }
625            simplify(pool.func(name, vec![inner]), pool).value
626        }
627        ExprData::Func { name, args } => {
628            let ys: Vec<ExprId> = args.iter().map(|x| fold_known_reals(*x, pool)).collect();
629            simplify(pool.func(name, ys), pool).value
630        }
631        _ => e,
632    }
633}
634
635fn flatten_mul(expr: ExprId, pool: &ExprPool) -> Vec<ExprId> {
636    match pool.get(expr) {
637        ExprData::Mul(xs) => xs.iter().flat_map(|a| flatten_mul(*a, pool)).collect(),
638        _ => vec![expr],
639    }
640}
641
642fn numerator_denominator(expr: ExprId, pool: &ExprPool) -> (ExprId, ExprId) {
643    let fac = flatten_mul(expr, pool);
644    let mut nums = Vec::new();
645    let mut dens = Vec::new();
646    for f in fac {
647        match pool.get(f) {
648            ExprData::Pow { base, exp } => {
649                if let ExprData::Integer(n) = pool.get(exp) {
650                    let nn = &n.0;
651                    if *nn == 0 {
652                        nums.push(pool.integer(1_i32));
653                    } else if *nn > 0 {
654                        nums.push(f);
655                    } else {
656                        let m = nn
657                            .clone()
658                            .abs()
659                            .to_u64()
660                            .and_then(|u| u32::try_from(u).ok())
661                            .map(|mag| pool.pow(base, pool.integer(mag as i64)));
662                        if let Some(p) = m {
663                            dens.push(p);
664                        } else {
665                            nums.push(f);
666                        }
667                    }
668                } else {
669                    nums.push(f);
670                }
671            }
672            _ => nums.push(f),
673        }
674    }
675    let n = if nums.is_empty() {
676        pool.integer(1_i32)
677    } else if nums.len() == 1 {
678        nums[0]
679    } else {
680        pool.mul(nums)
681    };
682    let d = if dens.is_empty() {
683        pool.integer(1_i32)
684    } else if dens.len() == 1 {
685        dens[0]
686    } else {
687        pool.mul(dens)
688    };
689    (n, d)
690}
691
692fn try_expansion_limit(
693    expr: ExprId,
694    var: ExprId,
695    point: ExprId,
696    direction: LimitDirection,
697    pool: &ExprPool,
698    order: u32,
699) -> Result<Option<ExprId>, LimitError> {
700    let exp = match local_expansion(expr, var, point, order, pool) {
701        Ok(e) => e,
702        Err(_) => return Ok(None),
703    };
704    expansion_to_limit(exp, pool, direction)
705}
706
707fn expansion_to_limit(
708    exp: LocalExpansion,
709    pool: &ExprPool,
710    direction: LimitDirection,
711) -> Result<Option<ExprId>, LimitError> {
712    let LocalExpansion {
713        valuation,
714        coeffs,
715        h_expr: _,
716    } = exp;
717
718    let mut idx = 0usize;
719    while idx < coeffs.len() && is_zero_like(coeffs[idx], pool) {
720        idx += 1;
721    }
722    if idx >= coeffs.len() {
723        // Truncation hit all zeros — indeterminate within this order.
724        return Ok(None);
725    }
726    let power = valuation + idx as i32;
727    let coeff = coeffs[idx];
728
729    if power > 0 {
730        return Ok(Some(pool.integer(0_i32)));
731    }
732    if power == 0 {
733        return Ok(Some(coeff));
734    }
735
736    // Polar — power < 0
737    let pole_order = (-power) as u32;
738    let sgn_c = structural_sign(coeff, pool).unwrap_or(1);
739    if pole_order % 2 == 0 {
740        return Ok(Some(signed_infinity(pool, sgn_c)));
741    }
742    let Some(hdir) = sign_from_h(direction, power) else {
743        return Err(LimitError::NeedsOneSided);
744    };
745    Ok(Some(signed_infinity(pool, sgn_c * hdir)))
746}
747
748/// For odd pole: sign of `h^power` with `power < 0` as `h → 0` from one side.
749fn sign_from_h(direction: LimitDirection, power: i32) -> Option<i8> {
750    if power >= 0 {
751        return Some(1);
752    }
753    let odd = (-power) % 2 != 0;
754    if !odd {
755        return Some(1);
756    }
757    match direction {
758        LimitDirection::Plus => Some(1),
759        LimitDirection::Minus => Some(-1),
760        LimitDirection::Bidirectional => None,
761    }
762}
763
764fn signed_infinity(pool: &ExprPool, sign: i8) -> ExprId {
765    if sign < 0 {
766        neg_infinity(pool)
767    } else {
768        pool.pos_infinity()
769    }
770}
771
772fn structural_sign(e: ExprId, pool: &ExprPool) -> Option<i8> {
773    match pool.get(e) {
774        ExprData::Integer(n) => {
775            if n.0 > 0 {
776                Some(1)
777            } else if n.0 < 0 {
778                Some(-1)
779            } else {
780                None
781            }
782        }
783        ExprData::Rational(r) => {
784            if r.0 == 0 {
785                None
786            } else if r.0 > 0 {
787                Some(1)
788            } else {
789                Some(-1)
790            }
791        }
792        ExprData::Mul(xs) => {
793            let mut s = 1i8;
794            for a in xs {
795                let sa = structural_sign(a, pool)?;
796                s *= sa;
797            }
798            Some(s)
799        }
800        ExprData::Pow { base: _, exp } if matches!(pool.get(exp), ExprData::Integer(n) if n.0.clone() % 2 == 0) => {
801            Some(1)
802        }
803        _ => None,
804    }
805}
806
807#[cfg(test)]
808mod tests {
809    use super::*;
810    use crate::kernel::Domain;
811
812    #[test]
813    fn limit_sin_over_x_zero() {
814        let p = ExprPool::new();
815        let x = p.symbol("x", Domain::Real);
816        let ex = simplify(
817            p.mul(vec![p.func("sin", vec![x]), p.pow(x, p.integer(-1_i32))]),
818            &p,
819        )
820        .value;
821        let r = limit(ex, x, p.integer(0_i32), LimitDirection::Bidirectional, &p).unwrap();
822        assert_eq!(r, p.integer(1_i32));
823    }
824
825    #[test]
826    fn limit_x_log_x_zero_plus() {
827        let p = ExprPool::new();
828        let x = p.symbol("x", Domain::Real);
829        let ex = simplify(p.mul(vec![x, p.func("log", vec![x])]), &p).value;
830        let r = limit(ex, x, p.integer(0_i32), LimitDirection::Plus, &p).unwrap();
831        assert_eq!(r, p.integer(0_i32));
832    }
833
834    #[test]
835    fn limit_exp_inf() {
836        let p = ExprPool::new();
837        let x = p.symbol("x", Domain::Real);
838        let ex = p.func("exp", vec![x]);
839        let r = limit(ex, x, p.pos_infinity(), LimitDirection::Bidirectional, &p).unwrap();
840        assert_eq!(r, p.pos_infinity());
841    }
842
843    #[test]
844    fn limit_x_squared_at_positive_infinity() {
845        let p = ExprPool::new();
846        let x = p.symbol("x", Domain::Real);
847        let ex = simplify(p.pow(x, p.integer(2_i32)), &p).value;
848        let r = limit(ex, x, p.pos_infinity(), LimitDirection::Bidirectional, &p).unwrap();
849        assert_eq!(r, p.pos_infinity(), "{}", p.display(r));
850    }
851
852    #[test]
853    fn rational_x_over_x_plus_one_after_inf_subst() {
854        let p = ExprPool::new();
855        let t = p.symbol("__lt_inf", Domain::Real);
856        let inv = p.pow(t, p.integer(-1));
857        let ex = p.mul(vec![
858            inv,
859            p.pow(p.add(vec![p.integer(1), inv]), p.integer(-1)),
860        ]);
861        let folded = flatten_nested_integer_pow(ex, &p);
862        let canon = canonical_polynomial_quotient_in_var(folded, t, &p);
863        let r = simplify(canon, &p).value;
864        let mut m = HashMap::new();
865        m.insert(t, p.integer(0));
866        let sub = fold_known_reals(simplify(subs(r, &m, &p), &p).value, &p);
867        assert_eq!(sub, p.integer(1), "canonical={}", p.display(canon));
868    }
869}
alkahest_cas/calculus/limits.rs

alkahest_cas/calculus/
limits.rs
