Opened 8 years ago

Closed 19 months ago

#10720 closed enhancement (fixed)

nth_root for (Laurent) power series

Reported by: pernici Owned by: pernici
Priority: minor Milestone: sage-8.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) Commit: eddd45db6e20c602db8c539fe6752ba259b2af9d
Dependencies: Stopgaps:

Description (last modified by vdelecroix)

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 n-th 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 sage-devel thread). When extend=True the method will simply raise a NotImplementedError while waiting for Puiseux series in Sage (see #4618).

On multi-variate 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)

trac_10720_power_series_nth_root.patch (3.8 KB) - added by pernici 8 years ago.
trac_10720_power_series_nth_root_2.patch (4.9 KB) - added by pernici 8 years ago.
trac_10720_power_series_nth_root_3.patch (2.3 KB) - added by pernici 8 years ago.
trac_10720_power_series_nth_root_4.patch (2.0 KB) - added by pernici 8 years ago.

Download all attachments as: .zip

Change History (49)

Changed 8 years ago by pernici

comment:1 Changed 8 years ago by pernici

  • Authors set to mario pernici
  • Component changed from PLEASE CHANGE to commutative algebra
  • Description modified (diff)
  • Keywords power series added
  • Milestone set to sage-4.6.2
  • Owner changed from tbd to pernici
  • Priority changed from major to minor

Changed 8 years ago by pernici

comment:2 Changed 8 years ago by pernici

  • Description modified (diff)

comment:3 Changed 8 years ago by pernici

The nth-root of power series y = x^n is computed using the Newton method for y = x^-n

x' = (1+1/n)*x - y*x^(n+1)/n
sage: R.<t> = QQ[]
sage: p = (1 + 2*t + 5*t^2 + 7*t^3 + O(t^4))^3
sage: p.nth_root(3)
1 + 2*t + 5*t^2 + 7*t^3 + O(t^4)
sage: p = (1 + 2*t + 5*t^2 + 7*t^3 + O(t^4))^-3
sage: p.nth_root(-3)
1 + 2*t + 5*t^2 + 7*t^3 + O(t^4)

With this patch nth_root works also for series on ZZ

sage: R.<t> = ZZ[[]]
sage: p = 1 + 4*t^2 + 3*t^3 + O(t^4)
sage: p.nth_root(2)
1 + 2*t^2 + 3/2*t^3 + O(t^4)
sage: p.sqrt(2)
1 + 2*t^2 + 3/2*t^3 + O(t^4)
sage: p.nth_root(3)
1 + 4/3*t^2 + t^3 + O(t^4)
sage: p = 4*t^2 + 3*t^3 + O(t^4)
sage: p.nth_root(2)
2*t + 3/4*t^2 + O(t^3)
sage: p.sqrt(2)
Traceback (most recent call last):
...
TypeError: no conversion of this rational to integer

there is a bug in sqrt in this case.

nth_root can be used to compute the square root of series which currently sqrt does not support

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 = 32*t^5 + x*t^6 + y*t^7 + O(t^8)
sage: p.nth_root(5)
2*t + 1/80*x*t^2 + (-1/6400*x^2 + 1/80*y)*t^3 + O(t^4)
sage: p.sqrt()
Traceback (most recent call last):
...
AttributeError: 'sage.rings.polynomial.multi_polynomial_libsingular.MPolynomial_libsingular' object has no attribute 'sqrt'

nth-root iteration formula is division-free, which is convenient in this case; in this case sqrt could call nth_root.

nth_root can be used in the multivariate series considered in ticket #1956 .

comment:4 Changed 8 years ago by pernici

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 8 years ago by robertwb

  • Status changed from new to needs_review
  • Type changed from PLEASE CHANGE to defect

comment:6 Changed 8 years ago by robertwb

Apply only trac_10720_power_series_nth_root_2.patch

Changed 8 years ago by pernici

comment:7 Changed 8 years ago by pernici

  • 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 8 years ago by pernici

comment:8 Changed 8 years ago by pernici

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_multiplication-alternative.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_multiplication-alternative.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 8 years ago by pernici

  • Description modified (diff)

comment:10 Changed 8 years ago by johanbosman

  • Type changed from defect to enhancement

This is an enhancement rather than a defect.

comment:11 Changed 8 years ago by vbraun

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 7 years ago by robertwb

  • Cc robertwb added

comment:13 Changed 6 years ago by saraedum

  • 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 follow-ups: Changed 6 years ago by 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.

comment:15 Changed 6 years ago by jdemeyer

  • Milestone changed from sage-5.11 to sage-5.12

comment:16 Changed 5 years ago by vbraun_spam

  • Milestone changed from sage-6.1 to sage-6.2

comment:17 in reply to: ↑ 14 Changed 5 years ago by rws

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 5 years ago by vbraun_spam

  • Milestone changed from sage-6.2 to sage-6.3

comment:19 Changed 5 years ago by vbraun_spam

  • Milestone changed from sage-6.3 to sage-6.4

comment:20 Changed 3 years ago by vdelecroix

  • Cc bruno added
  • Milestone changed from sage-6.4 to sage-7.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 3 years ago by vdelecroix

Moreover, polynomials should provide a nth_root_series_trunc (the same way there exists a inverse_series_trunc).

comment:22 Changed 19 months ago by vdelecroix

  • Description modified (diff)
  • Milestone changed from sage-7.3 to sage-8.2
  • Summary changed from nth_root in power series to nth_root for (Laurent) polynomials and power series

recycling...

comment:23 Changed 19 months ago by vdelecroix

  • Description modified (diff)

comment:24 Changed 19 months ago by vdelecroix

  • Description modified (diff)

comment:25 Changed 19 months ago by vdelecroix

  • Description modified (diff)

comment:26 Changed 19 months ago by vdelecroix

  • Description modified (diff)

comment:27 Changed 19 months ago by vdelecroix

  • Description modified (diff)

comment:28 Changed 19 months ago by vdelecroix

  • Authors changed from mario pernici to Vincent Delecroix
  • 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

New commits:

5bb007410720: _nth_root_series for univariate polynomials
1c3d03810720: nth_root for (Laurent) power series

comment:29 Changed 19 months ago by vdelecroix

  • Description modified (diff)

comment:30 Changed 19 months ago by vdelecroix

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 19 months ago by vdelecroix

  • 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 19 months ago by vdelecroix

  • Description modified (diff)

comment:33 Changed 19 months ago by git

  • Commit changed from 1c3d0388dd2b193a07bcf6a13d6cd7b8fd445160 to 2e432e1ffb6134b2262b8a85268b3656de1cbb8d

Branch pushed to git repo; I updated commit sha1. New commits:

2e432e110720: use Newton method for x^-n instead

comment:34 Changed 19 months ago by vdelecroix

  • Status changed from needs_work to needs_review

comment:35 Changed 19 months ago by vdelecroix

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 19 months ago by vdelecroix

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 19 months ago by git

  • Commit changed from 2e432e1ffb6134b2262b8a85268b3656de1cbb8d to bcade8aa53391ef8930e474281dfc9b51599bd1b

Branch pushed to git repo; I updated commit sha1. New commits:

bcade8a10720: more examples

comment:38 Changed 19 months ago by vdelecroix

  • Authors changed from Vincent Delecroix to Mario Pernici, Vincent Delecroix

comment:39 Changed 19 months ago by slabbe

  • Status changed from needs_review to needs_work
  1. def _nth_root_series(self, long n, long prec, start=None): is missing a INPUT block, in which, in particular, the goal of input start should be explained.
  1. At two places the input block says - ``prec`` -- integer (optional) - precision of the result, but the code does prec = min(self.prec(), prec). I suggest to add a sentence advertising it saying something like If prec is larger than self.prec(), then self.prec() is used.

comment:40 Changed 19 months ago by slabbe

  1. Check that the result are consistent -> results

comment:41 Changed 19 months ago by slabbe

  • 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 19 months ago by git

  • Commit changed from bcade8aa53391ef8930e474281dfc9b51599bd1b to eddd45db6e20c602db8c539fe6752ba259b2af9d

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

0d7545d10720: _nth_root_series for univariate polynomials
fc2bb6010720: nth_root for (Laurent) power series
c83d4c710720: use Newton method for x^-n instead
3a5d99210720: more examples
eddd45d10720: doc fixes

comment:43 Changed 19 months ago by vdelecroix

  • 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 19 months ago by vdelecroix

  • Status changed from needs_review to positive_review

Merci Sébastien!

comment:45 Changed 19 months ago by vbraun

  • Branch changed from u/vdelecroix/10720 to eddd45db6e20c602db8c539fe6752ba259b2af9d
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.