Inversion and fraction fields for power series rings

[simon-king-jena]

This ticket is about at least three bugs related with inversion of elements of power series rings.

Here is the first:

sage: R.<x> = ZZ[[]]
sage: (1/x).parent()
Laurent Series Ring in x over Integer Ring
sage: (x/x).parent()
Power Series Ring in x over Integer Ring

Both parents are wrong. Usually, the parent of a/b is the fraction field of the parent of a,b, even if a==b. And neither above parent is a field.

Next bug:

sage: (1/(2*x)).parent()
ERROR: An unexpected error occurred while tokenizing input
The following traceback may be corrupted or invalid
The error message is: ('EOF in multi-line statement', (919, 0))

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
... very long traceback
TypeError: no conversion of this rational to integer

And the third:

sage: F = FractionField(R)
sage: 1/x in F
False

With the proposed patch, we solve all bugs except

sage: (x/x).parent()
Power Series Ring in x over Integer Ring

---

Apply (old):

• attachment: trac-8972_fraction_of_power_series_combined.patch

CC: @tscrim @mkoeppe @nbruin @jhpalmieri @dkrenn @bgrenet @sagetrac-tmonteil @videlec @DaveWitteMorris

Component: algebra

Keywords: power series ring, fraction field

Work Issues: doctest failures

Author: Simon King, Michael Jung

Branch/Commit: 5aca068

Reviewer: Robert Bradshaw, Travis Scrimshaw

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

[simon-king-jena]

Work Issues: segfault of div of Laurent series

[simon-king-jena]

comment:1

OK, here is a patch.

Idea:

• Introduce a fraction_field method to power series rings that returns the Laurent series ring over the fraction field of the base ring. I hope that the base_extend method for Laurent polynomials works stably enough, otherwise it should be rewritten.
• In the division method of power serise, construct the fraction field of the parent, move numerator and denominator into the fraction field, and divide there. The original form of the method is preserved as a fail safe.

With the patch, I get:

sage: P.<t> = ZZ[]
sage: R.<x> = P[[]]
sage: FractionField(R)
Laurent Series Ring in x over Fraction Field of Univariate Polynomial Ring in t over Integer Ring

O dear. I just realise that there will be more work. There is a segmentation fault, as follows:

sage: R1.<x> = ZZ[[]]
sage: F = FractionField(R1)
sage: F
Laurent Series Ring in x over Rational Field
sage: F(x)
x
sage: ~F(x)
/home/king/SAGE/sage-4.3.1/local/bin/sage-sage: line 206:  4437 Segmentation fault      sage-ipython "$@" -i

So, the division of elements of a Laurent series ring fails with a segmentation fault. By consequence, division of power series segfaults as well, with the patch. "Needs work", I presume.

[simon-king-jena]

comment:2

Strangely, without the patch, the construction works:

sage: R.<x> = ZZ[[]]
sage: F = LaurentSeriesRing(FractionField(R.base()),R.variable_names())
sage: 1/F(x)
x^-1

The patch does this construction as well -- but segfaults. There seems to be a nasty side effect.

[simon-king-jena]

comment:3

The segfault problem seems to come from the fact that the div method for Laurent series relies on the div method for power series -- and with my patch, the div method for power series uses the div method for Laurent series. Not good. But that seems solvable.

[simon-king-jena]

Attachment: 8972_power_series_inverses.patch.gz

Bugfixes for fraction field and inverses of power series over non-fields

[simon-king-jena]

Changed work issues from segfault of div of Laurent series to none

[simon-king-jena]

comment:4

OK, I replaced the old patch, and now it seems to work!

For example:

sage: P.<t> = ZZ[]
sage: R.<x> = P[[]]
sage: 1/(t*x)
1/t*x^-1
sage: (1/x).parent() is FractionField(R)
True
sage: Frac(R)
Laurent Series Ring in x over Fraction Field of Univariate Polynomial Ring in t over Integer Ring

The doc tests for the two modified files pass. So, ready for review!

[simon-king-jena]

Author: Simon King

[simon-king-jena]

comment:5

PS: Even the following works:

sage: 1/(t+x)
1/t - 1/t^2*x + (-1/-t^3)*x^2 + (1/-t^4)*x^3 + 1/t^5*x^4 - 1/t^6*x^5 + 1/t^7*x^6 - 1/t^8*x^7 + 1/t^9*x^8 + (1/-t^10)*x^9 + 1/t^11*x^10 + (1/-t^12)*x^11 + 1/t^13*x^12 + (1/-t^14)*x^13 + (-1/-t^15)*x^14 + (1/-t^16)*x^15 + (-1/-t^17)*x^16 + (1/-t^18)*x^17 + 1/t^19*x^18 - 1/t^20*x^19 + O(x^20)

[robertwb]

comment:6

I'm wary of such a big change until we have some timing tests in place. In particular, I'm worried about potential slowdowns in the Monsky-Washnitzer code.

[simon-king-jena]

comment:7

Replying to @robertwb:

I'm wary of such a big change until we have some timing tests in place. In particular, I'm worried about potential slowdowns in the Monsky-Washnitzer code.

You are right, timing should be taken into account, and this is something my patch doesn't provide. Without the patch:

sage: R.<x> = ZZ[[]]
sage: timeit('a=1/(1+x)')
625 loops, best of 3: 1.08 ms per loop

With the patch

sage: R.<x> = ZZ[[]]
sage: timeit('a=1/(1+x)')
125 loops, best of 3: 6 ms per loop

So, it is slower by more than a factor five.

[simon-king-jena]

Work Issues: improve timings

[simon-king-jena]

comment:8

Concerning timings, I see a couple of things that might help improve the div method:

1. My patch calls the fraction_field method; in addition the original code calls the laurent_series_ring method. But the purpose of both is the same. So, only one should be done.
2. The fraction_field method should better be cached; this would save a lot of time, since the fradtion field construction involves the construction of a Laurent series ring and the construction of the fraction field of the base ring.
3. Currently the div methods of power series and of Laurent series call each other; I don't know if this is done efficiently (e.g., avoiding Python calls). I could imagine that this can be stratified.

So, something to do for tomorrow. And I guess the ticket is again "needs work".

[simon-king-jena]

comment:9

Replying to @simon-king-jena:

Concerning timings, I see a couple of things that might help improve the div method:

One more thing: The old code is quick if the result actually belongs to the power series ring, which is quite often the case; if this is not the case then often an error results. And I guess the parent should always be the fraction field, eventually.

What I just tested (but I really should get some sleep now...) is to cache the fraction_field method, and to try to use the old code if the valuation of the denominator is not bigger than the valuation of the numerator; if this fails, then put numerator and denominator into the fraction field, and try again.

Doing so brings the above timing to about 2ms, which is still a loss of factor two, but two is less than 5 or 6. I'll submit another patch after trying to lift the dependency of the two div methods on each other.

[simon-king-jena]

Attachment: 8972_power_series_timing.patch.gz

Improving the timings, to be applied after the bug fix patch

[simon-king-jena]

Changed work issues from improve timings to none

[simon-king-jena]

comment:10

The second patch does several things:

1. It turned out that in several places it is assumed that the quotient of two power series is a power series, not a Laurent series. I tried to take care of this.

2. I improved the timings, so that (according to the timings below) the performance is competitive (sometimes clearly better, sometimes very little worse) then the original performance.

The reason why the old timings were good was some kind of lazyness: The result of a division was not in the fraction field (as it should be!), but in a simpler ring, so that subsequent computations became easier. On the other hand, transformation into a Laurent polynomial was quite expensive, since there always was a transformation into the Laurent Series Ring's underlying Power Series Ring -- even if the given data already belong to it.

Solution (as usual): Add an optional argument "check" to the init method of Laurent Series.

And then I tried to benefit from keeping the results as simple as possible (like the old code did), but in a more proper way. Let f be a Laurent series in a Laurent series ring L. In the old code, one always had L.power_series_ring() is f.valuation_zero_part().parent(). I suggest to relax this condition: It is sufficient to have a coercion:

sage: R.<x> = ZZ[[]]
sage: f = 1/(1+2*x)
sage: f.parent()
Laurent Series Ring in x over Rational Field
sage: f.parent().power_series_ring()
Power Series Ring in x over Rational Field
sage: f.power_series().parent()
Power Series Ring in x over Integer Ring

Timings

Without the two patches:

sage: R.<x> = ZZ[[]]
sage: timeit('a=1/x')
625 loops, best of 3: 291 µs per loop
sage: timeit('a=(1+x)/x')
625 loops, best of 3: 295 µs per loop
sage: timeit('a=1/(1+x)')
625 loops, best of 3: 1.07 ms per loop
sage: timeit('a=(1+x)/(1-x)')
625 loops, best of 3: 1.07 ms per loop
sage: y = (3*x+2)/(1+x)
sage: y=y/x
sage: z = (x+x^2+1)/(1+x^4+2*x^3+4*x)
sage: z=y/x
sage: timeit('y+z')
625 loops, best of 3: 213 µs per loop
sage: timeit('y*z')
625 loops, best of 3: 118 µs per loop
sage: timeit('y-z')
625 loops, best of 3: 212 µs per loop
sage: timeit('y/z')
ERROR: An unexpected error occurred while tokenizing input
The following traceback may be corrupted or invalid
The error message is: ('EOF in multi-line statement', (27, 0))

---------------------------------------------------------------------------
ArithmeticError                           Traceback (most recent call last)
...
ArithmeticError: division not defined

With the two patches:

sage: R.<x> = ZZ[[]]
sage: timeit('a=1/x')
625 loops, best of 3: 220 µs per loop
sage: timeit('a=(1+x)/x')
625 loops, best of 3: 228 µs per loop
sage: timeit('a=1/(1+x)')
625 loops, best of 3: 1.25 ms per loop
sage: timeit('a=(1+x)/(1-x)')
625 loops, best of 3: 1.26 ms per loop
sage: y = (3*x+2)/(1+x)
sage: y=y/x
sage: z = (x+x^2+1)/(1+x^4+2*x^3+4*x)
sage: z=y/x
sage: timeit('y+z')
625 loops, best of 3: 191 µs per loop
sage: timeit('y*z')
625 loops, best of 3: 92.6 µs per loop
sage: timeit('y-z')
625 loops, best of 3: 191 µs per loop
sage: timeit('y/z')
25 loops, best of 3: 9.35 ms per loop

So, my patches not only fix some bugs, but they also slightly improve the performance.

I tested all files that I have changed, but I can not exclude errors in other files.

I forgot to insert my name as an author, but presumably there will be a revision needed anyway, after comments of referees...

[simon-king-jena]

comment:11

I forgot one remark:

I don't know how this happens, but the truncate method of Laurent series behaves different than before, although the truncate method of power series did not change:

sage: A.<t> = QQ[[]]
sage: f = (1+t)^100
sage: f.truncate(5)
3921225*t^4 + 161700*t^3 + 4950*t^2 + 100*t + 1
sage: f.truncate(5).parent()
Univariate Polynomial Ring in t over Rational Field
sage: g = 1/f
sage: g.truncate(5)
1 - 100*t + 5050*t^2 - 171700*t^3 + 4421275*t^4 + O(t^5)
sage: g.truncate(5).parent()
Laurent Series Ring in t over Rational Field

In other words, g.truncate(5) is now returning a Laurent series, but without the patch it used to return return a univariate polynomial, similar to f.truncate(5).

I don't know if this is acceptable, and also I don't understand how that happened. Shall I change it?

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

f0334c2	Trac #8972: case of non integral domain
e4b0977	Trac #8972: adapt doctests

[sagetrac-git]

Changed commit from 2602f8d to e4b0977

[sagetrac-git]

Changed commit from e4b0977 to 7399856

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

7399856

Trac #8972: use fraction_field instead of laurent_series

[mjungmath]

comment:68

This should do.

[mjungmath]

Changed author from Simon King to Simon King, Michael Jung

[pjbruin]

comment:70

Replying to @mjungmath:

Approach might be too naive, but it seems to solve the first and the last bug without compromising the timing.

Does this also address the dependency #21283?

[mjungmath]

Changed dependencies from #21283 to none

[mjungmath]

comment:73

Replying to @pjbruin:

Replying to @mjungmath:

Approach might be too naive, but it seems to solve the first and the last bug without compromising the timing.

Does this also address the dependency #21283?

Nope. I just wanted to see what the patchbot has to say before I modify the ticket description. Here we go.

[mjungmath]

comment:75

Yay! Patchbot is green! :)

[mjungmath]

comment:78

Would be nice to have this ticket reviewed. Thanks. :)

[mjungmath]

comment:81

Patchbot seems still be satisfied with beta7. It'd be nice to have this change in version 9.3 though. Otherwise, a rebase might become necessary soon...

[tscrim]

comment:82

There are a few details that still fixing:

This seems strange:

:doc:`laurent_series_ring_element`.

Wouldn't it be better to be a specific reference to the class? Although I am not sure it is worth the direct reference.

This is wrong:

        If we try a non-zero, non-unit leading coefficient, we end up in the in
        the fraction field, i.e. Laurent series ring::

Even for units, you still end up in the fraction field.

-However, this must fail if the underlying ring is no integral domain::
+However, this must fail if the underlying ring is not an integral domain::

[mjungmath]

comment:83

Thank you Travis for taking a look!

Replying to @tscrim:

There are a few details that still fixing:

This seems strange:

:doc:`laurent_series_ring_element`.

Wouldn't it be better to be a specific reference to the class? Although I am not sure it is worth the direct reference.

Alright, I'll change it.

This is wrong:

        If we try a non-zero, non-unit leading coefficient, we end up in the in
        the fraction field, i.e. Laurent series ring::

Even for units, you still end up in the fraction field.

Oh right. I think the previous author already meant constant terms instead of leading terms. At least that would make more sense.

-However, this must fail if the underlying ring is no integral domain::
+However, this must fail if the underlying ring is not an integral domain::

Will do.

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

5aca068

Trac #8972: documentation improved

[sagetrac-git]

Changed commit from 7399856 to 5aca068

[mjungmath]

comment:85

That should be better.

[tscrim]

comment:86

Thanks. LGTM.

[tscrim]

Changed reviewer from Robert Bradshaw to Robert Bradshaw, Travis Scrimshaw

[vbraun]

Changed branch from public/ticket/8972 to 5aca068