implement multivariate truncated power series arithmetic

[williamstein]

Multivariate truncated power series arithmetic has been requested a lot.

Apply:

1. attachment: trac_1956_combined_6.patch

This is a single patch combining all of the partial patches now attached to this ticket.

CC: @mwhansen @simon-king-jena @malb @unzvfu

Component: commutative algebra

Keywords: multivariate power series

Author: Niles Johnson

Reviewer: Martin Albrecht, Simon King

Merged: sage-4.7.1.alpha2

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

[sagetrac-mabshoff]

comment:2

Hi Mike,

hasn't something happened in that direction already or am I confused?

Cheers,

Michael

[nilesjohnson]

Author: niles

[nilesjohnson]

comment:3

I just attached a first version of multivariate power series. In order to pass all tests, one also needs the patches at #9457 and #9443. Without these, everything still works though.

The implementation is based on an idea mentioned in this sage-devel thread:
use univariate power series over multivariate polynomial ring to do arithmetic and track total-degree precision, but display as multivariate power series.

Although there are some limitations of this approach, I think there's a useful amount of functionality here and I'm happy with how this first version works. Having said that, here are some of the issues that are still unresolved. I'd be happy to hear recommendations about how to prioritize them, or others :)

• sparse: I tried to include sparse=True in the right places, but I haven't really tested to see whether it's actually better.

• inherited functions: some of them (e.g. 'variable') don't make sense; is there a way to block inheritance of certain attributes? (del or delattr?) For now, they just return "NotImplementedError"s

• inherited functions 2: for multi_power_series_ring_element, there is a pile of functions that could be implemented but aren't (e.g. sqrt, derivative, generating functions, etc). They're noted in the documentation and lumped together in the code.

• (rich) comparison: comparisons seem to work, but someone who understands the past and future of sage/python comparisons should take a look at it. (The classes I'm inheriting from seem to use an older system.)

• _l_action_ for multivariate power series: the function works, but doesn't get called when it should; I don't understand how to fix that. (This only affects cases where the scalar isn't already in the base ring, like when you want to multiply a power series over ZZ by a rational and have the result automatically move to a power series over QQ.)

• division of multivariate power series: works fine when the denominator is a unit, but not in other cases (even when the result would not be a Laurent series).

• un/pickling: I don't understand this so I haven't done anything with it.

• cython where appropriate: I don't know C or python, so I've ignored this.

• unifying multivariate and univariate power series, as has been done for polynomials: seems like a project for another time.

[nilesjohnson]

comment:4

Uploaded a revised patch which now includes doctests for all* functions. And yes, I did find a fix a number of bugs while doing that. So thanks, whoever decided to require 100% doctest coverage :)

• note: multi_power_series_ring.py has 100% doctest coverage,
  multi_power_series_ring_element.py has 61% doctest coverage, but the functions with missing doctests are all inherited functions that (as described above) raise "NotImplementedError"s.

further note: I commented out the function _im_gens_ which was supposed to be based on the same function for multivariate polynomial rings. I couldn't get it to work even for the polynomial case, so I let it go (after little effort, I admit).

[malb]

comment:5

Hi, I took a quick look, here are some examples:

• You don't need #random for random output, the seed of the pseudo random number generator is reset for each doctest plot
• It would be nice if you could add line breaks around 80 for the doc strings
• the arithmetic functions doctests do not test correctness of the output

I know very little about multivariate power series, thus I'm clueless about whether your implementation strategy is (a) standard and/or (b) clever. Is it something that is usually done?

[nilesjohnson]

comment:6

Replying to @malb:

Thanks!

• You don't need #random for random output, the seed of the pseudo random number generator is reset for each doctest plot

fixed now.

• It would be nice if you could add line breaks around 80 for the doc strings

sorry, not sure what you mean by 80 here; I'm happy to add linebreaks though :)

• the arithmetic functions doctests do not test correctness of the output

I've added comparisons with polynomial arithmetic; do you think I need to also have the output printed?

I know very little about multivariate power series, thus I'm clueless about whether your implementation strategy is (a) standard and/or (b) clever. Is it something that is usually done?

I don't know much about this either, but it seems like a natural enough idea. After a quick look around google, here's one example:

http://algo.inria.fr/seminars/sem00-01/lecerf.html

I'll take a look for more complete references too.

[malb]

comment:7

Replying to @nilesjohnson:

• It would be nice if you could add line breaks around 80 for the doc strings

sorry, not sure what you mean by 80 here; I'm happy to add linebreaks though :)

Sorry, I meant around 80 characters, which is the standard UNIX terminal width.

• the arithmetic functions doctests do not test correctness of the output

I've added comparisons with polynomial arithmetic; do you think I need to also have the output printed?

I'd prefer that to be honest.

I know very little about multivariate power series, thus I'm clueless about whether your implementation strategy is (a) standard and/or (b) clever. Is it something that is usually done?

I don't know much about this either, but it seems like a natural enough idea. After a quick look around google, here's one example: http://algo.inria.fr/seminars/sem00-01/lecerf.html I'll take a look for more complete references too.

Have you compared the speed with, say, Magma?

[nilesjohnson]

comment:8

Replying to @malb:

Replying to @nilesjohnson:

• It would be nice if you could add line breaks around 80 for the doc strings

fixed now

• the arithmetic functions doctests do not test correctness of the output

arithmetic tests now printed, and confirmed with polynomial tests

Have you compared the speed with, say, Magma?

Here are five simple comparisons, three with 2 variables and two with 4 variables, all over QQ. This sage patch performs substantially better in 3/4 cases, with magma performing substantially better in the other case. Using polynomial multiplication easily beats magma's restriction on the number of coefficients it will compute for a lazy power series (test 3).

Processor: 6-core: model name: Intel(R) Xeon(R) CPU X7460 @ 2.66GHz

sage:

sage: L.<a,b,c,d> = MPowerSeriesRing(QQ,4)
sage: %timeit s = (1 + 2*a + 3*b + 4*d + L(0).O(16))^-5
5 loops, best of 3: 3.18 s per loop

sage: time (1 + a^3 + b^3 + c^4 + d^4 + L(0).O(15))^-20
CPU times: user 0.02 s, sys: 0.00 s, total: 0.02 s
Wall time: 0.02 s
<SNIP>

sage: f = -1/6*b^6*d^14 - 1/9*a^4*b^10*c^4*d^12 + b^10*c^11*d^21 - a*b^14*c^24*d^4 - 1/66*a^16*b^11*c^17 + L(0).O(50)
sage: time f^20
CPU times: user 0.16 s, sys: 0.00 s, total: 0.16 s
Wall time: 0.16 s
1/3656158440062976*b^120*d^280 
+ 5/1371059415023616*a^4*b^124*c^4*d^278 
+ 95/4113178245070848*a^8*b^128*c^8*d^276 
- 5/152339935002624*b^124*c^11*d^287 
+ 5/152339935002624*a*b^128*c^24*d^270 
+ 5/10054435710173184*a^16*b^125*c^17*d^266 
+ O(a, b, c, d)^430

magma:

> L<a, b, c, d> := LazyPowerSeriesRing(Rationals(), 4);
> s := (1 + 2*a + 3*b + 4*d)^-5;
> time Coefficients(s, 15);
<SNIP>
Time: 29.020

> s := (1 + a^3 + b^3 + c^4 + d^4)^-20;
> time Coefficients(s, 14);
<SNIP>
Time: 48.850

> s := (-1/6*b^6*d^14 - 1/9*a^4*b^10*c^4*d^12 + b^10*c^11*d^21 - a*b^14*c^24*d^4 - 1/66*a^16*b^11*c^17)^20;
> time Coefficients(s, 14); 
<SNIP> 
Time: 16.340
> time Coefficients(s, 429);                                                                               

>> time Coefficients(s, 429);
                    ^
Runtime error in 'Coefficients': The number of coefficients to be returned must be small

sage:

sage: Z.<a,b> = MPowerSeriesRing(QQ,2)
sage: time (1 + 1/2*a + 3*b + 4*a*b + 1/5*a^2 + 5/6*b^2 + Z(0).O(30))^-20
CPU times: user 12.67 s, sys: 0.01 s, total: 12.68 s
Wall time: 12.68 s
<SNIP>

sage: f = 1 + a^14*b^9 + 8/3*a^22*b^11 - 1/13*a^23*b^19 + a^15*b^28 + Z(0).O(51)
sage: time f^-20
CPU times: user 0.05 s, sys: 0.00 s, total: 0.05 s
Wall time: 0.05 s
1 - 20*a^14*b^9 - 160/3*a^22*b^11 + 20/13*a^23*b^19 
- 20*a^15*b^28 + 210*a^28*b^18 + O(a, b)^51

magma:

> Z<a, b> := LazyPowerSeriesRing(Rationals(), 2);
> t := (1 + 1/2*a + 3*b + 4*a*b + 1/5*a^2 + 5/6*b^2)^-20;
> time Coefficients(t, 29);  
<SNIP>                            
Time: 1.580

> s := (1 + a^14*b^9 + 8/3*a^22*b^11 - 1/13*a^23*b^19 + a^15*b^28)^-20;
> time Coefficients(s, 50);
<SNIP>
Time: 90.840

In the case where this package performs badly, it seems to be caused by the rational coefficients:

sage: Z.<a,b> = MPowerSeriesRing(QQ,2)

sage: time (1 + a + b + a*b + a^2 + b^2 + Z(0).O(30))^-20;
CPU times: user 1.90 s, sys: 0.00 s, total: 1.90 s
Wall time: 1.90 s

sage: time (1 + a + 2*b - a*b + 3*a^2 - b^2 + Z(0).O(30))^-20;
CPU times: user 2.78 s, sys: 0.00 s, total: 2.78 s
Wall time: 2.78 s

sage: time (1 + a + 12*b - 10*a*b + 13*a^2 - 7*b^2 + Z(0).O(30))^-20;
CPU times: user 3.14 s, sys: 0.00 s, total: 3.14 s
Wall time: 3.14 s

[nilesjohnson]

comment:9

Replying to @nilesjohnson:

This sage patch performs substantially better in 3/4 cases...

oops; should be "This patch performs substantially better in 4/5 cases..."

[simon-king-jena]

Changed keywords from none to multivariate power series

[simon-king-jena]

comment:11

First of all: It would certainly be good to have multivariate power series in Sage. Also the performance seems to be pretty good. Thank you!

I have some criticism, though.

1. __contains__ and coercion

I was reading part of the patch, and I see that your multivariate power series rings have a custom __contains__ method. I think it should be possible to simply inherit it.

Moreover, the custom method seems to provide a non-standard behaviour: Usual in Sage is to have a in R if and only if a==R(a) is True (of course, if R(a) raises an error then a is not in R).

Another formulation of the rule is: If a.parent() is P and R has a coercion map from P then a is in R: Since there is a coercion map, R(a) is supposed to work, and a==R(a) is equivalent to R(a)==R(a), because this is how comparison with the coercion model works.

But you have in your doctests

sage: a in R
False
sage: R(a) in R
True

in a situation were apparently a coercion from a.parent() to R exists.

2. Use of double-underscore methods

In your __cmp__ method, you do
self._bg_value.__cmp__(other._bg_value)
Generally, one should avoid to call double-underscore methods directly. So, better do
cmp(self._bg_value, other._bg_value)

The same applies to doctests. So, instead of

sage: R.<x,y> = MPowerSeriesRing(GF(17)) 
sage: R   
Multivariate Power Series Ring in x, y over Finite Field of size 17 
sage: R.__repr__() 
'Multivariate Power Series Ring in x, y over Finite Field of size 17'

you should do

sage: R     # indirect doctest
Multivariate Power Series Ring in x, y over Finite Field of size 17

and instead of

sage: M = MPowerSeriesRing(ZZ,5,'t');
sage: N = M.remove_var(M.gens()[2])
sage: M.__contains__(N.random_element(10))
False
sage: M.__contains__(M.random_element(10))
True

you should have

sage: M.random_element(10) in M
True

I think it should even be

sage: N.random_element(10) in M
True   # not False!

because, if I understand correctly, there is a coercion from N to M -- see point 1.

I am now running make ptestlong, and will then have a closer look at the code - so, no review yet. But I think you should address the points above.

[simon-king-jena]

comment:12

With make ptestlong, I obtained:

The following tests failed:

        sage -t  -long devel/sage/sage/rings/multi_power_series_ring.py # 15 doctests failed
        sage -t  -long devel/sage/sage/rings/multi_power_series_ring_element.py # 28 doctests failed

These are quite many failures. So, this needs work.

[simon-king-jena]

comment:13

PS:

It would be nice if the double square bracket notation would work in the multivariate case:

sage: QQ[['x']]
Power Series Ring in x over Rational Field
sage: QQ[['x','y']]
---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)

/home/king/SAGE/<ipython console> in <module>()

/home/king/SAGE/sage-4.4.2/local/lib/python2.6/site-packages/sage/rings/ring.so in sage.rings.ring.Ring.__getitem__ (sage/rings/ring.c:2612)()

NotImplementedError: Power series rings only implemented in 1 variable

[simon-king-jena]

Work Issues: Fix doctests; fix contains; add syntactical sugar

[nilesjohnson]

comment:14

Thanks for taking a look at this; I'll address the other work issues soon, but for now I wanted to check the doctests. There aren't any tests flagged as #long, so I didn't think sage -t -long should be different from sage -t. And in any case, I get

sage -t -long "devel/sage/sage/rings/multi_power_series_ring.py"
	 [2.7 s]
 
----------------------------------------------------------------------
All tests passed!
Total time for all tests: 2.7 seconds

sage -t -long "devel/sage/sage/rings/multi_power_series_ring_element.py"
	 [3.3 s]
 
----------------------------------------------------------------------
All tests passed!
Total time for all tests: 3.3 seconds

Do you get all tests passed without -long? If not, maybe you forgot to add the patches at #9457 and #9443 (mentioned in my first trac comment)? If you still have doctest failures, could you send me a little more about what's failing?

[simon-king-jena]

comment:15

Replying to @nilesjohnson:

Do you get all tests passed without -long? If not, maybe you forgot to add the patches at #9457 and #9443 (mentioned in my first trac comment)?

Outch! Sorry, I missed that. I think I'll not be able to resume work before Monday, but then I'll try again with #9457 and #9433.

Cheers, Simon

[nilesjohnson]

Attachment: trac_1956_combined_5.patch.gz

updated for new repr of completion functor

[nilesjohnson]

comment:128

A new patch, updated according to the repr string of the completion functor. No other changes.

[nilesjohnson]

comment:129

Patchbot: apply trac_1956_combined_5.patch

[jdemeyer]

comment:130

Please rebase properly to sage-4.7.rc2:

sage -t  -force_lib devel/sage/sage/rings/multi_power_series_ring_element.py
**********************************************************************
File "/mnt/usb1/scratch/jdemeyer/merger/sage-4.7.1.alpha1/devel/sage-main/sage/rings/multi_power_series_ring_element.py", line 1299:
    sage: f.sqrt()
Expected:
    Traceback (most recent call last):
    ...
    AttributeError: 'sage.rings.polynomial.multi_polynomial_libsingular.MPolynomial_libsingular' object has no attribute 'sqrt'
Got:
    Traceback (most recent call last):
      File "/mnt/usb1/scratch/jdemeyer/merger/sage-4.7.1.alpha1/local/bin/ncadoctest.py", line 1231, in run_one_test
        self.run_one_example(test, example, filename, compileflags)
      File "/mnt/usb1/scratch/jdemeyer/merger/sage-4.7.1.alpha1/local/bin/sagedoctest.py", line 38, in run_one_example
        OrigDocTestRunner.run_one_example(self, test, example, filename, compileflags)
      File "/mnt/usb1/scratch/jdemeyer/merger/sage-4.7.1.alpha1/local/bin/ncadoctest.py", line 1172, in run_one_example
        compileflags, 1) in test.globs
      File "<doctest __main__.example_39[4]>", line 1, in <module>
        f.sqrt()###line 1299:
    sage: f.sqrt()
      File "/mnt/usb1/scratch/jdemeyer/merger/sage-4.7.1.alpha1/local/lib/python/site-packages/sage/rings/multi_power_series_ring_element.py", line 1304, in sqrt
        return self.parent(self._bg_value.sqrt())
      File "power_series_ring_element.pyx", line 1205, in sage.rings.power_series_ring_element.PowerSeries.sqrt (sage/rings/power_series_ring_element.c:8919)
        s = u[0].sqrt(extend=False)
      File "element.pyx", line 2003, in sage.structure.element.CommutativeRingElement.sqrt (sage/structure/element.c:14899)
      File "element.pyx", line 1906, in sage.structure.element.CommutativeRingElement.is_square (sage/structure/element.c:14740)
    NotImplementedError: is_square() not implemented for elements of Multivariate Polynomial Ring in a, b over Integer Ring
**********************************************************************

[nilesjohnson]

comment:131

Replying to @jdemeyer:

Please rebase properly to sage-4.7.rc2:

sage -t  -force_lib devel/sage/sage/rings/multi_power_series_ring_element.py
**********************************************************************

I've addressed this problem, but I see other issues with sage/interfaces/maxima.py that I don't understand since this patch does not touch maxima . . . are those spurious failures?

[nilesjohnson]

Attachment: trac_1956_combined_6.patch.gz

[nilesjohnson]

comment:132

Ah; problems with interfaces/maxima.py are independent of this patch, since they occur in my newly-installed 4.7.rc2. Patch 6 here fixes the error raised by sqrt to be the same as that raised by square_root, thus passing the doctest as it did before. All tests in rings/multi_power_series* now pass.

Patchbot: apply trac_1956_combined_6.patch

[nilesjohnson]

comment:133

Can we get this switched back to positive review? Sage 4.7.rc3 fixes the problems with interfaces/maxima.py; this patch applies cleanly and passes all tests under that release.

[jdemeyer]

comment:134

Replying to @nilesjohnson:

Can we get this switched back to positive review?

bump

Any reviewers?

[malb]

comment:135

okay

[jdemeyer]

Merged: sage-4.7.1.alpha2