Opened 11 years ago
Closed 7 years ago
#3979 closed defect (fixed)
Power series composition messes up precision
Reported by: | kedlaya | Owned by: | somebody |
---|---|---|---|
Priority: | critical | Milestone: | sage-5.4 |
Component: | algebra | Keywords: | power series, composition, precision |
Cc: | chapoton | Merged in: | sage-5.4.rc0 |
Authors: | Francis Clarke | Reviewers: | Kiran Kedlaya, Luis Felipe Tabera Alonso, Frédéric Chapoton |
Report Upstream: | N/A | Work issues: | |
Branch: | Commit: | ||
Dependencies: | Stopgaps: | #12783 |
Description (last modified by )
The composition of two power series is sometimes returned with the wrong precision. A trivial example:
sage: pow.<u> = PowerSeriesRing(Rationals()); print (1 + O(u^4))(u) 1
where the return value should have precision 4 rather than infinity. A more nontrivial example:
sage: pow.<u> = PowerSeriesRing(Rationals()); print (1 + u^2 + O(u^4))(u^2) 1 + u^4 + O(u^10)
where the return value should have precision 8 instead of 10.
Apply
Attachments (8)
Change History (37)
comment:1 follow-up: ↓ 4 Changed 10 years ago by
comment:2 Changed 10 years ago by
- Component changed from calculus to basic arithmetic
- Owner changed from burcin to somebody
Changed 8 years ago by
comment:3 Changed 8 years ago by
- Component changed from basic arithmetic to algebra
- Priority changed from major to critical
- Report Upstream set to N/A
- Status changed from new to needs_review
In the attached patch I have completely rewritten sage.rings.power_series_poly.__call__
. Several errors in the old version have been corrected. The new version more closely follows the corresponding function for polynomials, in particular referring to variables by name is now possible.
In order to make the __call__
function work correctly it was necessary to change the behaviour of sage.rings.power_series_poly.valuation
. At the moment
sage: R.<x> = QQ[] sage: O(x^3).valuation() +Infinity
If we interpret O(x^3)
as x^3
times an unknown power series, then the valuation could be anywhere between 3 and infinity, but 3 is a much better, and more cautious, estimate than infinity. It is also very strange to have a series whose valuation is greater than its precision. The new convention is also consistent with what happens for p-adic integers:
sage: O(7^3).valuation() 3
In the course of checking the power series code, a minor mistake in the polynomial code has been found and corrected.
A doctest in sage/rings/morphism.pyx
needed adjusting.
I have also deleted the doctest in sage.rings.power_series_mpoly.__call__
for two reason's : (1) it doesn't use this function; (2) it makes no sense anyway. Besides the first line of the file is
# NOT ready to be used -- possibly should be deleted.
comment:4 in reply to: ↑ 1 Changed 8 years ago by
comment:5 follow-up: ↓ 6 Changed 8 years ago by
- Status changed from needs_review to needs_work
The patch looks reasonable on its own. However, changing the call syntax for power series generates quite a number of doctest failures elsewhere, by triggering the error message
ValueError: Cannot substitute this value
Here are the examples I found in the rings and schemes directories; there may be more elsewhere that I didn't find. (This used 4.7.1.rc1, but I don't think the version much matters.)
sage -t "devel/sage/sage/rings/multi_power_series_ring.py" sage -t "devel/sage/sage/rings/laurent_series_ring_element.pyx" sage -t "devel/sage/sage/rings/multi_power_series_ring_element.py" sage -t "devel/sage/sage/rings/power_series_ring.py" sage -t "devel/sage/sage/schemes/hyperelliptic_curves/hyperelliptic_generic.py" sage -t "devel/sage/sage/schemes/hyperelliptic_curves/hyperelliptic_padic_field.py" sage -t "devel/sage/sage/schemes/elliptic_curves/ell_wp.py" sage -t "devel/sage/sage/schemes/elliptic_curves/formal_group.py" sage -t "devel/sage/sage/schemes/elliptic_curves/ell_rational_field.py" sage -t "devel/sage/sage/schemes/elliptic_curves/padic_lseries.py" sage -t "devel/sage/sage/schemes/elliptic_curves/padics.py"
This patch can't receive a positive review with all these broken doctests. The best thing would be to fix them all now, but if that is infeasible, I would propose the following.
- Deprecate the old syntax: accept it while raising a
DeprecationWarning
.
- Once this ticket is closed, open a second ticket to modify the syntax in the other doctests.
- Once the second ticket is closed, open a third ticket to remove the old syntax.
comment:6 in reply to: ↑ 5 Changed 8 years ago by
Replying to kedlaya:
Yes, I'm sorry to have missed those failures.
I've been able to fix most of them. The code in schemes/elliptic_curves/formal_group.py
was problematic since iterated univariate power series were used to approximate a power series ring in two variables, but this should be simple to fix with the multiple variable power series available in 4.7.1. More difficult to deal with may be the failures in schemes/hyperelliptic_curves/hyperelliptic_padic_field.py
since the problem here is with substitutions in power series with a known p-adic radius of convergence.
I'll get back to this in a couple of weeks.
comment:7 Changed 8 years ago by
- Description modified (diff)
- Status changed from needs_work to needs_review
I have attached a revised patch. All the previous failures have been dealt with. Some changes were essentially trivial, but more major were:
- Formal groups for elliptic curves have been rewritten to exploit the multi-variable power series code available since 4.7.1
local_coordinates_at_weierstrass
insage/schemes/hyperelliptic_curves/hyperelliptic_generic.py
has been substantially simplified.
- At several places in
sage/schemes/hyperelliptic_curves/hyperelliptic_padic_field.py
substitution in a power series has had to be replaced by substitution in the underlying polynomial. This works for now because in these instances the p-adic radius of convergence is known.
comment:8 follow-ups: ↓ 9 ↓ 12 Changed 8 years ago by
- Reviewers set to Kiran Kedlaya, Luis Felipe Tabera Alonso
- Status changed from needs_review to needs_work
Some problems I have found
This should work
sage: x=polygen(QQ) sage: f = 1 + 3*x + O(x^2) sage: f(x) ... ValueError
This should raise an error:
sage: x = LaurentSeriesRing(QQ,'x').gen() sage: f = x + O(x^2) sage: f(~x) O(x^-2)
You cannot substitute x by 1/x on a power series unless it is a Laurent polynomial.
Suggestions, comments:
On file laurent_series_ring_element.pyx
@446 def laurentpolynomial(...
Improve the documentation. By what is written it seems that the output should be a Laurent polynomial but the method actually returns a Laurent power series.
@1141 call documentation, specify that x needs to have a valuation at least 1.
@1165 raise ValueError?, "must not specify %s keyword and positional argument" % name
Add a doctest to the call method with both keywords and positional arguments one that works (name!= keyword) and one that raises the error, other possibilities welcomed.
On file multi_power_series_ring
@964,989 improve documentation, not clear if the input can be polynomials, powerseries, powerseries + bigoh or in which ring is the result. If we can use big_Oh in the input etc. Maybe for another ticket.
On file local_generic_element.pyx
@140 I would write: Returns self up to reduced precision prec
.
On file polynomial_element.pyx
@461-467 doctest should go in the TESTS section.
@567 raise ValueError?, "must not specify %s keyword and positional argument" % name
On file power_series_mpoly
@74 Add documentation and a valid example to the call method. Each method or function that is modified need to have a correct documentation and doctest.
On file power_series_poly
@290 raise ValueError?, "must not specify %s keyword and positional argument" % name
@936 This is a bug and should be fixed. This is a regression since this works without the patch (except for the incorrect precision).
On file scheme.py
@178 temp2 = temp.exp().change_ring(ZZ)
Is there a reason you want a powerseries in ZZ instead of QQ?
comment:9 in reply to: ↑ 8 Changed 8 years ago by
Replying to lftabera:
Thanks for your careful look at the patch. Most of the problems can be fixed quite easily, though the one at line 936 in power_series_poly
could be more difficult. I hope to submit a new patch soon.
On file power_series_mpoly @74 Add documentation and a valid example to the call method. Each method or function that is modified need to have a correct documentation and doctest.
The reason I removed rather than corrected the doctest in power_series_mpoly
was that it does not test this __call__
function but the one in power_series_poly
(I've seen too many such doctests). In fact I'm not sure that this file is used at all. In fact the first line is
# NOT ready to be used -- possibly should be deleted.
However I didn't have the confidence to delete it myself, and add other issues to an already complicated patch.
Adding documentation to this function would be hard to do since the whole file is so poorly documented that I can't understand what it's for.
comment:10 Changed 7 years ago by
Apologies for hijacking this ticket for this: Francis, what is your email address? F.Clarke@… is bouncing! John (john.cremona@…)
comment:11 Changed 7 years ago by
- Stopgaps set to 12783
comment:12 in reply to: ↑ 8 Changed 7 years ago by
- Description modified (diff)
- Milestone changed from sage-5.0 to sage-5.1
- Status changed from needs_work to needs_review
Replying to lftabera:
Finally I have a revised patch. I have applied all your suggestions with two exceptions. Comments on some of them follow:
On file laurent_series_ring_element.pyx @446
Rather than changing the documentation I have changed the code, so it does now return a Laurent polynomial.
On file multi_power_series_ring @964,989
I've left this for another ticket, as you suggested.
On file power_series_mpoly @74
I've left this unchanged, for reasons explained above.
On file scheme.py @178
I've undone this change. It belongs in another ticket.
comment:13 Changed 7 years ago by
- Cc chapoton added
Note that the patch corrects a bug raised in #12931.
comment:14 Changed 7 years ago by
- Description modified (diff)
Trying to help the bot:
Apply only trac_3979_power_series_substitution_rev2.patch
comment:15 follow-up: ↓ 16 Changed 7 years ago by
- Status changed from needs_review to needs_work
The patch must be rebased on a recent version.
comment:16 in reply to: ↑ 15 Changed 7 years ago by
- Description modified (diff)
- Status changed from needs_work to needs_review
Replying to chapoton:
The patch must be rebased on a recent version.
I've attached new patch.
I hope it can be reviewed before this has to be done again.
comment:17 Changed 7 years ago by
Apply only trac_3979_power_series_substitution_rev3.patch
comment:18 Changed 7 years ago by
- Description modified (diff)
comment:19 Changed 7 years ago by
- Description modified (diff)
comment:20 Changed 7 years ago by
Apply only trac_3979_power_series_substitution_rev4.patch
comment:21 follow-up: ↓ 22 Changed 7 years ago by
I would like to see, when possible, a more specific error instead of
raise ValueError, "Cannot substitute this value"
In particular, when this is because of negative valuation, one should say it.
comment:22 in reply to: ↑ 21 Changed 7 years ago by
- Description modified (diff)
Replying to chapoton:
I would like to see, when possible, a more specific error instead of
raise ValueError, "Cannot substitute this value"In particular, when this is because of negative valuation, one should say it.
A good point. The new patch (to be applied after trac_3979_power_series_substitution_rev4.patch) gives a more explicit error message.
Changed 7 years ago by
comment:23 Changed 7 years ago by
- Description modified (diff)
- Status changed from needs_review to positive_review
This looks good to me. The patches applies smoothly on 5.4beta1. All tests pass. This ticket solves some embarassing problems and is much wanted. Positive review !
comment:24 Changed 7 years ago by
- Reviewers changed from Kiran Kedlaya, Luis Felipe Tabera Alonso to Kiran Kedlaya, Luis Felipe Tabera Alonso, Frédéric Chapoton
comment:25 follow-up: ↓ 26 Changed 7 years ago by
The bot is back and is unhappy because the patch removes one test in rings/power_series_mpoly.pyx
This should be easy to correct, if really required to close the ticket.
comment:26 in reply to: ↑ 25 Changed 7 years ago by
Replying to chapoton:
The bot is back and is unhappy because the patch removes one test in rings/power_series_mpoly.pyx
This should be easy to correct, if really required to close the ticket.
As I said in July last year:
I have also deleted the doctest in
sage.rings.power_series_mpoly.__call__
for two reasons : (1) it doesn't use this function; (2) it makes no sense anyway. Besides the first line of the file is# NOT ready to be used -- possibly should be deleted.
And I explained this more fully a year ago:
The reason I removed rather than corrected the doctest in
power_series_mpoly
was that it does not test this__call__
function but the one inpower_series_poly
(I've seen too many such doctests). In fact I'm not sure that this file is used at all. In fact the first line is# NOT ready to be used -- possibly should be deleted.However I didn't have the confidence to delete it myself, and add other issues to an already complicated patch.
Adding documentation to this function would be hard to do since the whole file is so poorly documented that I can't understand what it's for.
I have now understood how to create an element of the relevant type (something that isn't done anywhere else, as far as I can see):
sage: S.<x> = QQ[] sage: R = sage.rings.power_series_ring.PowerSeriesRing_generic(S, 't', use_lazy_mpoly_ring=True) sage: t = R.gen() sage: f = 3 - x*t^3 + O(t^5) sage: type(f) <type 'sage.rings.power_series_mpoly.PowerSeries_mpoly'> sage: f(2) -2*t^3 + 3 sage: f(2, t^2) 3 - 2*t^6
The final answers are wrong (they shouldn't have infinite precision), inconsistent (the first is a polynomial, the second a power series), and the syntax is non-standard. Compare
sage: T.<u> = S[[]] sage: g = 3 - x*u^3 + O(u^5) sage: g(u^2, 2) 3 - 2*u^6 + O(u^10) sage: g(u^2) 3 - x*u^6 + O(u^10)
Of course g(2)
raises an error.
However, I have made a supplementary patch which reinserts into power_series_mpoly.pyx
a doctest based on the above. This should appease the patchbot, and though problematic it is less bad than before. The alternative would be to rewrite much of the code in power_series_mpoly.pyx
, which would seem to be a waste of time if it is destined for deletion.
comment:27 Changed 7 years ago by
- Description modified (diff)
comment:28 Changed 7 years ago by
- Description modified (diff)
The patchbot tried (and failed) to apply only patches 2 and 3 out of three
So I have merged them all into one patch. Hope this works.
comment:29 Changed 7 years ago by
- Merged in set to sage-5.4.rc0
- Resolution set to fixed
- Status changed from positive_review to closed
- Stopgaps changed from 12783 to #12783
A closely related issue is #5075.