Opened 12 years ago

Closed 5 years ago

# nth_root for (Laurent) power series

Reported by: Owned by: mario pernici mario pernici minor sage-8.2 commutative algebra power series Robert Bradshaw, Bruno Grenet Mario Pernici, Vincent Delecroix Sébastien Labbé N/A eddd45d eddd45db6e20c602db8c539fe6752ba259b2af9d

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.

### comment:1 Changed 12 years ago by mario pernici

Authors: → mario pernici PLEASE CHANGE → commutative algebra modified (diff) power series added → sage-4.6.2 changed from tbd to mario pernici major → minor

### comment:2 Changed 12 years ago by mario pernici

Description: modified (diff)

### comment:3 Changed 12 years ago by mario 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 12 years ago by mario 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 12 years ago by Robert Bradshaw

Status: new → needs_review PLEASE CHANGE → defect

### comment:6 Changed 12 years ago by Robert Bradshaw

Apply only trac_10720_power_series_nth_root_2.patch

### comment:7 Changed 12 years ago by mario 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

### comment:8 Changed 12 years ago by mario 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 12 years ago by mario pernici

Description: modified (diff)

### comment:10 Changed 11 years ago by Johan Bosman

Type: defect → enhancement

This is an enhancement rather than a defect.

### comment:11 Changed 11 years ago by Volker Braun

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:13 Changed 10 years ago by Julian Rüth

Description: modified (diff) 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 follow-ups:  17  36 Changed 10 years ago by Frédéric 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 9 years ago by Jeroen Demeyer

Milestone: sage-5.11 → sage-5.12

### comment:16 Changed 9 years ago by For batch modifications

Milestone: sage-6.1 → sage-6.2

### comment:17 in reply to:  14 Changed 9 years ago by Ralf Stephan

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 For batch modifications

Milestone: sage-6.2 → sage-6.3

### comment:19 Changed 8 years ago by For batch modifications

Milestone: sage-6.3 → sage-6.4

### comment:20 Changed 6 years ago by Vincent Delecroix

Cc: Bruno Grenet added sage-6.4 → 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 6 years ago by Vincent Delecroix

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 Vincent Delecroix

Description: modified (diff) sage-7.3 → sage-8.2 nth_root in power series → nth_root for (Laurent) polynomials and power series

recycling...

### comment:23 Changed 5 years ago by Vincent Delecroix

Description: modified (diff)

### comment:24 Changed 5 years ago by Vincent Delecroix

Description: modified (diff)

### comment:25 Changed 5 years ago by Vincent Delecroix

Description: modified (diff)

### comment:26 Changed 5 years ago by Vincent Delecroix

Description: modified (diff)

### comment:27 Changed 5 years ago by Vincent Delecroix

Description: modified (diff)

### comment:28 Changed 5 years ago by Vincent Delecroix

Authors: mario pernici → Vincent Delecroix → u/vdelecroix/10720 → 1c3d0388dd2b193a07bcf6a13d6cd7b8fd445160 modified (diff) needs_work → needs_review nth_root for (Laurent) polynomials and power series → nth_root for (Laurent) power series

New commits:

 ​5bb0074 `10720: _nth_root_series for univariate polynomials` ​1c3d038 `10720: nth_root for (Laurent) power series`

### comment:29 Changed 5 years ago by Vincent Delecroix

Description: modified (diff)

### comment:30 Changed 5 years ago by Vincent Delecroix

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 Vincent Delecroix

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 Vincent Delecroix

Description: modified (diff)

### comment:33 Changed 5 years ago by git

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 Vincent Delecroix

Status: needs_work → needs_review

### comment:35 Changed 5 years ago by Vincent Delecroix

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

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 git

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

 ​bcade8a `10720: more examples`

### comment:38 Changed 5 years ago by Vincent Delecroix

Authors: Vincent Delecroix → Mario Pernici, Vincent Delecroix

### comment:39 Changed 5 years ago by Sébastien Labbé

Status: needs_review → 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 5 years ago by Sébastien Labbé

1. `Check that the result are consistent` -> `results`

### comment:41 Changed 5 years ago by Sébastien Labbé

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 git

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 Vincent Delecroix

Status: needs_work → needs_review

rebased and fixed. I am only setting to needs review to have some patchbot reviews.

### comment:44 Changed 5 years ago by Vincent Delecroix

Status: needs_review → positive_review

Merci Sébastien!

### comment:45 Changed 5 years ago by Volker Braun

Branch: u/vdelecroix/10720 → eddd45db6e20c602db8c539fe6752ba259b2af9d → fixed positive_review → closed
Note: See TracTickets for help on using tickets.