def rs_series_reversion(p, x, n, y):
    """
    Reversion of a series.

    ``p`` is a series with ``O(x**n)`` of the form `p = a*x + f(x)`
    where `a` is a number different from 0.

    `f(x) = sum( a\_k*x\_k, k in range(2, n))`

      a_k : Can depend polynomially on other variables, not indicated.
      x : Variable with name x.
      y : Variable with name y.

    Solve `p = y`, that is, given `a*x + f(x) - y = 0`,
    find the solution x = r(y) up to O(y**n)

    Algorithm:

    If `r\_i` is the solution at order i, then:
    `a*r\_i + f(r\_i) - y = O(y**(i + 1))`

    and if r_(i + 1) is the solution at order i + 1, then:
    `a*r\_(i + 1) + f(r\_(i + 1)) - y = O(y**(i + 2))`

    We have, r_(i + 1) = r_i + e, such that,
    `a*e + f(r\_i) = O(y**(i + 2))`
    or `e = -f(r\_i)/a`

    So we use the recursion relation:
    `r\_(i + 1) = r\_i - f(r\_i)/a`
    with the boundary condition: `r\_1 = y`

    Examples
    ========

    >>> from sympy.polys.domains import QQ
    >>> from sympy.polys.rings import ring
    >>> from sympy.polys.ring_series import rs_series_reversion, rs_trunc
    >>> R, x, y, a, b = ring('x, y, a, b', QQ)
    >>> p = x - x**2 - 2*b*x**2 + 2*a*b*x**2
    >>> p1 = rs_series_reversion(p, x, 3, y); p1
    -2*y**2*a*b + 2*y**2*b + y**2 + y
    >>> rs_trunc(p.compose(x, p1), y, 3)
    y
    """
    if rs_is_puiseux(p, x):
        raise NotImplementedError
    R = p.ring
    nx = R.gens.index(x)
    y = R(y)
    ny = R.gens.index(y)
    if _has_constant_term(p, x):
        raise ValueError("p must not contain a constant term in the series "
                         "variable")
    a = _coefficient_t(p, (nx, 1))
    zm = R.zero_monom
    assert zm in a and len(a) == 1
    a = a[zm]
    r = y/a
    for i in range(2, n):
        sp = rs_subs(p, {x: r}, y, i + 1)
        sp = _coefficient_t(sp, (ny, i))*y**i
        r -= sp/a
    return r
dataset/ETHPy150Open sympy/sympy/sympy/polys/ring_series.py/rs_series_reversion/original