r"""Puiseux Series :mod:`abelfunctions.puiseux`
===========================================

Tools for computing Puiseux series. A necessary component for computing
integral bases and with Riemann surfaces.

Classes
-------

.. autosummary::

    PuiseuxTSeries
    PuiseuxXSeries

Functions
---------

.. autosummary::

    puiseux
    newton_iteration
    newton_iteration_step

References
----------

.. [Duval] D. Duval, "Rational puiseux expansions", Compositio
   Mathematica, vol. 70, no. 2, pp. 119-154, 1989.

.. [Poteaux] A. Poteaux, M. Rybowicz, "Towards a Symbolic-Numeric Method
   to Compute Puiseux Series: The Modular Part", preprint

Examples
--------

Contents
--------

"""
import numpy
import sympy

from abelfunctions.puiseux_series_ring import PuiseuxSeriesRing

from sage.all import xgcd, I
from sage.functions.log import log
from sage.functions.other import ceil, real_part, imag_part
from sage.rings.big_oh import O
from sage.rings.infinity import infinity
from sage.rings.laurent_series_ring import LaurentSeriesRing
from sage.rings.qqbar import QQbar, AlgebraicNumber
from sage.rings.rational_field import QQ
from sympy import Point, Segment


def newton_polygon_exceptional(H):
    r"""Computes the exceptional Newton polygon of `H`."""
    R = H.parent()
    x,y = R.gens()
    d = H(0,y).degree(y)
    return [[(0,0),(d,0)]]

def newton_polygon(H, additional_points=[]):
    r"""Computes the Newton polygon of `H`.

    It's assumed that the first generator of `H` here is the "dependent
    variable". For example, if `H = H(x,y)` and we are aiming to compute
    a `y`-covering of the complex `x`-sphere then each monomial of `H`
    is of the form

    .. math::

        a_{ij} x^j y^i.

    Parameters
    ----------
    H : bivariate polynomial

    Returns
    -------
    list
        Returns a list where each element is a list, representing a side
        of the polygon, which in turn contains tuples representing the
        points on the side.

    Note
    ----
    This is written using Sympy's convex hull algorithm for legacy purposes. It
    can certainly be rewritten to use Sage's Polytope but do so *very
    carefully*! There are a number of subtle things going on here due to the
    fact that boundary points are ignored.

    """
    # because of the way sympy.convex_hull computes the convex hull we
    # need to remove all points of the form (0,j) and (i,0) where j > j0
    # and i > i0, the points on the axes closest to the origin
    R = H.parent()
    x, y = R.gens()
    monomials = H.monomials()
    points = [(monom.degree(y), monom.degree(x)) for monom in monomials]
    support = [Point(pt) for pt in points] + additional_points
    i0 = min(P.x for P in support if P.y == 0)
    j0 = min(P.y for P in support if P.x == 0)
    support = [P for P in support if P.x <= i0 and P.y <= j0]
    convex_hull = sympy.convex_hull(*support)

    # special treatment when the hull is just a point or a segment
    if isinstance(convex_hull, Point):
        P = (convex_hull.x, convex_hull.y)
        return [[P]]
    elif isinstance(convex_hull, Segment):
        P = convex_hull.p1
        convex_hull = generalized_polygon_side(convex_hull)
        support.remove(P)
        support.append(convex_hull.p1)
        sides = [convex_hull]
    else:
        # recursive call with generalized point if a generalized newton
        # polygon is needed.
        sides = convex_hull.sides
        first_side = generalized_polygon_side(sides[0])
        if first_side != sides[0]:
            P = first_side.p1
            return newton_polygon(H, additional_points=[P])

    # convert the sides to lists of points
    polygon = []
    for side in sides:
        polygon_side = [P for P in support if P in side]
        polygon_side = sorted(map(lambda P: (int(P.x),int(P.y)), polygon_side))
        polygon.append(polygon_side)

        # stop the moment we hit the i-axis. despite the filtration at
        # the start of this function we need this condition to prevent
        # returning to the starting point of the newton polygon.
        #
        # (See test_puiseux.TestNewtonPolygon.test_multiple)
        if side.p2.y == 0: break

    return polygon

def generalized_polygon_side(side):
    r"""Returns the generalization of a side on the Newton polygon.

    A generalized Newton polygon is one where every side has slope no
    less than -1.

    Parameters
    ----------
    side : sympy.Segment

    Returns
    -------
    side
    """
    if side.slope < -1:
        p1,p2 = side.points
        p1y = p2.x + p2.y
        side = Segment((0,p1y),p2)
    return side

def bezout(q,m):
    r"""Returns :math:`u,v` such that :math:`uq+mv=1`.

    Parameters
    ----------
    q,m : integer
        Two coprime integers with :math:`q > 0`.

    Returns
    -------
    tuple of integers

    """
    if q == 1:
        return (1,0)
    g,u,v = xgcd(q,-m)
    return (u,v)

def transform_newton_polynomial(H, q, m, l, xi):
    r"""Recenters a Newton polynomial at a given singular term.

    Given the Puiseux data :math:`x=\mu x^q, y=x^m(\beta+y)` this
    function returns the polynomial

    .. math::

        \tilde{H} = H(\xi^v x^q, x^m(\xi^u+y)) / x^l.

    where :math:`uq+mv=1`.

    Parameters
    ----------
    H : polynomial in `x` and `y`
    q, m, l, xi : constants
        See above for the definitions of these parameters.

    Returns
    -------
    polynomial
    """
    R = H.parent()
    x,y = R.gens()

    u,v = bezout(q,m)
    newx = (xi**v)*(x**q)
    newy = (x**m)*(xi**u + y)
    newH = H(newx,newy)

    # divide by x**l
    R = newH.parent()
    x,y = R.gens()
    exponents, coefficients = zip(*(newH.dict().items()))
    exponents = [(e[0] - l, e[1]) for e in exponents]
    newH = R(dict(zip(exponents, coefficients)))
    return newH


def newton_data(H, exceptional=False):
    r"""Determines the "newton data" associated with each side of the polygon.

    For each side :math:`\Delta` of the Newton polygon of `H` we
    associate the data :math:`(q,m,l,`phi)` where

    .. math::

        \Delta: qj + mi = l \\
        \phi_{\Delta}(t) = \sum_{(i,j) \in \Delta} a_{ij} t^{(i-i_0)/q}

    Here, :math:`a_ij x^j y_i` is a term in the polynomial :math:`H` and
    :math:`i_0` is the smallest value of :math:`i` belonging to the
    polygon side :math:`\Delta`.

    Parameters
    ----------
    H : sympy.Poly
        Polynomial in `x` and `y`.

    Returns
    -------
    list
        A list of the tuples :math:`(q,m,l,\phi)`.
    """
    R = H.parent()
    x,y = R.gens()

    if exceptional:
        newton = newton_polygon_exceptional(H)
    else:
        newton = newton_polygon(H)

    # special case when the newton polygon is a single point
    if len(newton[0]) == 1:
        return []

    # for each side determine the corresponding newton data: side slope
    # information and corresponding side characteristic polynomial, phi
    result = []
    for side in newton:
        i0, j0 = side[0]
        i1, j1 = side[1]
        slope = QQ(j1 - j0) / QQ(i1 - i0)
        q = slope.denom()
        m = -slope.numer()
        l = min(q*j0 + m*i0, q*j1 + m*i1)
        phi = sum(H.coefficient({y:i, x:j})*x**((i - i0)//int(q)) for i, j in side)
        phi = phi.univariate_polynomial()
        result.append((q, m, l, phi))
    return result


def newton_iteration(G, n):
    r"""Returns a truncated series `y = y(x)` satisfying

    .. math::

        G(x,y(x)) \equiv 0 \bmod{x^r}

    where $r = \ceil{\log_2{n}}$. Based on the algorithm in [XXX].

    Parameters
    ----------
    G, x, y : polynomial
        A polynomial in `x` and `y`.
    n : int
        Requested degree of the series expansion.

    Notes
    -----
    This algorithm returns the series up to order :math:`2^r > n`. Any
    choice of order below :math:`2^r` will return the same series.

    """
    R = G.parent()
    x,y = R.gens()
    if n < 0:
        raise ValueError('Number of terms must be positive. (n=%d'%n)
    elif n == 0:
        return R(0)

    phi = G
    phiprime = phi.derivative(y)
    try:
        pi = R(x).polynomial(x)
        gi = R(0)
        si = R(phiprime(x,gi)).polynomial(x).inverse_mod(pi)
    except NotImplementedError:
        raise ValueError('Newton iteration for computing regular part of '
                         'Puiseux expansion failed. Curve is most likely '
                         'not regular at center.')

    r = ceil(log(n,2))
    for i in range(r):
        gi,si,pi = newton_iteration_step(phi,phiprime,gi,si,pi)
    return R(gi)


def newton_iteration_step(phi, phiprime, g, s, p):
    r"""Perform a single step of the newton iteration algorithm.

    Parameters
    ----------
    phi, phiprime : sympy.Poly
        Equation and its `y`-derivative.
    g, s : sympy.Poly
        Current solution and inverse (conjugate) modulo `p`.
    p : sympy.Poly
        The current modulus. That is, `g` is the Taylor series solution
        to `phi(t,g) = 0` modulo `p`.
    x,y : sympy.Symbol
        Dependent and independent variables, respectively.

    Returns
    -------
    gnext,snext,pnext

    """
    R = phi.parent()
    x,y = R.gens()
    g = R(g).univariate_polynomial()
    s = R(s).univariate_polynomial()
    p = R(p).univariate_polynomial()

    pnext = p**2
    gnext = g - phi(x,g).univariate_polynomial()*s
    gnext = gnext % pnext
    snext = 2*s - phiprime(x,gnext).univariate_polynomial()*s**2
    snext = snext % pnext

    gnext = R(gnext)
    snext = R(snext)
    pnext = R(pnext)
    return gnext,snext,pnext


def puiseux_rational(H, recurse=False):
    r"""Puiseux data for the curve :math:`H` above :math:`(x,y)=(0,0)`.

    Given a polynomial :math:`H = H(x,y)` :func:`puiseux_rational`
    returns the singular parts of all of the Puiseux series centered at
    :math:`x=0, y=0`.

    Parameters
    ----------
    H : polynomial
        A plane curve in `x` and `y`.
    recurse : boolean
        (Default: `True`) A flag used internally to keep track of which
        term in the singular expansion is being computed.

    Returns
    -------
    list of `(G,P,Q)`
        List of tuples where `P` and `Q` are the x- and y-parts of the
        Puiseux series, respectively, and `G` is a polynomial used in
        :func:`newton_iteration` to generate additional terms in the
        y-series.
    """
    R = H.parent()
    x,y = R.gens()

    # when recurse is true, return if the leading order of H(0,y) is y
    if recurse:
        IH = H(0,y).polynomial(y).ord()
        if IH == 1:
            return [(H,x,y)]

    # for each newton polygon side branch out a new puiseux series
    data = newton_data(H, exceptional=(not recurse))
    singular_terms = []
    for q,m,l,phi in data:
        u,v = bezout(q,m)
        for psi,k in phi.squarefree_decomposition():
            roots = psi.roots(ring=QQbar, multiplicities=False)
            map(lambda x: x.exactify(), roots)
            for xi in roots:
                Hprime = transform_newton_polynomial(H, q, m, l, xi)
                next_terms = puiseux_rational(Hprime, recurse=True)
                for (G,P,Q) in next_terms:
                    singular_term = (G, xi**v*P**q, P**m*(xi**u + Q))
                    singular_terms.append(singular_term)

    return singular_terms


def almost_monicize(f):
    r"""Transform `f` to an "almost monic" polynomial.

    Perform a sequence of substitutions of the form

    .. math::

        f \mapsto x^d f(x,y/x)

    such that :math:`l(0) \neq 0` where :math:`l=l(x)` is the leading
    order coefficient of :math:`f`.

    Parameters
    ----------
    f,x,y : sympy.Expr
        An algebraic curve in `x` and `y`.

    Returns
    -------
    g, transform
        A new, almost monic polynomial `g` and a polynomial `transform`
        such that `y -> y/transform`.
    """
    R = f.parent()
    x,y = R.gens()
    transform = R(1)
    monic = False
    while not monic:
        if f.polynomial(y).leading_coefficient()(0) == 0:
            # the denominator is always of the form x**d. Sage, however, has
            # trouble reducing the expression to simplest terms. the following
            # is a manual version
            r = f(x,y/x)
            n = r.numerator().polynomial(x)
            d = r.denominator().degree(x)
            shift = min(n.exponents() + [d])
            n = n.shift(-shift)
            f = R(n(x,y)) # XXX numerator evaluation is important!
            transform *= x
        else:
            monic = True
    return f, transform

def puiseux(f, alpha, beta=None, order=None, parametric=True):
    r"""Singular parts of the Puiseux series above :math:`x=\alpha`.

    Parameters
    ----------
    f : polynomial
        A plane algebraic curve in `x` and `y`.
    alpha : complex
        The x-point over which to compute the Puiseux series of `f`.
    t : variable
        Variable used in the Puiseux series expansions.
    beta : complex
        (Optional) The y-point at which to compute the Puiseux series.
    order : int
        (Default: `None`) If provided, returns Puiseux series expansions
        up the the specified order.

    Returns
    -------
    list of PuiseuxTSeries

    """
    R = f.parent()
    x,y = R.gens()

    # recenter the curve at x=alpha
    if alpha in [infinity,'oo']:
        alpha = infinity
        d = f.degree(x)
        F = f(1/x,y)*x**d
        n,d = F.numerator(), F.denominator()
        falpha,_ = n.polynomial(x).quo_rem(d.univariate_polynomial())
        falpha = falpha(x).numerator()
    else:
        if not isinstance(alpha, AlgebraicNumber):
            alpha = QQbar(QQ(real_part(alpha))+I*QQ(imag_part(alpha)))
        falpha = f(x+alpha,y)

    # determine the points on the curve lying above x=alpha
    R = falpha.parent()
    x,y = R.gens()
    g, transform = almost_monicize(falpha)
    galpha = R(g(0,y)).univariate_polynomial()
    betas = galpha.roots(ring=QQbar, multiplicities=False)

    # filter for requested value of beta. raise error if not found
    if not beta is None:
        betas = [b for b in betas if b == beta]
        if not betas:
            raise ValueError('The point ({0}, {1}) is not on the '
                             'curve {2}.'.format(alpha, beta, f))

    # for each (alpha, beta) determine the corresponding singular parts of the
    # Puiseux series expansions. note that there may be multiple, distinct
    # places above the same point.
    singular_parts = []
    for beta in betas:
        H = g(x,y+beta)
        singular_part_ab = puiseux_rational(H)

        # recenter the result back to (alpha, beta) from (0,0)
        for G,P,Q in singular_part_ab:
            Q += beta
            Q = Q/transform.univariate_polynomial()(P)
            if alpha == infinity:
                P = 1/P
            else:
                P += alpha

            # append to list of singular data
            singular_parts.append((G,P,Q))

    # instantiate PuiseuxTSeries from the singular data
    series = [PuiseuxTSeries(f, alpha, singular_data, order=order)
              for singular_data in singular_parts]
    return series


class PuiseuxTSeries(object):
    r"""A Puiseux t-series about some place :math:`(\alpha, \beta) \in X`.

    A parametric Puiseux series :math:`P(t)` centered at :math:`(x,y) =
    (\alpha, \beta)` is given in terms of a pair of functions

    .. math::

        x(t) = \alpha + \lambda t^e, \\
        y(t) = \sum_{h=0}^\infty \alpha_h t^{n_h},

    where :math:`x(0) = \alpha, y(0) = \beta`.

    The primary reference for the notation and computational method of
    these Puiseux series is D. Duval.


    Attributes
    ----------
    f, x, y : polynomial
    x0 : complex
        The x-center of the Puiseux series expansion.
    ramification_index : rational
        The ramification index :math:`e`.
    terms : list
        A list of exponent-coefficient pairs representing the y-series.
    order : int
        The order of the Puiseux series expansion.

    Methods
    -------
    xseries
    extend
    eval_x
    eval_y

    """
    @property
    def xdata(self):
        return (self.center, self.xcoefficient, self.ramification_index)
    @xdata.setter
    def xdata(self, value):
        self.center, self.xcoefficient, self.ramification_index = value

    @property
    def is_symbolic(self):
        return self._is_symbolic
    @property
    def is_numerical(self):
        return not self._is_symbolic

    @property
    def terms(self):
        terms = self.ypart.laurent_polynomial().dict().items()
        # note that the following greatly affects singularities() and Int()
        if not terms:
            terms = [(0,0)]
        return terms
    @property
    def xdatan(self):
        if self.is_numerical:
            return self.xdata
        else:
            return (complex(self.center),
                    complex(self.xcoefficient),
                    int(self.ramification_index))
    @property
    def order(self):
        return self._singular_order + self._regular_order

    @property
    def nterms(self):
        """Returns the number of non-zero computed terms.

        Parameters
        ----------
        None

        Returns
        -------
        int

        """
        terms = self.ypart.laurent_polynomial().dict().items()
        return len(terms)


    def __init__(self, f, x0, singular_data, order=None):
        r"""Initialize a PuiseuxTSeries using a set of :math:`\pi = \{\tau\}`
        data.

        Parameters
        ----------
        f, x, y : polynomial
            A plane algebraic curve.
        x0 : complex
            The x-center of the Puiseux series expansion.
        singular_data : list
            The output of :func:`singular`.
        t : variable
            The variable in which the Puiseux t series is represented.

        """
        R = f.parent()
        x,y = R.gens()
        extension_polynomial, xpart, ypart = singular_data
        L = LaurentSeriesRing(ypart.base_ring(), 't')
        t = L.gen()

        self.f = f
        self.t = t
        self._xpart = xpart
        self._ypart = ypart

        # store x-part attributes. handle the centered at infinity case
        self.x0 = x0
        if x0 == infinity:
            x0 = QQ(0)
        self.center = x0

        # extract and store information about the x-part of the puiseux series
        xpart = xpart(t,0)
        xpartshift = xpart - x0
        ramification_index, xcoefficient = xpartshift.laurent_polynomial().dict().popitem()
        self.xcoefficient = xcoefficient
        self.ramification_index = QQ(ramification_index).numerator()
        self.xpart = xpart

        # extract and store information about the y-part of the puiseux series
        self.ypart = L(ypart(t,0))
        self._initialize_extension(extension_polynomial)

        # determine the initial order. See the order property
        val = L(ypart(t,O(t))).prec()
        self._singular_order = 0 if val == infinity else val
        self._regular_order = self._p.degree(x)

        # extend to have at least two elements
        self.extend(nterms=1)

        # the curve, x-part, and terms output by puiseux make the puiseux
        # series unique. any mutability only adds terms
        self.__parent = self.ypart.parent()
        self._hash = hash((self.f, self.xpart, self.ypart))

    def parent(self):
        return self.__parent

    def _initialize_extension(self, extension_polynomial):
        r"""Set up regular part extension machinery.

        RootOfs in expressions are not preserved under this
        transformation. (that is, actual algebraic representations are
        calculated.) each RootOf is temporarily replaced by a dummy
        variable

        Parameters
        ----------
        extension_polynomial, x, y : polynomial

        Returns
        -------
        None : None
            Internally sets hidden regular extension attributes.
        """
        R = extension_polynomial.parent()
        x,y = R.gens()

        # store attributes
        _phi = extension_polynomial
        _p = R(x)
        _g = R(0)
        self._phi = _phi
        self._phiprime = _phi.derivative(y)
        self._p = _p
        self._g = _g

        # compute inverse of phi'(g) modulo x and store
        _g = _g.univariate_polynomial()
        _p = _p.univariate_polynomial()
        ppg = self._phiprime.subs({y:_g}).univariate_polynomial()
        _s = ppg.inverse_mod(_p)
        self._s = _s

    def __repr__(self):
        """Print the x- and y-parts of the Puiseux series."""
        s = '('
        s += str(self.xpart)
        s += ', '
        s += str(self.ypart)
        s += ' + O(%s^%s))'%(self.t,self.order)
        return s

    def __hash__(self):
        return self._hash

    def __eq__(self, other):
        r"""Check equality.

        A `PuiseuxTSeries` is uniquely identified by the curve it's
        defined on, its center, x-part terms, and the singular terms of
        the y-part.

        Parameters
        ----------
        other : PuiseuxTSeries

        Returns
        -------
        boolean

        """
        if isinstance(other, PuiseuxTSeries):
            if self._hash == other._hash:
                return True
        return False

    def xseries(self, all_conjugates=True):
        r"""Returns the corresponding x-series.

        Parameters
        ----------
        all_conjugates : bool
            (default: True) If ``True``, returns all conjugates
            x-representations of this Puiseux t-series. If ``False``,
            only returns one representative.

        Returns
        -------
        list
            List of PuiseuxXSeries representations of this PuiseuxTSeries.

        """
        # obtain relevant rings:
        #   o R = parent ring of curve
        #   o L = parent ring of T-series
        #   o S = temporary polynomial ring over base ring of T-series
        #   o P = Puiseux series ring
        L = self.ypart.parent()
        t = L.gen()
        S = L.base_ring()['z']
        z = S.gen()

        R = self.f.parent()
        x,y = R.gens()
        P = PuiseuxSeriesRing(L.base_ring(), str(x))
        x = P.gen()

        # given x = alpha + lambda*t^e solve for t. this involves finding an
        # e-th root of either (1/lambda) or of lambda, depending on e's sign
        ## (A sign on a ramification index ? hm)
        e = self.ramification_index
        abse = abs(e)
        lamb = S(self.xcoefficient)
        order = self.order
        if e > 0:
            phi = lamb*z**e - 1
        else:
            phi = z**abse - lamb
        mu = phi.roots(QQbar, multiplicities=False)[0]

        if all_conjugates:
            zeta_e=QQbar.zeta(abse)
            conjugates = [mu*zeta_e**k for k in range(abse)]
        else:
            conjugates = [mu]
        map(lambda x: x.exactify(), conjugates)

        # determine the resulting x-series
        xseries = []
        for c in conjugates:
            t = self.ypart.parent().gen()
            fconj = self.ypart(c*t)
            p = P(fconj(x**(QQ(1)/e)))
            p = p.add_bigoh(QQ(order+1)/abse)
            xseries.append(p)
        return xseries

    def add_term(self, order=None):
        r"""Extend the y-series terms in-place using Newton iteration.

        The modular Newtion iteration algorithm in
        :func:`newton_iteration` efficiently computes the series up to
        order :math:`t^{2^n}` where :math:`2^n` is the smallest power of
        two greater than the current order.

        """
        g,s,p = newton_iteration_step(
            self._phi, self._phiprime, self._g, self._s, self._p)

        self._g = g
        self._s = s
        self._p = p

        # operation below: yseries = ypart(y=g)(y=0)
        t = self.t
        L = self.ypart.parent()
        g = g.univariate_polynomial()(t)
        self.ypart = L(self._ypart(t,g))
        self._regular_order = self._p.degree()

    def extend(self, order=None, nterms=None):
        r"""Extends the series in place.

        Computes additional terms in the Puiseux series up to the
        specified `order` or with `nterms` number of non-zero terms. If
        neither `degree` nor `nterms` are provided then the next
        non-zero term will be added to this t-series.

        Remember that :meth:`add_term` updates `self.order` in-place.

        Parameters
        ----------
        order : int, optional
            The desired degree to extend the series to.
        nterms : int, optional
            The desired number of non-zero terms to extend the series to.

        Returns
        -------
        None

        """
        # order takes precedence
        if order:
            while self.order < order:
                self.add_term()
        elif nterms:
            while self.nterms < nterms:
                self.add_term()
        else:
            # if neither order or nterms is given, just call add_term
            self.add_term()

    def extend_to_t(self, t, curve_tol=1e-8):
        r"""Extend the series to accurately determine the y-values at `t`.

        Add terms to the t-series until the the regular place
        :math:`(x(t), y(t))` is within a particular tolerance of the
        curve that the Puiseux series is approximating.

        Parameters
        ----------
        t : complex
        eps : double
        curve_tol : double
            The tolerance for the corresponding point to lie on the curve.

        Returns
        -------
        none
            The PuiseuxTSeries is modified in-place.

        Note
        ----
        This doesn't work well in the infinite case. (Puiseux series centered
        at x=oo.)
        """
        num_iter = 0
        max_iter = 16
        while num_iter < max_iter:
            xt = self.eval_x(t)
            yt = self.eval_y(t)
            n,a = max(self.terms)
            curve_error = abs(self.f(xt,yt))
            if (curve_error < curve_tol):
                break
            else:
                self.add_term()
                num_iter += 1

    def extend_to_x(self, x, curve_tol=1e-8):
        r"""Extend the series to accurately determine the y-values at `x`.

        Add terms to the t-series until the the regular place :math:`(x,
        y)` is within a particular tolerance of the curve that the
        Puiseux series is approximating.

        Parameters
        ----------
        x : complex
        curve_tol : double
            The tolerance for the corresponding point to lie on the curve.

        Returns
        -------
        none
            The PuiseuxTSeries is modified in-place.
        """
        # simply convert to t and pass to extend. choose any conjugate since
        # the convergence rates between each conjugate is equal
        center, xcoefficient, ramification_index = self.xdata
        t = numpy.power((x-center)/xcoefficient, 1.0/ramification_index)
        self.extend_to_t(t, curve_tol=curve_tol)

    def eval_x(self, t):
        r"""Evaluate the x-part of the Puiseux series at `t`.

        Parameters
        ----------
        t : sympy.Expr or complex

        Returns
        -------
        val = complex

        """
        try:
            center, xcoefficient, ramification_index = self.xdata
            val = center + xcoefficient*t**ramification_index
        except ZeroDivisionError:
            val = infinity
        return val

    def eval_dxdt(self, t):
        r"""Evaluate the derivative of the x-part of the Puiseux series at 't'.

        Parameters
        ----------
        t : complex

        Returns
        -------
        val : complex
        """
        try:
            center, xcoefficient, ramification_index = self.xdata
            val = xcoefficient*ramification_index*t**(ramification_index-1)
        except ZeroDivisionError:
            val = infinity
        return val

    def eval_y(self, t, order=None):
        r"""Evaluate of the y-part of the Puiseux series at `t`.

        The y-part can be evaluated up to a certain order or with a
        certain number of terms.

        Parameters
        ----------
        t : complex
        nterms : int, optional
            If provided, only evaluates using `nterms` in the y-part of
            the series.  If set to zero, will evaluate the principal
            part of the series: the terms in the series which
            distinguishes places with the same x-projection.
        order : int, optional
            If provided, only evaluates up to `order`.

        Returns
        -------
        complex

        Notes
        -----

        This can be sped up using a Holder-like fast exponent evaluation
        trick.

        """
        if order:
            self.extend(order=order)

        # set which terms will be used for evaluation
        if order is not None and order >= 0:
            terms = [(n,alpha) for n,alpha in self.terms if n < order]
        else:
            terms = self.terms

        try:
            val = sum(alpha*t**n for n,alpha in terms)
        except ZeroDivisionError:
            val = infinity
        return val
