Mrvasympt

sympy/core/expr.py

@@ -2468,6 +2468,9 @@ def series(self, x=None, x0=0, n=6, dir="+", logx=None):
                 if (s1 + o).removeO() == s1:
                     o = S.Zero
 
+            # Try asymptotic expansion
+            if s1.removeO() == 0:
+                return self.subs(x, 1/x).aseries(x, n, logx, False).subs(x, 1/x)
             try:
                 return collect(s1, x) + o
             except NotImplementedError:
@@ -2646,6 +2649,101 @@ def _eval_nseries(self, x, n, logx):
                      nseries calls it.""" % self.func)
                      )
 
+    def aseries(self, x, n=6, logx=None, bound=0, hir=False):
+        """
+        This algorithm is directly induced from the limit computational algorithm
+        provided by Gruntz. It majorly uses the mrv and rewrite sub-routines.
+        The overall idea of this algorithm is first to look for the most
+        rapidly varying subexpression w of a given expression f and then expands f
+        in a series in w. Then same thing is recursively done on the leading coefficient
+        till we get constant coefficients.
+
+        Use the ``hir`` parameter to produce hierarchical series. It stops the recursion
+        at an early level and may provide nicer and more useful results.
+
+        If the most rapidly varrying subexpression of a given expression f is f itself,
+        the algorithm tries to find a normalised representation of the mrv set and rewrites f
+        using this normalised representation.
+        Use the ``bound`` parameter to give limit on rewriting coefficients in its normalised form.
+
+        Examples
+        ========
+
+        >>> from sympy.abc import x, y
+        >>> e = sin(1/x + exp(-x)) - sin(1/x)
+        >>> e.aseries(x)
+        (1/(24*x**4) - 1/(2*x**2) + 1 + O(x**(-6), (x, oo)))*exp(-x)
+        >>> e.aseries(x, n=3, hir=True)
+        -exp(-2*x)*sin(1/x)/2 + exp(-x)*cos(1/x) + O(exp(-3*x), (x, oo))
+
+        >>> e = exp(exp(x)/(1 - 1/x))
+        >>> e.aseries(x, bound=3)
+        exp(exp(x)/x**2)*exp(exp(x)/x)*exp(-exp(x) + exp(x)/(1 - 1/x) - exp(x)/x - exp(x)/x**2)*exp(exp(x))
+        >>> e.aseries(x)
+        exp(exp(x)/(1 - 1/x))
+
+
+        References
+        ==========
+
+        .. [1] A New Algorithm for Computing Asymptotic Series - Dominik Gruntz
+        .. [2] Gruntz thesis - p90
+        .. [3] http://en.wikipedia.org/wiki/Asymptotic_expansion
+
+        """
+        from sympy.series.gruntz import mrv, rewrite, mrv_leadterm
+        from sympy.functions import exp, log
+
+        omega, exps = mrv(self, x)
+        if x in omega:
+            return self.subs(x, exp(x)).aseries(x, n, logx, hir).subs(x, log(x))
+        d = C.Dummy('d', positive=True)
+        f, logw = rewrite(exps, omega, x, d)
+
+        if self in omega:
+            # Need to find a canonical representative
+            if bound <= 0:
+                return self
+            a = self.args[0]
+            s = a.aseries(x, bound=bound)
+            s = s.func(*[t.removeO() for t in s.args])
+            rep = exp(s.subs(x, 1/x).as_leading_term(x).subs(x, 1/x))
+            f = exp(self.args[0] - rep.args[0]) / d
+            logw = log(1/rep)
+
+        s = f.series(d, 0, n, logx=logx)
+        # Hierarchical series: break after first recursion
+        if hir:
+            return s.subs(d, exp(logw))
+
+        o = s.getO()
+        s = s.removeO()
+        if s.func is Add:
+            terms = s.args
+        else:
+            terms = [s]
+        def pow_cmp(x, y):
+            return int(x.as_coeff_exponent(d)[1] - y.as_coeff_exponent(d)[1])
+        terms = sorted(terms, cmp=pow_cmp)
+        s = 0
+        gotO = False
+
+        for t in terms:
+            coeff, expo = t.as_coeff_exponent(d)
+            if coeff.has(x):
+                s1 = coeff.aseries(x, n, logx, bound=bound-1)
+                if gotO and s1.getO():
+                    break
+                elif s1.getO():
+                    gotO = True
+                s += (s1 * d**expo)
+            else:
+                s += t
+        if not o or gotO:
+            return s.subs(d, exp(logw))
+        else:
+            return (s + o).subs(d, exp(logw))
+
     def limit(self, x, xlim, dir='+'):
         """ Compute limit x->xlim.
         """

sympy/series/tests/test_aseries.py

@@ -0,0 +1,48 @@
+from sympy import (Symbol, Rational, ln, exp, log, sqrt, E, O, pi, I, oo, sinh,
+    sin, cosh, cos, tanh, coth, asinh, acosh, atanh, acoth, tan, cot, Integer,
+    PoleError, floor, ceiling, asin, symbols, limit, Piecewise, Eq, sign,
+    Derivative)
+from sympy.abc import x, y, z
+
+from sympy.utilities.pytest import raises, XFAIL
+
+def test_simple():
+    e = sin(1/x + exp(-x)) - sin(1/x)
+    assert e.aseries(x) == (1/(24*x**4) - 1/(2*x**2) + 1 + O(x**(-6), (x, oo)))*exp(-x)
+    assert e.series(x, oo) == (1/(24*x**4) - 1/(2*x**2) + 1 + O(x**(-6), (x, oo)))*exp(-x)
+
+    e = exp(exp(x)) * (exp(sin(1/x + 1/exp(exp(x)))) - exp(sin(1/x)))
+    assert e.aseries(x, n=4) == -1/(2*x**3) + 1/x + 1 + O(x**(-4), (x, oo))
+    e = exp(x) * (exp(1/x + exp(-x)) - exp(1/x))
+    assert e.aseries(x, n=4) == 1/(6*x**3) + 1/(2*x**2) + 1/x + 1 + O(x**(-4), (x, oo))
+    e = exp(sin(1/x + exp(-exp(x)))) - exp(sin(1/x))
+    assert e.aseries(x, n=4) == (-1/(2*x**3) + 1/x + 1 + O(x**(-4), (x, oo)))*exp(-exp(x))
+
+    e = exp(exp(x)/(1 - 1/x))
+    assert e.aseries(x) == exp(exp(x)/(1 - 1/x))
+    assert e.aseries(x, bound=3) == exp(exp(x)/x**2)*exp(exp(x)/x)*exp(-exp(x) + exp(x)/(1 - 1/x) - \
+            exp(x)/x - exp(x)/x**2)*exp(exp(x))
+
+    n = Symbol('n', integer=True)
+    e = (sqrt(n)*log(n)**2*exp(sqrt(log(n))*log(log(n))**2*exp(sqrt(log(log(n)))*log(log(log(n)))**3)))/n
+    assert e.aseries(n) == exp(exp(sqrt(log(log(n)))*log(log(log(n)))**3)*sqrt(log(n))*log(log(n))**2)*log(n)**2/sqrt(n)
+
+
+# TODO: add test cases involving special functions
+def test_special():
+    pass
+
+
+def test_hierarchical():
+    e = sin(1/x + exp(-x))
+    assert e.aseries(x, n=3, hir=True) == -exp(-2*x)*sin(1/x)/2 + \
+            exp(-x)*cos(1/x) + sin(1/x) + O(exp(-3*x), (x, oo))
+    e = sin(x) * cos(exp(-x))
+    assert e.aseries(x, hir=True) == exp(-4*x)*sin(x)/24 - \
+            exp(-2*x)*sin(x)/2 + sin(x) + O(exp(-6*x), (x, oo))
+    raises(PoleError, lambda: e.aseries(x))
+
+    a, b = symbols('a b', integer=True)
+    e = exp(1/x + exp(-x**2) * (exp(a*x) - exp(b*x))) - exp(1/x)
+    assert e.aseries(x, n=3, hir=True) == (exp(2*a*x + 1/x)/2 + exp(2*b*x + 1/x)/2 - \
+            exp(a*x + b*x + 1/x))*exp(-2*x**2) + (exp(a*x + 1/x) - exp(b*x + 1/x))*exp(-x**2) + O(exp(-3*x**2), (x, oo))