Feature/asymptotics

ramanujantools/asymptotics/__init__.py

@@ -0,0 +1,10 @@
+from .series_matrix import SeriesMatrix
+from .growth_rate import GrowthRate
+from .reducer import PrecisionExhaustedError, Reducer
+
+__all__ = [
+    "PrecisionExhaustedError",
+    "GrowthRate",
+    "SeriesMatrix",
+    "Reducer",
+]

ramanujantools/asymptotics/growth_rate.py

@@ -0,0 +1,235 @@
+from __future__ import annotations
+
+import sympy as sp
+from sympy.abc import t, n
+
+from functools import total_ordering
+
+
+@total_ordering
+class GrowthRate:
+    r"""
+    Represents the formal asymptotic growth rate of a solution to a linear difference equation.
+
+    The asymptotic behavior is captured by the Birkhoff-Trjitzinsky formal series representation,
+    which defines a canonical basis solution $E(n)$ as:
+    $$E(n) = (n!)^{d} \cdot \lambda^n \cdot e^{Q(n)} \cdot n^{D} \cdot (\log n)^{m}$$
+
+    This class mathematically isolates the distinct orders of infinity (the exponents and bases)
+    into a structured object. It acts as an element in a Tropical Semiring, where addition (`+`)
+    filters for the strictly dominant growth rate, and multiplication (`*`) algebraically
+    combines the formal exponents.
+
+    By default, an uninitialized `GrowthRate` represents the additive identity (a "zero" or
+    "dead" growth), as $\lambda = 0$ collapses the entire expression to zero.
+
+    Args:
+        factorial_power: The exponent $d$ applied to the factorial term $(n!)^d$.
+            Strictly dominates all other growth components.
+        exp_base: The base $\lambda$ of the primary exponential growth $\lambda^n$.
+            Evaluated by its absolute magnitude. A value of $0$ nullifies the entire solution.
+        sub_exp: The fractional exponent $Q(n)$ applied to $e^{Q(n)}$.
+            Must be strictly sub-linear (e.g., contains fractional powers like $n^{1/p}$).
+        polynomial_degree: The exponent $D$ applied to the polynomial term $n^D$.
+            Can be a fractional rational shift introduced by gauge transformations.
+        log_power: The exponent $m$ applied to the logarithmic term $(\log n)^m$.
+            Typically represents the Jordan block depth of degenerate eigenvalues.
+    """
+
+    def __init__(
+        self,
+        factorial_power: sp.Rational = sp.S.Zero,
+        exp_base: sp.Expr = sp.S.Zero,
+        sub_exp: sp.Expr = sp.S.Zero,
+        polynomial_degree: sp.Expr = sp.S.Zero,
+        log_power: sp.Rational = sp.S.Zero,
+    ):
+        self.factorial_power: sp.Expr = sp.S(factorial_power)
+        self.exp_base: sp.Expr = sp.S(exp_base)
+        self.sub_exp: sp.Expr = sp.S(sub_exp)
+        self.polynomial_degree: sp.Expr = sp.S(polynomial_degree)
+        self.log_power: sp.Expr = sp.S(log_power)
+
+        if not self.factorial_power.is_number:
+            raise ValueError("Factorial power must be an integer.")
+
+        if not self.log_power.is_number:
+            raise ValueError("Factorial power must be an integer.")
+
+        syms = (
+            self.exp_base.free_symbols
+            | self.sub_exp.free_symbols
+            | self.polynomial_degree.free_symbols
+        )
+
+        if not syms.issubset({n}):
+            raise ValueError(
+                "Only 'n' can be a free symbol in the growth rate components."
+            )
+
+    def __add__(self, other: GrowthRate) -> GrowthRate:
+        """Addition acts as a max() filter, keeping only the dominant GrowthRate."""
+        if not isinstance(other, GrowthRate):
+            return NotImplemented
+        return self if self > other else other
+
+    def __radd__(self, other: GrowthRate) -> GrowthRate:
+        return self.__add__(other)
+
+    def __mul__(self, other: GrowthRate) -> GrowthRate:
+        """Strictly combines two GrowthRates by adding their formal exponents."""
+        if not isinstance(other, GrowthRate):
+            return NotImplemented
+
+        return GrowthRate(
+            factorial_power=sp.simplify(self.factorial_power + other.factorial_power),
+            exp_base=sp.simplify(self.exp_base * other.exp_base),
+            sub_exp=sp.simplify(self.sub_exp + other.sub_exp),
+            polynomial_degree=sp.simplify(
+                self.polynomial_degree + other.polynomial_degree
+            ),
+            log_power=self.log_power + other.log_power,
+        )
+
+    def __rmul__(self, other: GrowthRate) -> GrowthRate:
+        return self.__mul__(other)
+
+    def __eq__(self, other: GrowthRate) -> bool:
+        if not isinstance(other, GrowthRate):
+            return NotImplemented
+
+        return (
+            self.factorial_power == other.factorial_power
+            and self.exp_base == other.exp_base
+            and self.sub_exp == other.sub_exp
+            and self.polynomial_degree == other.polynomial_degree
+            and self.log_power == other.log_power
+        )
+
+    def __gt__(self, other: GrowthRate) -> bool:
+        if not isinstance(other, GrowthRate):
+            return NotImplemented
+
+        n_real = sp.Symbol(n.name, real=True, positive=True)
+
+        def is_greater(a, b):
+            diff = sp.simplify(a - b)
+            if diff.is_zero:
+                return None
+
+            diff_real = diff.subs(n, n_real)
+            diff_re = sp.re(diff_real)
+
+            lim = sp.limit(diff_re, n_real, sp.oo)
+
+            if lim == sp.oo or lim.is_positive:
+                return True
+            if lim == -sp.oo or lim.is_negative:
+                return False
+
+            if lim.is_number:
+                try:
+                    val = float(lim.evalf())
+                    if val > 0:
+                        return True
+                    if val < 0:
+                        return False
+                except TypeError:
+                    pass
+
+            return None
+
+        cmp_d = is_greater(self.factorial_power, other.factorial_power)
+        if cmp_d is not None:
+            return cmp_d
+
+        cmp_lam = is_greater(sp.Abs(self.exp_base), sp.Abs(other.exp_base))
+        if cmp_lam is not None:
+            return cmp_lam
+
+        cmp_Q = is_greater(self.sub_exp, other.sub_exp)
+        if cmp_Q is not None:
+            return cmp_Q
+
+        cmp_D = is_greater(self.polynomial_degree, other.polynomial_degree)
+        if cmp_D is not None:
+            return cmp_D
+
+        cmp_log = is_greater(self.log_power, other.log_power)
+        if cmp_log is not None:
+            return cmp_log
+
+        return sp.default_sort_key(self.as_expr(n)) > sp.default_sort_key(
+            other.as_expr(n)
+        )
+
+    def __repr__(self) -> str:
+        return (
+            f"GrowthRate(factorial_power={self.factorial_power}, exp_base={self.exp_base}, "
+            f"sub_exp={self.sub_exp}, polynomial_degree={self.polynomial_degree}, log_power={self.log_power})"
+        )
+
+    def __str__(self) -> str:
+        return str(self.as_expr(sp.Symbol("n")))
+
+    def as_expr(self, n: sp.Symbol) -> sp.Expr:
+        """Renders the formal growth as a SymPy expression."""
+        expr = (
+            (sp.factorial(n) ** self.factorial_power)
+            * (self.exp_base**n if self.exp_base != 0 else 0)
+            * sp.exp(self.sub_exp)
+            * (n**self.polynomial_degree)
+            * sp.log(n) ** self.log_power
+        )
+        return sp.simplify(expr).rewrite(sp.factorial)
+
+    def simplify(self) -> GrowthRate:
+        """Returns a new GrowthRate with all components simplified."""
+        return GrowthRate(
+            factorial_power=sp.simplify(self.factorial_power),
+            exp_base=sp.simplify(self.exp_base),
+            sub_exp=sp.simplify(self.sub_exp),
+            polynomial_degree=sp.simplify(self.polynomial_degree),
+            log_power=self.log_power,
+        )
+
+    @classmethod
+    def from_taylor_coefficients(
+        cls, coeffs: list[sp.Expr], p: int, log_power: int = 0, factorial_power: int = 0
+    ) -> GrowthRate:
+        """
+        Calculates the exact asymptotic bounds of a formal product by extracting
+        the sub-exponential and polynomial degrees via a logarithmic Maclaurin expansion.
+        """
+        exp_base = sp.cancel(sp.expand(coeffs[0]))
+        if exp_base == sp.S.Zero:
+            return cls(exp_base=sp.S.Zero)
+
+        precision = len(coeffs)
+
+        x = sum(
+            (coeffs[k] / exp_base) * (t**k) for k in range(1, min(precision, p + 1))
+        )
+
+        # Maclaurin series of ln(1+x) up to O(t^(p+1))
+        log_series = sp.expand(
+            sum(((-1) ** (j + 1) / sp.Rational(j)) * (x**j) for j in range(1, p + 1))
+        )
+        sub_exp, poly_deg = sp.S.Zero, sp.S.Zero
+
+        for k in range(1, p + 1):
+            c_k = sp.cancel(sp.expand(log_series.coeff(t, k)))
+            if c_k != sp.S.Zero:
+                if k < p:
+                    power = 1 - sp.Rational(k, p)
+                    sub_exp += (c_k / power) * (n**power)
+                else:
+                    poly_deg = c_k
+
+        return cls(
+            exp_base=exp_base,
+            sub_exp=sub_exp,
+            polynomial_degree=poly_deg,
+            log_power=log_power,
+            factorial_power=factorial_power,
+        )

ramanujantools/asymptotics/growth_rate_test.py

@@ -0,0 +1,144 @@
+import pytest
+
+import sympy as sp
+from sympy.abc import n
+
+from ramanujantools.asymptotics.growth_rate import GrowthRate
+
+
+def test_equality_type_gate():
+    growth_rate = GrowthRate()
+    assert growth_rate is not None
+
+    assert growth_rate.__eq__(1) == NotImplemented
+    assert growth_rate.__eq__("Some String") == NotImplemented
+
+
+def test_tropical_dot_product_simulation():
+    """Verify the combined workflow used in the final U_inv * CFM matrix multiplication."""
+    g_base = GrowthRate(polynomial_degree=2, exp_base=5)
+    g_shift = GrowthRate(polynomial_degree=7, exp_base=1)
+
+    # The Tropical Dot Product
+    result = g_base * g_shift
+
+    assert result.polynomial_degree == 9
+    assert result.exp_base == 5
+
+
+def test_simplification():
+    """Equality must survive mathematically identical but unsimplified expressions."""
+    growth_rate = GrowthRate(
+        polynomial_degree=n**2 - n**2,
+        sub_exp=sp.expand((n + 1) ** 2 - n**2 - 2 * n - 1),
+    )
+
+    assert GrowthRate(polynomial_degree=0) == growth_rate.simplify()
+
+
+def test_add_type_gate():
+    g = GrowthRate()
+    assert g.__add__(3) == NotImplemented
+
+
+def test_add_max_filter():
+    """Addition must act as a strict max() gatekeeper using __gt__."""
+    dominant = GrowthRate(factorial_power=2)
+    weak = GrowthRate(factorial_power=1)
+
+    assert dominant + weak == dominant
+    assert weak + dominant == dominant
+
+
+def test_add_zero_passthrough():
+    """Addition with 0 or None must pass the object through untouched."""
+    growth_rate = GrowthRate(factorial_power=2, exp_base=3)
+    assert growth_rate + GrowthRate() == growth_rate
+    assert GrowthRate() + growth_rate == growth_rate
+
+
+def test_mul_type_gate():
+    g = GrowthRate()
+    assert g.__mul__(n**2) == NotImplemented
+
+
+def test_mul_growth_combination():
+    """Multiplication of two GrowthRates must cleanly combine their formal exponents."""
+    g1 = GrowthRate(
+        factorial_power=1,
+        exp_base=2,
+        sub_exp=sp.sqrt(n),
+        polynomial_degree=2,
+        log_power=1,
+    )
+    g2 = GrowthRate(
+        factorial_power=2, exp_base=3, sub_exp=n, polynomial_degree=3, log_power=2
+    )
+
+    result = g1 * g2
+
+    assert result.factorial_power == 3
+    assert result.exp_base == 6
+    assert sp.simplify(result.sub_exp - (sp.sqrt(n) + n)) == 0
+    assert result.polynomial_degree == 5
+    assert result.log_power == 3
+
+
+def test_gt_level_1_factorial():
+    """factorial_power (Factorial) strictly dominates all other bounds."""
+    g1 = GrowthRate(factorial_power=2, exp_base=1, polynomial_degree=-100)
+    g2 = GrowthRate(factorial_power=1, exp_base=1000, polynomial_degree=100)
+
+    assert g1 > g2
+    assert not (g2 > g1)
+
+
+def test_gt_level_2_base_exponential():
+    """exp_base (Base Exp) strictly dominates sub_exp (Fractional Exp) and polynomials."""
+    # g2 has a larger base lambda, so it dominates g1 despite g1's massive fractional sub_exp
+    g1 = GrowthRate(exp_base=1, sub_exp=1000 * sp.sqrt(n))
+    g2 = GrowthRate(exp_base=2, sub_exp=0)
+
+    assert g2 > g1
+
+    # Complex magnitude check: |2i| = 2 > |1|
+    g3 = GrowthRate(exp_base=sp.I * 2)
+    g4 = GrowthRate(exp_base=1)
+    assert g3 > g4
+
+
+def test_gt_level_3_fractional_exponential():
+    """sub_exp (Fractional Exp) dominates polynomials."""
+    # Since lambda is tied at 1, sub_exp triggers
+    g1 = GrowthRate(exp_base=1, sub_exp=sp.sqrt(n), polynomial_degree=0)
+    g2 = GrowthRate(exp_base=1, sub_exp=0, polynomial_degree=1000)
+
+    assert g1 > g2
+
+
+def test_gt_level_4_polynomial():
+    """polynomial_degree (Polynomial) dominates logarithmic Jordan depth."""
+    # Since lambda and sub_exp are tied, polynomial_degree triggers
+    g1 = GrowthRate(polynomial_degree=sp.Rational(3, 2), log_power=0)
+    g2 = GrowthRate(polynomial_degree=1, log_power=10)
+
+    assert g1 > g2
+
+
+def test_gt_level_5_logarithmic():
+    """Jordan depth acts as the final logarithmic tie-breaker."""
+    g1 = GrowthRate(log_power=2)
+    g2 = GrowthRate(log_power=1)
+
+    assert g1 > g2
+
+
+def test_gt_complex_oscillation_fallthrough():
+    """
+    Pure imaginary terms in sub_exp (oscillation) have a real limit of 0.
+    The > operator must recognize the tie and fall through to the next level.
+    """
+    g1 = GrowthRate(sub_exp=sp.I * n, polynomial_degree=2)
+    g2 = GrowthRate(sub_exp=0, polynomial_degree=1)
+
+    assert g1 > g2