Ticket #20571: poly_nth_root.py

File poly_nth_root.py, 4.5 KB (added by Vincent Delecroix, 7 years ago)
Line 
1def nth_root(p, n):
2    r"""
3    Return a ``n``-th root of the polynomial ``p`` or raise a ``ValueError`` if
4    there is no such root.
5
6    EXAMPLES::
7
8        sage: R.<x> = QQ[]
9
10        sage: nth_root((x-2)**2, 2)
11        -x + 2
12        sage: nth_root((3*x+4)**2, 2)
13        3*x + 4
14        sage: nth_root(x^3 - 1, 2)
15        Traceback (most recent call last):
16        ...
17        ValueError: not a 2-th power
18        sage: nth_root(x^2 + 4, 2)
19        Traceback (most recent call last):
20        ...
21        ValueError: not a 2-th power
22
23        sage: p = (x+1)**20 + x^20
24        sage: nth_root(p, 20)
25        Traceback (most recent call last):
26        ...
27        ValueError: not a 20-th power
28
29        sage: nth_root(x^12, 6)
30        x^2
31        sage: nth_root(x^13, 6)
32        Traceback (most recent call last):
33        ...
34        ValueError: not a 6-th power
35        sage: nth_root(x^12 * (x+1)**6, 6)
36        x^3 + x^2
37
38        sage: for _ in range(100):
39        ....:     p = R.random_element(degree=randint(10,20))
40        ....:     n = ZZ.random_element(2,20)
41        ....:     r = nth_root(p**n, n)
42        ....:     assert r == p or r == -p, "r = {}\np = {}".format(r,p)
43
44    The algorithm works for any domain that implements a n-th root method. In
45    particular for multivariate polynomials::
46
47        sage: R1.<x> = ZZ[]
48        sage: R2.<y> = R1[]
49        sage: p = x*y + 1
50        sage: nth_root(p**5, 5) == p
51        True
52        sage: p = y**2 * (x**3 + x - 3) + 13*y + 1
53        sage: nth_root(p**13, 13) == p
54        True
55        sage: p = (x*y+1)^20 + y^20
56        sage: nth_root(p, 20)
57        Traceback (most recent call last):
58        ...
59        ValueError: not a 20-th power
60
61        sage: R.<x> = QQbar[]
62        sage: p = x**3 + QQbar(2).sqrt() * x - QQbar(3).sqrt()
63        sage: r = nth_root(p**5, 5)
64        sage: r * p[0] == p * r[0]
65        True
66        sage: p = (x+1)^20 + x^20
67        sage: nth_root(p, 20)
68        Traceback (most recent call last):
69        ...
70        ValueError: not a 20-th power
71
72        sage: z = GF(4).gen()
73        sage: R.<x> = GF(4)[]
74        sage: p = z*x**4 + 2*x - 1
75        sage: r = nth_root(p**15, 15)
76        sage: r * p[0] == p * r[0]
77        True
78        sage: nth_root((x+1)**2, 2)
79        x + 1
80        sage: nth_root((x+1)**4, 4)
81        x + 1
82        sage: nth_root((x+1)**12, 12)
83        x + 1
84        sage: nth_root(x^3 + 1, 2)
85        Traceback (most recent call last):
86        ...
87        ValueError: not a 2-th power
88        sage: p = (x+1)^17 + x^17
89        sage: r = nth_root(p, 17)
90        Traceback (most recent call last):
91        ...
92        ValueError: not a 17-th power
93
94        sage: for _ in range(100):
95        ....:     p = R.random_element(degree=randint(10,20))
96        ....:     n = ZZ.random_element(2,20)
97        ....:     r = nth_root(p**n, n)
98        ....:     assert r[0] * p == r * p[0], "r = {}\np = {}".format(r,p)
99    """
100    n = ZZ(n)
101    if p.is_zero() or p.is_one() or n.is_one():
102        return p
103    if n <= 0:
104        raise ValueError("n = {} must be a positive integer".format(n))
105    d = p.degree()
106    if d%n:
107        raise ValueError("not a {}-th power".format(n))
108
109    if p[0].is_zero():
110        # p = x^k q
111        # p^(1/n) = x^(k/n) q^(1/n)
112        i = 1
113        while p[i].is_zero():
114            i += 1
115        if i%n:
116            raise ValueError("not a {}-th power".format(n))
117        S = p.parent()
118        return S.gen()**(i//n) * nth_root(p>>i, n)
119
120    R = p.base_ring()
121
122    c = R.characteristic()
123    if c and not n%c:
124        # characteristic divides n
125        e = n.valuation(c)
126        q = c**e
127        ans = {}
128        for i in range(p.degree()+1):
129            if p[i]:
130                if i%q:
131                    raise ValueError("not a {}-th power")
132                ans[i//q] = p[i].nth_root(q)
133        p = p.parent()(ans)
134        n = n // q
135        if n.is_one():
136            return p
137
138    R = R.fraction_field()
139    Sorig = p.parent()
140    p = p.change_ring(R)
141    S = p.parent()
142    if p[0].is_one():
143        q = S.one()
144    else:
145        q = S(p[0].nth_root(n))
146
147    from sage.misc.misc import newton_method_sizes
148    for i in newton_method_sizes(p.degree()+1):
149        qi = (q**(n-1)).inverse_series_trunc(i)
150        q = ((n-1) * q + qi._mul_trunc_(p,i)) / n
151
152        # NOTE: if we knew that p was a n-th power we could remove the check
153        # below and just return q after the whole loop
154        r = q**n - p
155        if not r:
156            return Sorig(q)
157        elif not r.truncate(i).is_zero():
158            raise ValueError("not a {}-th power".format(n))
159
160    raise ValueError("not a {}-th power".format(n))