Opened 9 years ago
Closed 8 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 )
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
Attachments (4)
Change History (33)
comment:1 Changed 9 years ago by
- Description modified (diff)
comment:2 follow-up: ↓ 7 Changed 9 years ago by
comment:3 Changed 9 years ago by
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 9 years ago by
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 9 years ago by
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
- 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
- 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 8 years ago by
- Description modified (diff)
comment:9 Changed 8 years ago by
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 8 years ago by
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 8 years ago by
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 8 years ago by
- 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:
- Given a power series with infinite precision and degree
deg
, it's reversion is computed with precisiondeg + 1
. - 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.
- 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 andZmod(n)
, but fails over a base ring likeZmod(n)[x]
.
Each of these is demonstrated with an example.
comment:12 follow-up: ↓ 13 Changed 8 years ago by
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 8 years ago by
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 8 years ago by
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 precisiondeg + 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)
comment:15 in reply to: ↑ 14 Changed 8 years ago by
- Description modified (diff)
comment:16 follow-up: ↓ 17 Changed 8 years ago by
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 8 years ago by
Replying to fwclarke:
This is looking good. Some comments on the docstring:
...
Sorry for not catching these before; fixed now.
comment:18 Changed 8 years ago by
- 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 8 years ago by
- Description modified (diff)
- Status changed from needs_work to needs_review
Thanks for catching that -- fixed now.
comment:20 Changed 8 years ago by
- Reviewers set to Francis Clarke
- Status changed from needs_review to positive_review
Now it's perfect. Positive review.
comment:21 Changed 8 years ago by
- Milestone changed from sage-4.6.1 to sage-4.6.2
comment:22 Changed 8 years ago by
- 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.
comment:24 Changed 8 years ago by
- Status changed from needs_review to positive_review
commit message ok now
comment:25 Changed 8 years ago by
- 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 8 years ago by
My apologies -- I'll upload a fixed version soon.
comment:27 follow-up: ↓ 28 Changed 8 years ago by
- 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 8 years ago by
- 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 8 years ago by
- Merged in set to sage-4.6.2.alpha1
- Resolution set to fixed
- Status changed from positive_review to closed
Lagrange inversion is significantly faster. With
I found