Opened 6 years ago
Closed 6 years ago
#20571 closed enhancement (fixed)
Newton method for nth_root of polynomial
Reported by:  Vincent Delecroix  Owned by:  

Priority:  major  Milestone:  sage7.3 
Component:  algebra  Keywords:  
Cc:  Nils Bruin, Benjamin Hackl, Kiran Kedlaya  Merged in:  
Authors:  Vincent Delecroix  Reviewers:  Bruno Grenet 
Report Upstream:  N/A  Work issues:  
Branch:  24ffc9f (Commits, GitHub, GitLab)  Commit:  24ffc9fca48ab2ce4123ba83cefe5543a11114c8 
Dependencies:  Stopgaps: 
Description (last modified by )
In #20086 a nth_root
method for polynomial was implemented using factorization. But we can use Newton method for that!
It is faster than factorization even over ZZ where factor is highly optimized:
sage: x = polygen(ZZ) sage: p = x**14 + x**3  12 sage: q = p**13 sage: %timeit _ = q.nth_root(13) 1000 loops, best of 3: 416 µs per loop sage: %timeit _ = q.factor() 1000 loops, best of 3: 895 µs per loop sage: q = p**37 sage: %timeit _ = q.nth_root(37) 1000 loops, best of 3: 1.17 µs per loop sage: %timeit _ = q.factor() 100 loops, best of 3: 3.92 ms per loop
And the Newton method also works over polynomial when factorization is not implemented
sage: R1.<x> = QQ[] sage: R2.<y> = R1[] sage: R3.<z> = R2[] sage: q = (x+y+z)**3 sage: q.factor() Traceback (most recent call last): ... NotImplementedError: sage: q.nth_root(3) z + y + x
Attachments (1)
Change History (16)
comment:1 Changed 6 years ago by
comment:2 Changed 6 years ago by
Description:  modified (diff) 

comment:3 Changed 6 years ago by
Description:  modified (diff) 

comment:4 Changed 6 years ago by
We need to do something else when the characteristic divides n.
sage: R.<x>=GF(5)[] sage: nth_root((x+1)^5,5) ZeroDivisionError: Inverse does not exist.
Taking a pth root (with p the characteristic) is easy, though: all exponents that occur need to be a multiple of p and we simply divide them by p. Similarly, we take a pth root of all coefficients.
Changed 6 years ago by
Attachment:  poly_nth_root.py added 

comment:5 Changed 6 years ago by
Indeed! Thanks for having a look. I added some random tests in positive characteristic as well (and they pass).
comment:6 Changed 6 years ago by
Branch:  → u/vdelecroix/20571 

Commit:  → 530a5639eb44d9f6feccb1b4d382ca356faf3146 
Dependencies:  #20086 
Milestone:  sage7.2 → sage7.3 
Status:  new → needs_review 
comment:7 Changed 6 years ago by
Description:  modified (diff) 

comment:8 Changed 6 years ago by
Description:  modified (diff) 

comment:9 Changed 6 years ago by
Commit:  530a5639eb44d9f6feccb1b4d382ca356faf3146 → f1996ab8f035d548bec7d2e518da0081ab211b29 

Branch pushed to git repo; I updated commit sha1. New commits:
f1996ab  Trac 20571: use m.ordinal_str() and fix doctests

comment:10 Changed 6 years ago by
Cc:  Kiran Kedlaya added 

comment:11 Changed 6 years ago by
Status:  needs_review → needs_work 

Hi Vincent,
 For the Newton method, your loop ranges over
newton_method_sizes(p.degree()+1)
: if you are computing ann
th root, isn't it sufficient to loop up top.degree()//n + 1
? It would fasten the case whenp.nth_root(n)
raises anException
.
 Instead of the following code, you may use
i = self.valuation()
. The code would gain in legibility, and it is slightly faster if I ran my tests correctly:i = 1 while self[i].is_zero(): i += 1
 Two typos:
This is computed using Newton method in the ring of power series. This method works only when the base ring is an integral domain. Morever, for polynomial whose coefficient of lower degree is different from 1, the  elemehts of the base ring should have a method ``nth_root`` implemented. + elements of the base ring should have a method ``nth_root`` implemented.
 # begining of Newton method + # beginning of Newton method Sorig = p.parent()
comment:12 Changed 6 years ago by
Commit:  f1996ab8f035d548bec7d2e518da0081ab211b29 → 24ffc9fca48ab2ce4123ba83cefe5543a11114c8 

Branch pushed to git repo; I updated commit sha1. New commits:
24ffc9f  Trac 20571: review

comment:13 Changed 6 years ago by
Status:  needs_work → needs_review 

Hi Bruno,
Thanks for having a look. I implemented your suggestions.
Vincent
comment:14 Changed 6 years ago by
Reviewers:  → Bruno Grenet 

Status:  needs_review → positive_review 
Looks good to me!
comment:15 Changed 6 years ago by
Branch:  u/vdelecroix/20571 → 24ffc9fca48ab2ce4123ba83cefe5543a11114c8 

Resolution:  → fixed 
Status:  positive_review → closed 
Waiting for #20086 to be merged... in the mean time people can have a look at poly_nth_root.py.