Opened 13 years ago

Closed 12 years ago

# generic power series reversion

Reported by: Owned by: William Stein Alex Ghitza major sage-4.6.2 algebra lagrange, reversion Francis Clarke sage-4.6.2.alpha1 Niles Johnson Francis Clarke N/A

### Description (last modified by Niles Johnson)

From this sage-support thread: Make the following work over any base ring:

```sage: R.<x> = QQ[[]]
sage: f = 1/(1-x) - 1; f
x + x^2 + x^3 + x^4 + x^5 + x^6 + x^7 + x^8 + x^9 + x^10 + x^11 + x^12
+ x^13 + x^14 + x^15 + x^16 + x^17 + x^18 + x^19 + O(x^20)
sage: g = f.reversion(); g
x - x^2 + x^3 - x^4 + x^5 - x^6 + x^7 - x^8 + x^9 - x^10 + x^11 - x^12
+ x^13 - x^14 + x^15 - x^16 + x^17 - x^18 + x^19 + O(x^20)
sage: f(g)
x + O(x^20)
```

Matt Bainbridge says about power series reversion, which uses pari in some cases, and maybe isn't there in others:

```Its easy enough to code this in sage.  This seems to work over any
field:

def ps_inverse(f):
if f.prec() is infinity:
raise ValueError, "series must have finite precision for
reversion"
if f.valuation() != 1:
raise ValueError, "series must have valuation one for
reversion"
t = parent(f).gen()
a = 1/f.coefficients()[0]
g = a*t
for i in range(2, f.prec()):
g -=  ps_coefficient((f + O(t^(i+1)))(g),i)*a*t^i
g += O(t^f.prec())
return g

def ps_coefficient(f,i):
if i >= f.prec():
raise ValueError, "that coefficient is undefined"
else:
return f.padded_list(f.prec())[i]
```

## Apply

1. trac_7644_reversion_lagrange_6.patch

### comment:1 Changed 13 years ago by Minh Van Nguyen

Description: modified (diff)

### comment:2 follow-up:  7 Changed 13 years ago by Francis Clarke

Lagrange inversion is significantly faster. With

```def ps_inverse_Lagrange(f):
if f.valuation() != 1:
raise ValueError, "series must have valuation one for reversion"
if f.prec() is infinity:
raise ValueError, "series must have finite precision for reversion"
t = parent(f).gen()
h = t/f
k = 1
g = 0
for i in range(1, f.prec()):
k *= h
g += (1/i)*ps_coefficient(k, i - 1)*t^i
g += O(t^f.prec())
return g
```

I found

```sage: R.<x> = QQ[[]]
sage: f = exp(x) - 1
sage: timeit('ps_inverse(f)')
5 loops, best of 3: 552 ms per loop
sage: timeit('ps_inverse_Lagrange(f)')
5 loops, best of 3: 74 ms per loop
```

### comment:3 Changed 13 years ago by William Stein

```The code doesn't quite run, since it references some other function
(ps_coefficient); here's an updated version which uses only built-in
functions:

def ps_inverse_Lagrange(f):
if f.valuation() != 1:
raise ValueError, "series must have valuation one for
reversion"
if f.prec() is infinity:
raise ValueError, "series must have finite precision for
reversion"
t = parent(f).gen()
h = t/f
k = 1
g = 0
for i in range(1, f.prec()):
k *= h
g += (1/i)*k.padded_list(i)[i - 1]*t^i
g += O(t^f.prec())
return g
```

### comment:4 follow-up:  5 Changed 13 years ago by Alex Ghitza

What should we do with power series with coefficients in, say, ZZ? Raise an error, or return a power series over the fraction field?

### comment:5 in reply to:  4 Changed 13 years ago by Francis Clarke

Replying to AlexGhitza:

What should we do with power series with coefficients in, say, ZZ? Raise an error, or return a power series over the fraction field?

In general, return a power series over the fraction field. But if the leading coefficient is a unit, then despite the fact that Lagrange inversion involves division, the inverse series has coefficients in the same ring as the original series. E.g., with the function defined in was,

```sage: PS.<t> = ZZ[[]]
sage: f = t + t^2 + O(t^10)
sage: ps_inverse_Lagrange(f)
t - t^2 + 2*t^3 - 5*t^4 + 14*t^5 - 42*t^6 + 132*t^7 - 429*t^8 + 1430*t^9 + O(t^10)
```

though

```sage: ps_inverse_Lagrange(f).parent()
Power Series Ring in t over Rational Field
```

Over a ring of finite characteristic, to use Lagrange inversion, you have to lift to a ring of characteristic zero, invert, and then project down to the original ring.

### comment:6 Changed 12 years ago by Niles Johnson

Status: new → needs_review

I've uploaded a patch which implements ps_inverse_Lagrange from above. For simplicity, I didn't implement over rings of positive characteristic, or in the case that the leading coefficient is not a unit. Instead, I included examples of changing the base ring and carrying out the reversion there.

Note: the line

```g += (1/i)*k....
```

gave me some trouble, since I didn't realize .pyx files use python's integer division; someone who understands the implicit conversions between sage and python might want to check I dealt with that correctly. It now reads:

```g += k.padded_list(i)[i - 1]/i*t**i
```

and works fine in all my tests.

### comment:7 in reply to:  2 Changed 12 years ago by Niles Johnson

Authors: → niles Francis Clarke added lagrange reversion added

Replying to fwclarke:

Hi Francis,

You may be interested in reviewing this patch since it's based on your code. I believe it will be an easy review.

best, Niles

### comment:8 Changed 12 years ago by Niles Johnson

Authors: niles → Niles Johnson modified (diff)

### comment:9 Changed 12 years ago by Niles Johnson

The latest version of this patch applies cleanly to sage 4.6.1.alpha2 and passes all (-long) doctests. It also adds `power_series_poly` to the sage documentation, and makes some minor changes so that all documentation builds cleanly.

### comment:10 follow-up:  11 Changed 12 years ago by Francis Clarke

I've been looking at this. One concern is that for power series with rational coefficients the existing method using pari is a great deal faster, so that at least in this case the pari method should be retained. Strangely the existing method fails for integer power series, though pari handles them perfectly ok.

On investigation it turns out that it is the conversion to pari which is failing and that this problem is pointed out in modular/overconvergent/genus0.py, and raised as #4376.  I have produced a short patch (awaiting review) which extends the range of base rings over which pari conversion of a power series is possible.

What I propose then is that a revised patch (depending on the #4376 patch) should try to use pari and that only if this fails should the Lagrange inversion code be used.

I also think it is sensible to be able to revert infinite precision series, either by specifying the desired precision or by using the parent's default precision.

### Changed 12 years ago by Niles Johnson

Attachment: trac_7644_reversion_lagrange_3.patch​ added

apply only this patch; tries to use pari first, then falls back to Lagrange inversion

### comment:11 in reply to:  10 ; follow-up:  14 Changed 12 years ago by Niles Johnson

Description: modified (diff)

Replying to fwclarke:

I've been looking at this.

Thanks!

One concern is that for power series with rational coefficients the existing method using pari is a great deal faster, so that at least in this case the pari method should be retained.

Yes, especially if there is work in progress to support converting more rings to pari. I wrote a revised patch which first attempts to convert to pari and do reversion there, and then tries the Lagrange inversion if conversion to pari fails. I think that implementation means that this patch can be independent of #4376 -- it will simply perform better when that patch is included.

I also think it is sensible to be able to revert infinite precision series, either by specifying the desired precision or by using the parent's default precision.

Yes, upon further consideration I agree. I've made this and two other improvements:

1. Given a power series with infinite precision and degree `deg`, it's reversion is computed with precision `deg + 1`.
2. Given a power series whose leading coefficient is not a unit, the code tries to pass to the fraction field of the base ring and compute the reversion there.
3. Given a power series over a base ring of positive characteristic, the code tries to lift to a characteristic zero base (using `.lift()`), compute the reversion, and then push back down to the positive characteristic base. This works over finite fields and `Zmod(n)`, but fails over a base ring like `Zmod(n)[x]`.

Each of these is demonstrated with an example.

### comment:12 follow-up:  13 Changed 12 years ago by David Kirkby

Is it sensible to remove the first 3 patches - if so, I can do so, as I have admin rights.

Have all the results been verified by independent means? If so, it would be useful to state how. If not, then I personally don't agree with such "tests", but as a minimum it should be documented if the results have not been verified.

As for the actual maths, I can't review that - this is well outside my level of maths, since none of my degrees are in maths

Dave

### comment:13 in reply to:  12 Changed 12 years ago by Niles Johnson

Replying to drkirkby:

Is it sensible to remove the first 3 patches - if so, I can do so, as I have admin rights.

Sure; thanks :)

Have all the results been verified by independent means? If so, it would be useful to state how. If not, then I personally don't agree with such "tests", but as a minimum it should be documented if the results have not been verified.

The verifications are performed by checking that `f(g)` and `g(f)` both return `x + O(x^prec)`. Actually just one of these is necessary to determine the reversion of f, so I view this as two different methods for verifying the answer.

As for the actual maths, I can't review that - this is well outside my level of maths, since none of my degrees are in maths

Thanks for checking into it :)

### comment:14 in reply to:  11 ; follow-up:  15 Changed 12 years ago by Francis Clarke

Replying to niles:

Yes, especially if there is work in progress to support converting more rings to pari. I wrote a revised patch which first attempts to convert to pari and do reversion there, and then tries the Lagrange inversion if conversion to pari fails. I think that implementation means that this patch can be independent of #4376 -- it will simply perform better when that patch is included.

A good point (but it would be nice if #4376 could be reviewed; it's very short).

I also think it is sensible to be able to revert infinite precision series, either by specifying the desired precision or by using the parent's default precision.

Yes, upon further consideration I agree. I've made this and two other improvements: 1. Given a power series with infinite precision and degree `deg`, its reversion is computed with precision `deg + 1`.

I don't see the logic for this. I would suggest having a keyword `precision` with default `None`, and replacing

```    if f.prec() is infinity:
out_prec = f.degree() + 1
f = f.add_bigoh(out_prec)
else:
out_prec = f.prec()
```

with

```    if f.prec() is infinity and precision is None:
precision = f.parent().default_prec()
if precision:
f = f.add_bigoh(precision)
```

Then one could do (to get some Catalan numbers):

```sage: R.<x> = QQ[[]]
sage: (x - x^2).reversion()
x + x^2 + 2*x^3 + 5*x^4 + 14*x^5 + 42*x^6 + 132*x^7 + 429*x^8
+ 1430*x^9 + 4862*x^10 + 16796*x^11 + 58786*x^12 + 208012*x^13
+ 742900*x^14 + 2674440*x^15 + 9694845*x^16 + 35357670*x^17
+ 129644790*x^18 + 477638700*x^19 + O(x^20)
```

or

```sage: (x - x^2).reversion(precision=8)
x + x^2 + 2*x^3 + 5*x^4 + 14*x^5 + 42*x^6 + 132*x^7 + O(x^8)
```

rather than

```sage: (x - x^2).reversion()
x + x^2 + O(x^3)
```

### Changed 12 years ago by Niles Johnson

Attachment: trac_7644_reversion_lagrange_4.patch​ added

apply only this patch; choice for default precision improved

### comment:15 in reply to:  14 Changed 12 years ago by Niles Johnson

Description: modified (diff)

Replying to fwclarke:

Replying to niles:

A good point (but it would be nice if #4376 could be reviewed; it's very short).

agreed -- I'll add it to my list

I don't see the logic for this. I would suggest ...

that's a better idea; it's implemented in the new patch

### comment:16 follow-up:  17 Changed 12 years ago by Francis Clarke

This is looking good. Some comments on the docstring:

```If f has infinite precision, then the precision
of the reversion defaults to the default precision of
``f.parent()``.
```

should say "unless precision is set", or words to that effect.

The remark a few lines later that

```self must have finite precision (i.e. this cannot be done
for polynomials).
```

is no longer appropriate. Nor is

```Under the current implementation, the leading
coefficient must be a unit in the base ring, and the base ring must
have characteristic zero.
```

I'm not sure that

```In positive characteristic, attempt first to lift to characteristic
zero and perform the reversion there:
```

is in the right place; it's really to do with the algorithm. Perhaps better to say here that the function can handle some base rings of non-zero characteristic.

### comment:17 in reply to:  16 Changed 12 years ago by Niles Johnson

Replying to fwclarke:

This is looking good. Some comments on the docstring:

...

Sorry for not catching these before; fixed now.

### comment:18 Changed 12 years ago by Francis Clarke

Status: needs_review → needs_work

This is fine now, except for one little thing. In the patch file, the line

```+        most ``f.prec()).  If f has infinite precision, and the argument
```

should read

```+        most ``f.prec())``.  If ``f`` has infinite precision, and the argument
```

### comment:19 Changed 12 years ago by Niles Johnson

Description: modified (diff) needs_work → needs_review

Thanks for catching that -- fixed now.

### comment:20 Changed 12 years ago by Francis Clarke

Reviewers: → Francis Clarke needs_review → positive_review

Now it's perfect. Positive review.

### comment:21 Changed 12 years ago by Jeroen Demeyer

Milestone: sage-4.6.1 → sage-4.6.2

### comment:22 Changed 12 years ago by Jeroen Demeyer

Status: positive_review → needs_work

You should change the commit message of the patch (use `hg qrefresh -e` for that). Make sure you include the ticket number.

### Changed 12 years ago by Niles Johnson

Attachment: trac_7644_reversion_lagrange_5.patch​ added

apply only this patch

### comment:23 Changed 12 years ago by Niles Johnson

Status: needs_work → needs_review

done.

### comment:24 Changed 12 years ago by Francis Clarke

Status: needs_review → positive_review

commit message ok now

### comment:25 Changed 12 years ago by Jeroen Demeyer

Status: positive_review → needs_work

There are some formatting errors in the docstrings (a `::` should be followed by indented code, there are some places where this is violated). I can fix these, but later.

### comment:26 Changed 12 years ago by Niles Johnson

My apologies -- I'll upload a fixed version soon.

### Changed 12 years ago by Niles Johnson

Attachment: trac_7644_reversion_lagrange_6.patch​ added

apply only this patch

### comment:27 follow-up:  28 Changed 12 years ago by Niles Johnson

Description: modified (diff) needs_work → needs_review

The documentation issue is fixed now -- all documentation builds without warnings.

### comment:28 in reply to:  27 Changed 12 years ago by Jeroen Demeyer

Status: needs_review → positive_review

Replying to niles:

The documentation issue is fixed now -- all documentation builds without warnings.

Agreed.

### comment:29 Changed 12 years ago by Jeroen Demeyer

Merged in: → sage-4.6.2.alpha1 → fixed positive_review → closed
Note: See TracTickets for help on using tickets.