Ticket #6397: trac_6397-power-series-newton-ncalexan-1.pyx

File trac_6397-power-series-newton-ncalexan-1.pyx, 4.3 KB (added by ncalexan, 10 years ago)
Line 
1# -*- pyrex -*-
2r"""
3Newton's method for polynomials over power series.
4
5An initial guess `x_0` is refined by the rule `x_{k+1} = x_{k} -
6f(x_k)/f'(x_k)`, the prime denoting the derivative.
7
8AUTHORS:
9
10- Nick Alexander (original implementation)
11"""
12
13# from sage.rings.all import QQ, CDF
14
15from sage.rings.power_series_ring_element cimport PowerSeries
16# from sage.rings.all import infinity
17import sage.misc.misc
18
19def ps_roots(f, prec=5, all=False):
20    r"""
21    Find roots of polynomials over power series rings like `R[[t]][y]`.
22
23    ALGORITHM:
24
25    Newton's method.  An initial guess `x_0` is refined by the rule `x_{k+1} =
26    x_{k} - f(x_k)/f'(x_k)`, the prime denoting the derivative.
27
28    WARNING:
29
30    .. the algorithm is not known to be numerically stable.  Inexact rings
31    will suffer from the usual maladies of root finding: see :meth:`roots` for
32    more information.
33
34    AUTHORS:
35
36    - Nick Alexander (original implementation)
37
38    EXAMPLES::
39
40        sage: R.<t> = QQ[[]]
41        sage: S.<y> = R[]
42
43        sage: f = y^3 - (t^2 + t)*y + t^5 + t + 1; g = ps_roots(f, prec=5)
44        sage: g
45        -1 - 2/3*t - 1/9*t^2 - 1/81*t^3 + 1/27*t^4 + O(t^5)
46        sage: f(g)
47        O(t^5)
48        sage: g = ps_roots(f, prec=12)
49        sage: g
50        -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)
51        sage: f(g)
52        O(t^12)
53
54    This one should fail.  There's an approximation to order 0, namely `y=0 +
55    O(t)`, but there is no approximation to order 1.  Correspondingly, the
56    ``ValuationError`` witnesses that the derivative of `f` at 0 is 0 (to order
57    1) ::
58
59        sage: f = y^3 - (t^2 + t)*y + t^5 + t + 11; g = ps_roots(f, prec=5)
60        Traceback (most recent call last):
61        ...
62        ValueError: There exists no approximation to order 0 -- that is, no root of y^3 + 11 defined over Rational Field
63
64    With this one, the first order approximation `y=0 + O(t)` cannot be
65    improved, since `0 + O(t)` is a root of the derivative of `y^3`, namely `3
66    y^2`::
67
68        sage: f = y^3 - (t^2 + t)*y + t^5 + t; g = ps_roots(f, prec=5)
69        Traceback (most recent call last):
70        ...
71        ValueError: There exists no approximation to order 2
72
73    Now let's test the ``all`` parameter.  The first difference is that the
74    following doesn't raise an error, it returns an empty list::
75
76        sage: f = y^3 - (2*t^2 + t + 2)^2; ps_roots(f, prec=10, all=True)
77        []
78
79        sage: f = y^2 - (2*t^2 + t + 2)^2; gs = ps_roots(f, prec=10, all=True)
80        sage: gs
81        [2 + t + 2*t^2 + O(t^10), -2 - t - 2*t^2 + O(t^10)]
82        sage: [ f(g) for g in gs ]
83        [O(t^10), O(t^10)]
84
85        sage: f = (y - (2 + t + t^5))*(y - (1 + 15*t + t^2))*(y - (3 + 10*t + t^5))
86        sage: gs = ps_roots(f, prec=10, all=True); gs
87        [3 + 10*t + t^5 + O(t^10), 2 + t + t^5 + O(t^10), 1 + 15*t + t^2 + O(t^10)]
88        sage: [ f(g) for g in gs ]
89        [O(t^10), O(t^10), O(t^10)]
90
91    In this example, two roots agree to order 1.  That means Newton's method
92    doesn't work::
93
94        sage: f = (y - (2 + t + t^5))*(y - (2 + 15*t + t^2))*(y - (3 + 10*t + t^5))
95        sage: ps_roots(f, prec=10, all=True)
96        Traceback (most recent call last):
97        ...
98        ValueError: There exists no approximation to order 2
99    """
100    fp = f.derivative()
101    poly_ring, y = f.parent().objgen()
102    series_ring, t = poly_ring.base_ring().objgen()
103    base_ring = series_ring.base_ring()
104
105    # we start with a root of f mod t defined over R
106    g = base_ring['y']([ c[0] for c in f ]) # substitute t = 0
107    rts = g.roots(multiplicities=False)
108
109    cdef PowerSeries s
110    if all is False:
111        if not rts:
112            raise ValueError, "There exists no approximation to order 0 -- that is, no root of %s defined over %s" % (g, base_ring)
113        rts = [ rts[0] ] # just choose a random root
114
115    ss = []
116    for rt in rts:
117        s = series_ring(rt)
118        s._prec = 1
119
120        for cur_prec in sage.misc.misc.newton_method_sizes(prec)[1:]:
121            old_prec = s._prec
122            s._prec = cur_prec
123            D = fp(s)
124            if D.valuation() >= old_prec:
125                raise ValueError, "There exists no approximation to order %s" % cur_prec
126
127            s = s - f(s)/fp(s)
128        ss.append(s)
129    if all is False:
130        return ss[0]
131    return ss