Opened 12 years ago
Closed 5 years ago
#10720 closed enhancement (fixed)
nth_root for (Laurent) power series
Reported by:  mario pernici  Owned by:  mario pernici 

Priority:  minor  Milestone:  sage8.2 
Component:  commutative algebra  Keywords:  power series 
Cc:  Robert Bradshaw, Bruno Grenet  Merged in:  
Authors:  Mario Pernici, Vincent Delecroix  Reviewers:  Sébastien Labbé 
Report Upstream:  N/A  Work issues:  
Branch:  eddd45d (Commits, GitHub, GitLab)  Commit:  eddd45db6e20c602db8c539fe6752ba259b2af9d 
Dependencies:  Stopgaps: 
Description (last modified by )
There is a nth_root method defined on univariate polynomial (via Newton method)
sage: R.<x> = QQ[] sage: ((1 + x  x^2)**5).nth_root(5) x^2 + x + 1
We provide a more general implementation in a new method _nth_root_series
that compute the series expansion of the nth root for univariate polynomials. Using it we implement straightforward nth_root
for univariate (Laurent) power series.
This branch will not consider support for extend=True
(see this sagedevel thread). When extend=True
the method will simply raise a NotImplementedError
while waiting for Puiseux series in Sage (see #4618).
On multivariate polynomials there is also a nth_root
method but which is implemented via factorization (sic)! The multivariate case should just call the univariate case with appropriate variable ordering. This will be dealt with in another ticket.
Attachments (4)
Change History (49)
Changed 12 years ago by
Attachment:  trac_10720_power_series_nth_root.patch added 

comment:1 Changed 12 years ago by
Authors:  → mario pernici 

Component:  PLEASE CHANGE → commutative algebra 
Description:  modified (diff) 
Keywords:  power series added 
Milestone:  → sage4.6.2 
Owner:  changed from tbd to mario pernici 
Priority:  major → minor 
Changed 12 years ago by
Attachment:  trac_10720_power_series_nth_root_2.patch added 

comment:2 Changed 12 years ago by
Description:  modified (diff) 

comment:3 Changed 12 years ago by
comment:4 Changed 12 years ago by
in the above message there is a wrong example on where sqrt fails; here is the corrected example
sage: R.<x,y> = QQ[] sage: S.<t> = R[[]] sage: p = 4 + t*x + t^2*y + O(t^3) sage: p.nth_root(2) 2 + 1/4*x*t + (1/64*x^2 + 1/4*y)*t^2 + O(t^3) sage: p.sqrt() Traceback (most recent call last): ... AttributeError: 'sage.rings.polynomial.multi_polynomial_libsingular.MPolynomial_libsingular' object has no attribute 'sqrt' }
comment:5 Changed 12 years ago by
Status:  new → needs_review 

Type:  PLEASE CHANGE → defect 
Changed 12 years ago by
Attachment:  trac_10720_power_series_nth_root_3.patch added 

comment:7 Changed 12 years ago by
Description:  modified (diff) 

With trac_10720_power_series_nth_root_3.patch nth_root is much faster for n
large
sage: S.<t> = QQ[[]] sage: p = t.exp(30) sage: %time p1 = p.nth_root(1000) Wall time: 0.01 s
with the previous version it takes 84s
Changed 12 years ago by
Attachment:  trac_10720_power_series_nth_root_4.patch added 

comment:8 Changed 12 years ago by
With trac_10720_power_series_nth_root_4.patch nth_root and inversion become faster for series on polynomial rings, when used together with trac_10480_fast_PowerSeries_poly_multiplicationalternative.patch from ticket #10480 (otherwise these functions perform as without the patch).
In the following benchmarks apply trac_karatsuba_improvements.patch and trac_10480_fast_PowerSeries_poly_multiplicationalternative.patch, trac_10720_power_series_nth_root_2.patch and trac_10720_power_series_nth_root_3.patch; in (2) apply also the new patch
sage: N = 2 sage: R= PolynomialRing(QQ,N,'x') sage: x = R.gens() sage: S.<t> = R[[]] sage: sx = sum(x) sage: prec = 15 sage: p = (t.exp(prec)).subs({t:t*sx}) sage: %time p1 = p^1 sage: %time p1 = p.nth_root(2)
Inversion benchmark:
N prec (1) (2) 2 15 0.06 0.01 3 15 2.1 0.07 4 15 58 0.5 5 15 1440 3.4
Square root benchmark:
N prec (1) (2) 2 15 0.15 0.03 3 15 5.2 0.31 4 15 161 3.3
The speedup is large because the terms with high degree in t
are large polynomials, due to the combinatorial
possibilities with several variables; using _mul_ instead of _mul_trunc
several products involving these large polynomials are computed and
later truncated away.
comment:9 Changed 12 years ago by
Description:  modified (diff) 

comment:10 Changed 11 years ago by
Type:  defect → enhancement 

This is an enhancement rather than a defect.
comment:11 Changed 11 years ago by
Which patches need to be applied? Please put some information in the ticket description.
If trac_10720_power_series_nth_root_3.patch
is meant to be applied, please add a doctest to __pow_trunc
.
comment:12 Changed 11 years ago by
Cc:  Robert Bradshaw added 

comment:13 Changed 10 years ago by
Description:  modified (diff) 

Status:  needs_review → needs_work 
As Volker mentioned, __pow_trunc
lacks a docstring. (Just to remove this from the list of tickets needing review until that is fixed.)
comment:14 followups: 17 36 Changed 10 years ago by
One related remark:
Thanks to ticket #14115, it is now easy to implement a general rational power by
def rational_power(f,alpha): return (f.log()*alpha).exp()
I do not know how it compares in terms of speed with the method of this ticket.
comment:15 Changed 9 years ago by
Milestone:  sage5.11 → sage5.12 

comment:16 Changed 9 years ago by
Milestone:  sage6.1 → sage6.2 

comment:17 Changed 9 years ago by
Replying to chapoton:
Thanks to ticket #14115, it is now easy to implement a general rational power by
def rational_power(f,alpha): return (f.log()*alpha).exp()I do not know how it compares in terms of speed with the method of this ticket.
Additionally #15601 implements Pari series. If a speed comparison is made please include that implementation as third case to compare.
comment:18 Changed 9 years ago by
Milestone:  sage6.2 → sage6.3 

comment:19 Changed 8 years ago by
Milestone:  sage6.3 → sage6.4 

comment:20 Changed 6 years ago by
Cc:  Bruno Grenet added 

Milestone:  sage6.4 → sage7.3 
We should use (and possibly improve) #20571. Would be a good idea to use this ticket to implement a _pow_trunc_
following the __pow_trunc
from trac_10720_power_series_nth_root_3.patch.
comment:21 Changed 6 years ago by
Moreover, polynomials should provide a nth_root_series_trunc
(the same way there exists a inverse_series_trunc
).
comment:22 Changed 5 years ago by
Description:  modified (diff) 

Milestone:  sage7.3 → sage8.2 
Summary:  nth_root in power series → nth_root for (Laurent) polynomials and power series 
recycling...
comment:23 Changed 5 years ago by
Description:  modified (diff) 

comment:24 Changed 5 years ago by
Description:  modified (diff) 

comment:25 Changed 5 years ago by
Description:  modified (diff) 

comment:26 Changed 5 years ago by
Description:  modified (diff) 

comment:27 Changed 5 years ago by
Description:  modified (diff) 

comment:28 Changed 5 years ago by
Authors:  mario pernici → Vincent Delecroix 

Branch:  → u/vdelecroix/10720 
Commit:  → 1c3d0388dd2b193a07bcf6a13d6cd7b8fd445160 
Description:  modified (diff) 
Status:  needs_work → needs_review 
Summary:  nth_root for (Laurent) polynomials and power series → nth_root for (Laurent) power series 
comment:29 Changed 5 years ago by
Description:  modified (diff) 

comment:30 Changed 5 years ago by
Reproducing the example of 7 comment:7
sage: S.<t> = QQ[[]] sage: p = t.exp(30) sage: %time p1 = p.nth_root(1000) CPU times: user 11.7 ms, sys: 3.27 ms, total: 15 ms Wall time: 13.8 ms
comment:31 Changed 5 years ago by
Status:  needs_review → needs_work 

The initial proposal of using Newton method for y^n
is more clever since we don't need to use power series inversion...
comment:32 Changed 5 years ago by
Description:  modified (diff) 

comment:33 Changed 5 years ago by
Commit:  1c3d0388dd2b193a07bcf6a13d6cd7b8fd445160 → 2e432e1ffb6134b2262b8a85268b3656de1cbb8d 

Branch pushed to git repo; I updated commit sha1. New commits:
2e432e1  10720: use Newton method for x^n instead

comment:34 Changed 5 years ago by
Status:  needs_work → needs_review 

comment:35 Changed 5 years ago by
It is indeed better
sage: S.<t> = QQ[[]] sage: p = t.exp(30) sage: %time p1 = p.nth_root(1000) CPU times: user 3.8 ms, sys: 30 µs, total: 3.83 ms Wall time: 3.69 ms
comment:36 Changed 5 years ago by
Replying to chapoton:
One related remark:
Thanks to ticket #14115, it is now easy to implement a general rational power by
def rational_power(f,alpha): return (f.log()*alpha).exp()I do not know how it compares in terms of speed with the method of this ticket.
sage: R.<x> = PowerSeriesRing(QQ, default_prec=100) sage: p = (1 + 2*x  x^4)**200 sage: %timeit p1 = p.nth_root(1000) 100 loops, best of 3: 4.18 ms per loop sage: %timeit p2 = (p.log()/1000).exp() 100 loops, best of 3: 5.84 ms per loop
But the main advantage of the method proposed here is that it works in positive characteristic.
comment:37 Changed 5 years ago by
Commit:  2e432e1ffb6134b2262b8a85268b3656de1cbb8d → bcade8aa53391ef8930e474281dfc9b51599bd1b 

Branch pushed to git repo; I updated commit sha1. New commits:
bcade8a  10720: more examples

comment:38 Changed 5 years ago by
Authors:  Vincent Delecroix → Mario Pernici, Vincent Delecroix 

comment:39 Changed 5 years ago by
Status:  needs_review → needs_work 

def _nth_root_series(self, long n, long prec, start=None):
is missing a INPUT block, in which, in particular, the goal of inputstart
should be explained.
 At two places the input block says
 ``prec``  integer (optional)  precision of the result
, but the code doesprec = min(self.prec(), prec)
. I suggest to add a sentence advertising it saying something likeIf prec is larger than self.prec(), then self.prec() is used.
comment:41 Changed 5 years ago by
Reviewers:  → Sébastien Labbé 

Once, those three changes on the documentation aspect are done, feel free to change the status to positive review on my behalf.
comment:42 Changed 5 years ago by
Commit:  bcade8aa53391ef8930e474281dfc9b51599bd1b → eddd45db6e20c602db8c539fe6752ba259b2af9d 

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
0d7545d  10720: _nth_root_series for univariate polynomials

fc2bb60  10720: nth_root for (Laurent) power series

c83d4c7  10720: use Newton method for x^n instead

3a5d992  10720: more examples

eddd45d  10720: doc fixes

comment:43 Changed 5 years ago by
Status:  needs_work → needs_review 

rebased and fixed. I am only setting to needs review to have some patchbot reviews.
comment:45 Changed 5 years ago by
Branch:  u/vdelecroix/10720 → eddd45db6e20c602db8c539fe6752ba259b2af9d 

Resolution:  → fixed 
Status:  positive_review → closed 
The nthroot of power series
y = x^n
is computed using the Newton method fory = x^n
With this patch
nth_root
works also for series onZZ
there is a bug in
sqrt
in this case.nth_root
can be used to compute the square root of series which currentlysqrt
does not supportnthroot
iteration formula is divisionfree, which is convenient in this case; in this casesqrt
could callnth_root
.nth_root
can be used in the multivariate series considered in ticket #1956 .