Power series composition messes up precision

kedlaya commented 16 years ago

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

attachment: trac_3979_power_series_substitution_rev5.patch

CC: @fchapoton

Component: algebra

Keywords: power series, composition, precision

Stopgaps: #12783

Author: Francis Clarke

Reviewer: Kiran Kedlaya, Luis Felipe Tabera Alonso, Frédéric Chapoton

Merged: sage-5.4.rc0

Issue created by migration from https://trac.sagemath.org/ticket/3979

kedlaya commented 15 years ago

comment:1

A closely related issue is #5075.

5ef57b06-33a8-4a32-80d2-0f5ab07e80d9 commented 13 years ago

comment:3

Attachment: trac_3979_power_series_substitution.patch.gz

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.

5ef57b06-33a8-4a32-80d2-0f5ab07e80d9 commented 13 years ago

Author: Francis Clarke

5ef57b06-33a8-4a32-80d2-0f5ab07e80d9 commented 13 years ago

comment:4

Replying to @kedlaya:

A closely related issue is #5075.

Related, but I don't believe it's the same. The problem in #5075 is still there after the patch.

kedlaya commented 13 years ago

comment:5

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.

5ef57b06-33a8-4a32-80d2-0f5ab07e80d9 commented 13 years ago

comment:6

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.

5ef57b06-33a8-4a32-80d2-0f5ab07e80d9 commented 13 years ago

Apply only this file

5ef57b06-33a8-4a32-80d2-0f5ab07e80d9 commented 13 years ago

comment:7

Attachment: trac_3979_power_series_substitution_rev1.patch.gz

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 in sage/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.

5ef57b06-33a8-4a32-80d2-0f5ab07e80d9 commented 13 years ago

Description changed:

--- 
+++ 
@@ -12,3 +12,6 @@

where the return value should have precision 8 instead of 10.

+--- + +Apply trac_3979_power_series_substitution_rev1.patch

lftabera commented 13 years ago

Reviewer: Kiran Kedlaya, Luis Felipe Tabera Alonso

lftabera commented 13 years ago

comment:8

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?

5ef57b06-33a8-4a32-80d2-0f5ab07e80d9 commented 13 years ago

comment:9

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.

JohnCremona commented 13 years ago

comment:10

Apologies for hijacking this ticket for this: Francis, what is your email address? F.Clarke@swansea.ac.uk is bouncing! John (john.cremona@gmail.com)

kedlaya commented 12 years ago

Stopgaps: 12783

5ef57b06-33a8-4a32-80d2-0f5ab07e80d9 commented 12 years ago

Replaces previous revised patch

5ef57b06-33a8-4a32-80d2-0f5ab07e80d9 commented 12 years ago

Description changed:

--- 
+++ 
@@ -14,4 +14,4 @@

 ---

-Apply [trac_3979_power_series_substitution_rev1.patch](https://github.com/sagemath/sage-prod/files/10641829/trac_3979_power_series_substitution_rev1.patch.gz)
+Apply [trac_3979_power_series_substitution_rev2.patch](https://github.com/sagemath/sage-prod/files/10641830/trac_3979_power_series_substitution_rev2.patch.gz)

5ef57b06-33a8-4a32-80d2-0f5ab07e80d9 commented 12 years ago

comment:12

Attachment: trac_3979_power_series_substitution_rev2.patch.gz

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.

5ef57b06-33a8-4a32-80d2-0f5ab07e80d9 commented 12 years ago

comment:13

Note that the patch corrects a bug raised in #12931.

fchapoton commented 12 years ago

Description changed:

--- 
+++ 
@@ -14,4 +14,4 @@

 ---

-Apply [trac_3979_power_series_substitution_rev2.patch](https://github.com/sagemath/sage-prod/files/10641830/trac_3979_power_series_substitution_rev2.patch.gz)
+Apply only trac_3979_power_series_substitution_rev2.patch

fchapoton commented 12 years ago

comment:14

Trying to help the bot:

Apply only trac_3979_power_series_substitution_rev2.patch

fchapoton commented 12 years ago

comment:15

The patch must be rebased on a recent version.

5ef57b06-33a8-4a32-80d2-0f5ab07e80d9 commented 12 years ago

Rebased for 5.2

5ef57b06-33a8-4a32-80d2-0f5ab07e80d9 commented 12 years ago

Description changed:

--- 
+++ 
@@ -14,4 +14,4 @@

 ---

-Apply only trac_3979_power_series_substitution_rev2.patch
+Apply only trac_3979_power_series_substitution_rev3.patch

5ef57b06-33a8-4a32-80d2-0f5ab07e80d9 commented 12 years ago

comment:16

Attachment: trac_3979_power_series_substitution_rev3.patch.gz

Replying to @fchapoton:

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.

fchapoton commented 12 years ago

comment:17

Apply only trac_3979_power_series_substitution_rev3.patch

fchapoton commented 12 years ago

Description changed:

--- 
+++ 
@@ -14,4 +14,4 @@

 ---

-Apply only trac_3979_power_series_substitution_rev3.patch
+Apply only [attachment: trac_3979_power_series_substitution_rev3.patch](https://github.com/sagemath/sage-prod/files/10641831/trac_3979_power_series_substitution_rev3.patch.gz)

fchapoton commented 12 years ago

Description changed:

--- 
+++ 
@@ -14,4 +14,4 @@

 ---

-Apply only [attachment: trac_3979_power_series_substitution_rev3.patch](https://github.com/sagemath/sage-prod/files/10641831/trac_3979_power_series_substitution_rev3.patch.gz)
+Apply only [attachment: trac_3979_power_series_substitution_rev4.patch](https://github.com/sagemath/sage-prod/files/10641833/trac_3979_power_series_substitution_rev4.patch.gz)

fchapoton commented 12 years ago

comment:20

Apply only trac_3979_power_series_substitution_rev4.patch

fchapoton commented 12 years ago

comment:21

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.

5ef57b06-33a8-4a32-80d2-0f5ab07e80d9 commented 12 years ago

apply after trac_3979_power_series_substitution_rev4.patch

5ef57b06-33a8-4a32-80d2-0f5ab07e80d9 commented 12 years ago

Description changed:

--- 
+++ 
@@ -14,4 +14,4 @@

 ---

-Apply only [attachment: trac_3979_power_series_substitution_rev4.patch](https://github.com/sagemath/sage-prod/files/10641833/trac_3979_power_series_substitution_rev4.patch.gz)
+Apply [attachment: trac_3979_power_series_substitution_rev4.patch](https://github.com/sagemath/sage-prod/files/10641833/trac_3979_power_series_substitution_rev4.patch.gz) followed by  [attachment: trac_3979_power_series_substitution_rev4_extra.patch](https://github.com/sagemath/sage-prod/files/10641832/trac_3979_power_series_substitution_rev4_extra.patch.gz)

5ef57b06-33a8-4a32-80d2-0f5ab07e80d9 commented 12 years ago

comment:22

Attachment: trac_3979_power_series_substitution_rev4_extra.patch.gz

Replying to @fchapoton:

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.

fchapoton commented 12 years ago

Attachment: trac_3979_power_series_substitution_rev4.patch.gz

fchapoton commented 12 years ago

comment:23

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 !

fchapoton commented 12 years ago

Description changed:

--- 
+++ 
@@ -14,4 +14,6 @@

 ---

-Apply [attachment: trac_3979_power_series_substitution_rev4.patch](https://github.com/sagemath/sage-prod/files/10641833/trac_3979_power_series_substitution_rev4.patch.gz) followed by  [attachment: trac_3979_power_series_substitution_rev4_extra.patch](https://github.com/sagemath/sage-prod/files/10641832/trac_3979_power_series_substitution_rev4_extra.patch.gz)
+Apply 
+* [attachment: trac_3979_power_series_substitution_rev4.patch](https://github.com/sagemath/sage-prod/files/10641833/trac_3979_power_series_substitution_rev4.patch.gz)
+* [attachment: trac_3979_power_series_substitution_rev4_extra.patch](https://github.com/sagemath/sage-prod/files/10641832/trac_3979_power_series_substitution_rev4_extra.patch.gz)

fchapoton commented 12 years ago

Changed reviewer from Kiran Kedlaya, Luis Felipe Tabera Alonso to Kiran Kedlaya, Luis Felipe Tabera Alonso, Frédéric Chapoton

fchapoton commented 12 years ago

comment:25

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.

5ef57b06-33a8-4a32-80d2-0f5ab07e80d9 commented 12 years ago

comment:26

Replying to @fchapoton:

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

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.

5ef57b06-33a8-4a32-80d2-0f5ab07e80d9 commented 12 years ago

Apply after trac_3979_power_series_substitution_rev4_extra.patch

5ef57b06-33a8-4a32-80d2-0f5ab07e80d9 commented 12 years ago

comment:27

Attachment: trac_3979_power_series_substitution_rev4_supplementary.patch.gz

5ef57b06-33a8-4a32-80d2-0f5ab07e80d9 commented 12 years ago

Description changed:

--- 
+++ 
@@ -17,3 +17,4 @@
 Apply 
 * [attachment: trac_3979_power_series_substitution_rev4.patch](https://github.com/sagemath/sage-prod/files/10641833/trac_3979_power_series_substitution_rev4.patch.gz)
 * [attachment: trac_3979_power_series_substitution_rev4_extra.patch](https://github.com/sagemath/sage-prod/files/10641832/trac_3979_power_series_substitution_rev4_extra.patch.gz)
+* [attachment: trac_3979_power_series_substitution_rev4_supplementary.patch](https://github.com/sagemath/sage-prod/files/10641834/trac_3979_power_series_substitution_rev4_supplementary.patch.gz)

5ef57b06-33a8-4a32-80d2-0f5ab07e80d9 commented 12 years ago

Apply only this patch

5ef57b06-33a8-4a32-80d2-0f5ab07e80d9 commented 12 years ago

comment:28

Attachment: trac_3979_power_series_substitution_rev5.patch.gz

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.

5ef57b06-33a8-4a32-80d2-0f5ab07e80d9 commented 12 years ago

Description changed:

--- 
+++ 
@@ -15,6 +15,5 @@
 ---

 Apply 
-* [attachment: trac_3979_power_series_substitution_rev4.patch](https://github.com/sagemath/sage-prod/files/10641833/trac_3979_power_series_substitution_rev4.patch.gz)
-* [attachment: trac_3979_power_series_substitution_rev4_extra.patch](https://github.com/sagemath/sage-prod/files/10641832/trac_3979_power_series_substitution_rev4_extra.patch.gz)
-* [attachment: trac_3979_power_series_substitution_rev4_supplementary.patch](https://github.com/sagemath/sage-prod/files/10641834/trac_3979_power_series_substitution_rev4_supplementary.patch.gz)
+* [attachment: trac_3979_power_series_substitution_rev5.patch](https://github.com/sagemath/sage-prod/files/10641835/trac_3979_power_series_substitution_rev5.patch.gz)
+

jdemeyer commented 12 years ago

Merged: sage-5.4.rc0

jdemeyer commented 12 years ago

Changed stopgaps from 12783 to #12783

sagemath / sage

Power series composition messes up precision #3979