# -*- pyrex -*-
r"""
Newton's method for polynomials over power series.

An initial guess `x_0` is refined by the rule `x_{k+1} = x_{k} -
f(x_k)/f'(x_k)`, the prime denoting the derivative.

AUTHORS:

- Nick Alexander (original implementation)
"""

# from sage.rings.all import QQ, CDF

from sage.rings.power_series_ring_element cimport PowerSeries
# from sage.rings.all import infinity
import sage.misc.misc

def ps_roots(f, prec=5, all=False):
    r"""
    Find roots of polynomials over power series rings like `R[[t]][y]`.

    ALGORITHM:

    Newton's method.  An initial guess `x_0` is refined by the rule `x_{k+1} =
    x_{k} - f(x_k)/f'(x_k)`, the prime denoting the derivative.

    WARNING:

    .. the algorithm is not known to be numerically stable.  Inexact rings
    will suffer from the usual maladies of root finding: see :meth:`roots` for
    more information.

    AUTHORS:

    - Nick Alexander (original implementation)

    EXAMPLES::

        sage: R.<t> = QQ[[]]
        sage: S.<y> = R[]

        sage: f = y^3 - (t^2 + t)*y + t^5 + t + 1; g = ps_roots(f, prec=5)
        sage: g
        -1 - 2/3*t - 1/9*t^2 - 1/81*t^3 + 1/27*t^4 + O(t^5)
        sage: f(g)
        O(t^5)
        sage: g = ps_roots(f, prec=12)
        sage: g
        -1 - 2/3*t - 1/9*t^2 - 1/81*t^3 + 1/27*t^4 - 29/81*t^5 + 2284/6561*t^6 - 4624/19683*t^7 + 10693/59049*t^8 - 284333/1594323*t^9 + 50659/177147*t^10 - 2003644/4782969*t^11 + O(t^12)
        sage: f(g)
        O(t^12)

    This one should fail.  There's an approximation to order 0, namely `y=0 +
    O(t)`, but there is no approximation to order 1.  Correspondingly, the
    ``ValuationError`` witnesses that the derivative of `f` at 0 is 0 (to order
    1) ::

        sage: f = y^3 - (t^2 + t)*y + t^5 + t + 11; g = ps_roots(f, prec=5)
        Traceback (most recent call last):
        ...
        ValueError: There exists no approximation to order 0 -- that is, no root of y^3 + 11 defined over Rational Field

    With this one, the first order approximation `y=0 + O(t)` cannot be
    improved, since `0 + O(t)` is a root of the derivative of `y^3`, namely `3
    y^2`::

        sage: f = y^3 - (t^2 + t)*y + t^5 + t; g = ps_roots(f, prec=5)
        Traceback (most recent call last):
        ...
        ValueError: There exists no approximation to order 2

    Now let's test the ``all`` parameter.  The first difference is that the
    following doesn't raise an error, it returns an empty list::

        sage: f = y^3 - (2*t^2 + t + 2)^2; ps_roots(f, prec=10, all=True)
        []

        sage: f = y^2 - (2*t^2 + t + 2)^2; gs = ps_roots(f, prec=10, all=True)
        sage: gs
        [2 + t + 2*t^2 + O(t^10), -2 - t - 2*t^2 + O(t^10)]
        sage: [ f(g) for g in gs ]
        [O(t^10), O(t^10)]

        sage: f = (y - (2 + t + t^5))*(y - (1 + 15*t + t^2))*(y - (3 + 10*t + t^5))
        sage: gs = ps_roots(f, prec=10, all=True); gs
        [3 + 10*t + t^5 + O(t^10), 2 + t + t^5 + O(t^10), 1 + 15*t + t^2 + O(t^10)]
        sage: [ f(g) for g in gs ]
        [O(t^10), O(t^10), O(t^10)]

    In this example, two roots agree to order 1.  That means Newton's method
    doesn't work::

        sage: f = (y - (2 + t + t^5))*(y - (2 + 15*t + t^2))*(y - (3 + 10*t + t^5))
        sage: ps_roots(f, prec=10, all=True)
        Traceback (most recent call last):
        ...
        ValueError: There exists no approximation to order 2
    """
    fp = f.derivative()
    poly_ring, y = f.parent().objgen()
    series_ring, t = poly_ring.base_ring().objgen()
    base_ring = series_ring.base_ring()

    # we start with a root of f mod t defined over R
    g = base_ring['y']([ c[0] for c in f ]) # substitute t = 0
    rts = g.roots(multiplicities=False)

    cdef PowerSeries s
    if all is False:
        if not rts:
            raise ValueError, "There exists no approximation to order 0 -- that is, no root of %s defined over %s" % (g, base_ring)
        rts = [ rts[0] ] # just choose a random root

    ss = []
    for rt in rts:
        s = series_ring(rt)
        s._prec = 1

        for cur_prec in sage.misc.misc.newton_method_sizes(prec)[1:]:
            old_prec = s._prec
            s._prec = cur_prec
            D = fp(s)
            if D.valuation() >= old_prec:
                raise ValueError, "There exists no approximation to order %s" % cur_prec

            s = s - f(s)/fp(s)
        ss.append(s)
    if all is False:
        return ss[0]
    return ss