Opened 10 years ago

Closed 9 years ago

#7644 closed enhancement (fixed)

generic power series reversion

Reported by: was Owned by: AlexGhitza
Priority: major Milestone: sage-4.6.2
Component: algebra Keywords: lagrange, reversion
Cc: fwclarke Merged in: sage-4.6.2.alpha1
Authors: Niles Johnson Reviewers: Francis Clarke
Report Upstream: N/A Work issues:
Branch: Commit:
Dependencies: Stopgaps:

Description (last modified by niles)

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

Attachments (4)

trac_7644_reversion_lagrange_3.patch (17.3 KB) - added by niles 9 years ago.
apply only this patch; tries to use pari first, then falls back to Lagrange inversion
trac_7644_reversion_lagrange_4.patch (18.6 KB) - added by niles 9 years ago.
apply only this patch; choice for default precision improved
trac_7644_reversion_lagrange_5.patch (18.9 KB) - added by niles 9 years ago.
apply only this patch
trac_7644_reversion_lagrange_6.patch (18.8 KB) - added by niles 9 years ago.
apply only this patch

Download all attachments as: .zip

Change History (33)

comment:1 Changed 10 years ago by mvngu

  • Description modified (diff)

comment:2 follow-up: Changed 10 years ago by fwclarke

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 10 years ago by was

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: Changed 10 years ago by 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?

comment:5 in reply to: ↑ 4 Changed 10 years ago by fwclarke

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 9 years ago by niles

  • Status changed from new to 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 9 years ago by niles

  • Authors set to niles
  • Cc fwclarke added
  • Keywords 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 9 years ago by niles

  • Authors changed from niles to Niles Johnson
  • Description modified (diff)

comment:9 Changed 9 years ago by niles

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: Changed 9 years ago by fwclarke

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 9 years ago by niles

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

comment:11 in reply to: ↑ 10 ; follow-up: Changed 9 years ago by niles

  • 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: Changed 9 years ago by drkirkby

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 9 years ago by niles

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: Changed 9 years ago by fwclarke

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 9 years ago by niles

apply only this patch; choice for default precision improved

comment:15 in reply to: ↑ 14 Changed 9 years ago by niles

  • 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: Changed 9 years ago by fwclarke

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 9 years ago by niles

Replying to fwclarke:

This is looking good. Some comments on the docstring:

...

Sorry for not catching these before; fixed now.

comment:18 Changed 9 years ago by fwclarke

  • Status changed from needs_review to 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 9 years ago by niles

  • Description modified (diff)
  • Status changed from needs_work to needs_review

Thanks for catching that -- fixed now.

comment:20 Changed 9 years ago by fwclarke

  • Reviewers set to Francis Clarke
  • Status changed from needs_review to positive_review

Now it's perfect. Positive review.

comment:21 Changed 9 years ago by jdemeyer

  • Milestone changed from sage-4.6.1 to sage-4.6.2

comment:22 Changed 9 years ago by jdemeyer

  • Status changed from positive_review to 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 9 years ago by niles

apply only this patch

comment:23 Changed 9 years ago by niles

  • Status changed from needs_work to needs_review

done.

comment:24 Changed 9 years ago by fwclarke

  • Status changed from needs_review to positive_review

commit message ok now

comment:25 Changed 9 years ago by jdemeyer

  • Status changed from positive_review to 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 9 years ago by niles

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

Changed 9 years ago by niles

apply only this patch

comment:27 follow-up: Changed 9 years ago by niles

  • Description modified (diff)
  • Status changed from needs_work to needs_review

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

comment:28 in reply to: ↑ 27 Changed 9 years ago by jdemeyer

  • Status changed from needs_review to positive_review

Replying to niles:

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

Agreed.

comment:29 Changed 9 years ago by jdemeyer

  • Merged in set to sage-4.6.2.alpha1
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.