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

Priority:  minor  Milestone:  sage8.2 
Component:  commutative algebra  Keywords:  power series 
Cc:  robertwb, bruno  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 10 years ago by
comment:1 Changed 10 years ago by
 Component changed from PLEASE CHANGE to commutative algebra
 Description modified (diff)
 Keywords power series added
 Milestone set to sage4.6.2
 Owner changed from tbd to pernici
 Priority changed from major to minor
Changed 10 years ago by
comment:2 Changed 10 years ago by
 Description modified (diff)
comment:3 Changed 10 years ago by
comment:4 Changed 10 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 10 years ago by
 Status changed from new to needs_review
 Type changed from PLEASE CHANGE to defect
comment:6 Changed 10 years ago by
Apply only trac_10720_power_series_nth_root_2.patch
Changed 10 years ago by
comment:7 Changed 10 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 10 years ago by
comment:8 Changed 10 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 10 years ago by
 Description modified (diff)
comment:10 Changed 9 years ago by
 Type changed from defect to enhancement
This is an enhancement rather than a defect.
comment:11 Changed 9 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 9 years ago by
 Cc robertwb added
comment:13 Changed 8 years ago by
 Description modified (diff)
 Status changed from needs_review to 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 8 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 8 years ago by
 Milestone changed from sage5.11 to sage5.12
comment:16 Changed 7 years ago by
 Milestone changed from sage6.1 to sage6.2
comment:17 in reply to: ↑ 14 Changed 7 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 7 years ago by
 Milestone changed from sage6.2 to sage6.3
comment:19 Changed 7 years ago by
 Milestone changed from sage6.3 to sage6.4
comment:20 Changed 5 years ago by
 Cc bruno added
 Milestone changed from sage6.4 to 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 5 years ago by
Moreover, polynomials should provide a nth_root_series_trunc
(the same way there exists a inverse_series_trunc
).
comment:22 Changed 3 years ago by
 Description modified (diff)
 Milestone changed from sage7.3 to sage8.2
 Summary changed from nth_root in power series to nth_root for (Laurent) polynomials and power series
recycling...
comment:23 Changed 3 years ago by
 Description modified (diff)
comment:24 Changed 3 years ago by
 Description modified (diff)
comment:25 Changed 3 years ago by
 Description modified (diff)
comment:26 Changed 3 years ago by
 Description modified (diff)
comment:27 Changed 3 years ago by
 Description modified (diff)
comment:28 Changed 3 years ago by
 Branch set to u/vdelecroix/10720
 Commit set to 1c3d0388dd2b193a07bcf6a13d6cd7b8fd445160
 Description modified (diff)
 Status changed from needs_work to needs_review
 Summary changed from nth_root for (Laurent) polynomials and power series to nth_root for (Laurent) power series
comment:29 Changed 3 years ago by
 Description modified (diff)
comment:30 Changed 3 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 3 years ago by
 Status changed from needs_review to 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 3 years ago by
 Description modified (diff)
comment:33 Changed 3 years ago by
 Commit changed from 1c3d0388dd2b193a07bcf6a13d6cd7b8fd445160 to 2e432e1ffb6134b2262b8a85268b3656de1cbb8d
Branch pushed to git repo; I updated commit sha1. New commits:
2e432e1  10720: use Newton method for x^n instead

comment:34 Changed 3 years ago by
 Status changed from needs_work to needs_review
comment:35 Changed 3 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 in reply to: ↑ 14 Changed 3 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 3 years ago by
 Commit changed from 2e432e1ffb6134b2262b8a85268b3656de1cbb8d to bcade8aa53391ef8930e474281dfc9b51599bd1b
Branch pushed to git repo; I updated commit sha1. New commits:
bcade8a  10720: more examples

comment:38 Changed 3 years ago by
comment:39 Changed 3 years ago by
 Status changed from needs_review to 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:40 Changed 3 years ago by
Check that the result are consistent
>results
comment:41 Changed 3 years ago by
 Reviewers set to 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 3 years ago by
 Commit changed from bcade8aa53391ef8930e474281dfc9b51599bd1b to 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 3 years ago by
 Status changed from needs_work to needs_review
rebased and fixed. I am only setting to needs review to have some patchbot reviews.
comment:44 Changed 3 years ago by
 Status changed from needs_review to positive_review
Merci Sébastien!
comment:45 Changed 3 years ago by
 Branch changed from u/vdelecroix/10720 to eddd45db6e20c602db8c539fe6752ba259b2af9d
 Resolution set to fixed
 Status changed from positive_review to 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 .