# -*- 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. = QQ[[]]
sage: S. = 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