# Ticket #3979(closed defect: fixed)

Opened 5 years ago

Last modified 8 months ago

## Power series composition messes up precision

Reported by: Owned by: kedlaya somebody critical sage-5.4 algebra power series, composition, precision chapoton N/A Kiran Kedlaya, Luis Felipe Tabera Alonso, Frédéric Chapoton Francis Clarke sage-5.4.rc0 #12783

### Description (last modified by fwclarke) (diff)

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

## Change History

### comment:1 follow-up: ↓ 4 Changed 4 years ago by kedlaya

A closely related issue is #5075.

### comment:2 Changed 4 years ago by burcin

• Owner changed from burcin to somebody
• Component changed from calculus to basic arithmetic

### comment:3 Changed 22 months ago by fwclarke

• Priority changed from major to critical
• Status changed from new to needs_review
• Report Upstream set to N/A
• Component changed from basic arithmetic to algebra
• Authors set to Francis Clarke

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 22 months ago by fwclarke

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.

### comment:5 follow-up: ↓ 6 Changed 22 months ago by kedlaya

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

1. Deprecate the old syntax: accept it while raising a DeprecationWarning.
1. Once this ticket is closed, open a second ticket to modify the syntax in the other doctests.
1. Once the second ticket is closed, open a third ticket to remove the old syntax.

### comment:6 in reply to: ↑ 5 Changed 21 months ago by fwclarke

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.

### Changed 20 months ago by fwclarke

Apply only this file

### comment:7 Changed 20 months ago by fwclarke

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

I have attached a revised patch. All the previous failures have been dealt with. Some changes were essentially trivial, but more major were:

1. Formal groups for elliptic curves have been rewritten to exploit the multi-variable power series code available since 4.7.1
1. local_coordinates_at_weierstrass in sage/schemes/hyperelliptic_curves/hyperelliptic_generic.py has been substantially simplified.
1. 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 20 months ago by lftabera

• Status changed from needs_review to needs_work
• Reviewers set to Kiran Kedlaya, Luis Felipe Tabera Alonso

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 20 months ago by fwclarke

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 19 months ago by cremona

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

### comment:11 Changed 14 months ago by kedlaya

• Stopgaps set to 12783

### Changed 12 months ago by fwclarke

Replaces previous revised patch

### comment:12 in reply to: ↑ 8 Changed 12 months ago by fwclarke

• Status changed from needs_work to needs_review
• Description modified (diff)
• Milestone changed from sage-5.0 to sage-5.1

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 12 months ago by fwclarke

• Cc chapoton added

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

### comment:14 Changed 9 months ago by chapoton

• Description modified (diff)

Trying to help the bot:

Apply only trac_3979_power_series_substitution_rev2.patch

### comment:15 follow-up: ↓ 16 Changed 9 months ago by chapoton

• Status changed from needs_review to needs_work

The patch must be rebased on a recent version.

Rebased for 5.2

### comment:16 in reply to: ↑ 15 Changed 9 months ago by fwclarke

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

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 9 months ago by chapoton

Apply only trac_3979_power_series_substitution_rev3.patch

### comment:18 Changed 9 months ago by chapoton

• Description modified (diff)

### comment:19 Changed 9 months ago by chapoton

• Description modified (diff)

### comment:20 Changed 9 months ago by chapoton

Apply only trac_3979_power_series_substitution_rev4.patch

### comment:21 follow-up: ↓ 22 Changed 9 months ago by 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.

### Changed 9 months ago by fwclarke

apply after trac_3979_power_series_substitution_rev4.patch

### comment:22 in reply to: ↑ 21 Changed 9 months ago by fwclarke

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

### comment:23 Changed 8 months ago by chapoton

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

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 8 months ago by chapoton

• 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 8 months ago by 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.

### comment:26 in reply to: ↑ 25 Changed 8 months ago by fwclarke

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

### Changed 8 months ago by fwclarke

Apply after trac_3979_power_series_substitution_rev4_extra.patch

### comment:27 Changed 8 months ago by fwclarke

• Description modified (diff)

### Changed 8 months ago by fwclarke

Apply only this patch

### comment:28 Changed 8 months ago by fwclarke

• 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 8 months ago by jdemeyer

• Status changed from positive_review to closed
• Resolution set to fixed
• Merged in set to sage-5.4.rc0
• Stopgaps changed from 12783 to #12783
Note: See TracTickets for help on using tickets.