nth_root in power series

[pernici]

computation of an nth root of a power series using the Newton method

---

Apply

1. trac_10720_power_series_nth_root_2.patch
2. trac_10720_power_series_nth_root_3.patch
3. trac_10720_power_series_nth_root_4.patch

[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 .

[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'
}

[robertwb]

Apply only trac_10720_power_series_nth_root_2.patch

[pernici]

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

[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.

[johanbosman]

This is an enhancement rather than a defect.

[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.

[saraedum]

As Volker mentioned, __pow_trunc lacks a docstring. (Just to remove this from the list of tickets needing review until that is fixed.)

[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.