Opened 11 years ago
Closed 8 years ago
#1956 closed enhancement (fixed)
implement multivariate truncated power series arithmetic
Reported by: | was | Owned by: | pernici |
---|---|---|---|
Priority: | major | Milestone: | sage-4.7.1 |
Component: | commutative algebra | Keywords: | multivariate power series |
Cc: | mhansen, SimonKing, malb, hlaw | Merged in: | sage-4.7.1.alpha2 |
Authors: | Niles Johnson | Reviewers: | Martin Albrecht, Simon King |
Report Upstream: | N/A | Work issues: | |
Branch: | Commit: | ||
Dependencies: | Stopgaps: |
Description (last modified by )
Multivariate truncated power series arithmetic has been requested a *lot*.
Apply:
This is a single patch combining all of the partial patches now attached to this ticket.
Attachments (22)
Change History (158)
comment:1 Changed 11 years ago by
- Milestone changed from sage-2.10.1 to sage-feature
- Type changed from defect to enhancement
comment:2 Changed 11 years ago by
- Cc mhansen added
comment:3 Changed 9 years ago by
- Report Upstream set to N/A
- Status changed from new to needs_review
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.
comment:4 Changed 9 years ago by
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).
comment:5 follow-up: ↓ 6 Changed 9 years ago by
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?
comment:6 in reply to: ↑ 5 ; follow-up: ↓ 7 Changed 9 years ago by
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.
comment:7 in reply to: ↑ 6 ; follow-up: ↓ 8 Changed 9 years ago by
Replying to niles:
- 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?
comment:8 in reply to: ↑ 7 ; follow-up: ↓ 9 Changed 9 years ago by
Replying to malb:
Replying to niles:
- 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
comment:9 in reply to: ↑ 8 Changed 9 years ago by
Replying to niles:
This sage patch performs substantially better in 3/4 cases...
oops; should be "This patch performs substantially better in 4/5 cases..."
comment:10 Changed 9 years ago by
- Cc SimonKing added
comment:11 follow-up: ↓ 16 Changed 9 years ago by
- Keywords multivariate power series added
- Milestone changed from sage-feature to sage-5.0
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.
__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.
- 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.
comment:12 follow-up: ↓ 14 Changed 9 years ago by
- Status changed from needs_review to needs_work
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.
comment:13 Changed 9 years ago by
- Work issues set to Fix doctests; fix __contains__; add syntactical sugar
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
comment:14 in reply to: ↑ 12 ; follow-up: ↓ 15 Changed 9 years ago by
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?
comment:15 in reply to: ↑ 14 Changed 9 years ago by
Replying to niles:
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
comment:16 in reply to: ↑ 11 Changed 9 years ago by
- Work issues changed from Fix doctests; fix __contains__; add syntactical sugar to Fix doctests
Thanks for the tips!
Replying to SimonKing:
__contains__
and coercion
__contains__
is now inherited, and non-standard behaviour fixed (I wasn't aware of the standard).
- Use of double-underscore methods
I have fixed as many of these as I could find (probably all of them, but let me know if there are more issues).
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.
Indeed; I have found that my installation of sage fails some long doctests even before the patches are applied, so I will have to do a fresh installation before I can test this myself properly.
Also, I had some time to kill so I've added support for the double bracket method, and this is included in the documentation.
sage: ZZ[['s,t,u']] Multivariate Power Series Ring in s, t, u over Integer Ring
Finally, with a little more free time I added support for exponents
and a basic version of the verschiebung, V
. These were among the pile of not-yet-implemented methods inherited from univariate power series.
comment:17 Changed 9 years ago by
- Status changed from needs_work to needs_review
- Work issues Fix doctests deleted
Since the patch at #9457 seems to be running into trouble, I decided to see if I could make this patch independent from it. As it turns out, this has happened already somehow!
make ptestlong
just finished, passing all tests, with this patch and the one at #9443. To verify this, I recommend first doctesting just the files sage/rings/multi_power_series_ring.py
, sage/rings/multi_power_series_ring_element.py
, and sage/schemes/elliptic_curves/sha_tate.py
(with -long
); if these pass all tests, then make ptestlong
should too. I'm changing the status to "needs review" since all of the work issues have been addressed, but of course it can be switched back if there are further issues.
comment:18 Changed 9 years ago by
- Milestone changed from sage-5.0 to sage-4.5.3
Changing the milestone to 4.5.3 in an effort to generate attention. Since dependence on #9457 has been removed, there's no reason to wait until 5.0, is there?
comment:19 follow-up: ↓ 20 Changed 9 years ago by
I ran doctests on sage.math and those fail because (I assume) the random seed somehow differs. Niles, does it work for you on sage.math with 4.5.3?
comment:20 in reply to: ↑ 19 Changed 9 years ago by
Replying to malb:
I ran doctests on sage.math and those fail because (I assume) the random seed somehow differs. Niles, does it work for you on sage.math with 4.5.3?
Ah; I don't have an account on sage.math. Most of the uses of random_element could be avoided, so I'll submit a new patch without them.
Changed 9 years ago by
apply only this patch; removes needless uses of random_element, adds some doctests, cleans up some documentation
comment:21 Changed 9 years ago by
- Description modified (diff)
Attached a new, combined patch which removes some unnecessary uses of random_element
, cleans up some of the documentation, and adds doctests for NotImplemented
functions inherited from univariate power series (c.f. comment 4).
I now have (agian) the patch passing its own tests and documentation building without warning; make ptestlong
is currently in progress on sage.math.
comment:22 follow-up: ↓ 23 Changed 9 years ago by
- Status changed from needs_review to needs_work
I get two long doctest failures
File "/mnt/usb1/scratch/malb/sage-4.5.3/devel/sage/sage/rings/multi_power_series_ring.py", line 496: sage: W = MPowerSeriesRing(InfinitePolynomialRing(ZZ,'a'),2,'x,y') Exception raised: Traceback (most recent call last): ... TypeError: is_integral_domain() takes exactly 1 argument (2 given)
sage: W.is_noetherian() Exception raised: Traceback (most recent call last): ... NameError: name 'W' is not defined
I'll start reading the code now.
comment:23 in reply to: ↑ 22 Changed 9 years ago by
- Status changed from needs_work to needs_review
Replying to malb:
I get two long doctest failures
TypeError?: is_integral_domain() takes exactly 1 argument (2 given)
It is originally stated on this ticket that "one also needs the patches at #9457 and #9443". Did you take this into account?
The purpose of #9443 is to make is_integral_domain
accept optional arguments, and this should solve the problem. Therefore, I switch back to "needs review".
comment:24 Changed 9 years ago by
The code looks good modulo some tiny issues:
- I'll attach a referee patch in a sec which improves the Sphinx typesetting slightly
- I'm not happy about the docstring "This function works, but doesn't get called when it should." Either the function should get removed or the problem fixed.
Modulo those I give this patch a positive review.
PS:http://sage.math.washington.edu/home/malb/scratch_sage/sage-4.5.3/devel/sage/doc/output/html/en/reference/sage/rings/multi_power_series_ring_element.html to see how it looks like in the reference manual.
comment:25 Changed 9 years ago by
- Status changed from needs_review to needs_work
Sorry for not checking the original description carefully enough.
comment:26 Changed 9 years ago by
- Description modified (diff)
sorry; I should have mentioned the dependency in the description (fixed now).
While trying to write a detailed email explaining how I can't understand the coercion model, I realized that the problem with _l_action_
not being called is probably that the Completion
functor didn't know how to complete multivariate polynomial rings to multivariate power series rings . . . so that's fixed now, and it does fix the _l_action_
problem :)
After I check for new doctest failures, I'll upload the new patch.
Changed 9 years ago by
apply only this patch; fixes coercion problems coming from incomplete functorial construction
comment:27 Changed 9 years ago by
- Description modified (diff)
- Status changed from needs_work to needs_review
Ok, the latest patch trac_1956_multi_power_series_new_2.patch passes all tests in the sage/rings
directory, so I'm uploading it and starting ptestlong
.
At this point, I'd like to start thinking that any other limitations of this patch should be filed under a new ticket, but we'll see what comes up.
comment:28 Changed 9 years ago by
- Description modified (diff)
forgot to check the documentation! There's just one minor fix, and it reminded me to update the behavior of remove_var
to be consistent with the (new) behavior of remove_var
for multivariate polynomial rings. Here's the new patch: trac_1956_multi_power_series_new_3.patch.
sorry for the rapid-fire posts!
comment:29 Changed 9 years ago by
- Status changed from needs_review to positive_review
All looks good to me now.
comment:30 Changed 9 years ago by
- Reviewers set to Martin Albrecht, Simon King
comment:31 follow-up: ↓ 32 Changed 9 years ago by
- Status changed from positive_review to needs_work
I think there are still some problems with this. For example:
- Pickling does not work for these objects -- one needs to define a proper __reduce__ method.
- I don't think we should continue the use of __ attributes as they just cause problems. For example, the code in here has to set attributes like _PowerSeriesRing_generic__power_series_class in subclasses.
- There are formating issues in some of the docstrings like '_latex_'
- We shouldn't have a 'MPowerSeriesRing' -- that functionality should just be in PowerSeriesRing?. MPolynomialRing was deprecated for the same reason a long time ago.
comment:32 in reply to: ↑ 31 ; follow-up: ↓ 34 Changed 9 years ago by
- Status changed from needs_work to needs_review
Replying to mhansen:
I think there are still some problems with this. For example:
- Pickling does not work for these objects -- one needs to define a proper __reduce__ method.
done now; patch forthcoming.
- I don't think we should continue the use of __ attributes as they just cause problems. For example, the code in here has to set attributes like _PowerSeriesRing_generic__power_series_class in subclasses.
The use of __
attributes here is modeled after their use (for better or worse) in PowerSeriesRing_generic
(see 4. below); it is also in keeping with the python philosophy, that __
attributes or methods may be removed or changed in future versions. It is entirely possible that someone will eventually want to significantly rework multi/univariate power series. Having said that, they can be removed without too much trouble if sage development policy is strongly against them (in this case, perhaps the developer guide should mention it clearly somewhere).
- There are formating issues in some of the docstrings like '_latex_'
_latex_
is fixed now, and I looked for other problems but didn't see any; my best known method for checking docstring formatting is using sage -docbuild
, which only shows warnings/errors for non-underscore methods . . . can you recommend a better way?
- We shouldn't have a 'MPowerSeriesRing' -- that functionality should just be in PowerSeriesRing?. MPolynomialRing was deprecated for the same reason a long time ago.
I agree, but I think that issue is beyond the scope of this ticket. Implementing multivariate power series arithmetic is already long and complex enough, without having to also worry about integrating them with the univariate power series code (which is heavily specialized for the single-variable case). I have spent time thinking about how to do this, and I think it will be a subtle problem to solve; there is a sizable body of code that depends on PowerSeriesRing
as it stands, and one will have to be careful about extending its functionality without breaking any of the current uses. So I suggest that we merge this ticket and then open a new ticket for Unify construction of uni- and multi-variate power series rings, as has been done for polynomial rings.
It would probably be good if the other people involved in this ticket could comment on this last point.
comment:33 Changed 9 years ago by
- Description modified (diff)
comment:34 in reply to: ↑ 32 ; follow-up: ↓ 35 Changed 9 years ago by
Replying to niles:
I agree, but I think that issue is beyond the scope of this ticket. Implementing multivariate power series arithmetic is already long and complex enough, without having to also worry about integrating them with the univariate power series code (which is heavily specialized for the single-variable case). I have spent time thinking about how to do this, and I think it will be a subtle problem to solve; there is a sizable body of code that depends on
PowerSeriesRing
as it stands, and one will have to be careful about extending its functionality without breaking any of the current uses. So I suggest that we merge this ticket and then open a new ticket for Unify construction of uni- and multi-variate power series rings, as has been done for polynomial rings. It would probably be good if the other people involved in this ticket could comment on this last point
I think you may have misunderstood what I meant. In the patch, you've made a function "MPowerSeriesRing" which is basically a copy of the function "PowerSeriesRing?" defined in power_series_ring.py. It is just responsible for caching and returning the correct type. There is an assert statement in PowerSeriesRing? that makes sure that you're only trying to make a univariate power series ring. This function should just be changed to return the appropriate thing in the multivariate case. This is maybe like 10 lines of code or so (and changing doctests).
comment:35 in reply to: ↑ 34 Changed 9 years ago by
Replying to mhansen:
I think you may have misunderstood what I meant. In the patch, you've made a function "MPowerSeriesRing" which is basically a copy of the function "PowerSeriesRing?" defined in power_series_ring.py. It is just responsible for caching and returning the correct type. There is an assert statement in PowerSeriesRing? that makes sure that you're only trying to make a univariate power series ring. This function should just be changed to return the appropriate thing in the multivariate case. This is maybe like 10 lines of code or so (and changing doctests).
Well, I'm slightly embarrassed to say it, but I think did understand what you meant. I have tried a few times now, but each time I look at the code for PowerSeriesRing
, I find it confusing and a bit overwhelming. I thought perhaps looking at the code for PolynomialRing
would offer some guidance, but that only made things worse. It seems that each of these go to great lengths to allow (positional) arguments to be given in a variety of orders, and somehow I find this intimidating. (I have a hard time telling how various input cases are handled, for example.)
I'll spend a couple more days on it, but I won't be heartbroken if someone reassigns it to a new patch, or decides that it's so easy they can do it themselves in under an hour ;)
comment:36 follow-ups: ↓ 37 ↓ 38 Changed 9 years ago by
I'll update the patch to take care of this.
comment:37 in reply to: ↑ 36 Changed 9 years ago by
comment:38 in reply to: ↑ 36 Changed 9 years ago by
- Status changed from needs_review to needs_work
- Work issues set to remove `MPowerSeriesRing`, implement through `PowerSeriesRing`
For what it's worth, here is the specific issue I'm stuck on: PowerSeries
allows the argument default_prec
to be specified using either the first or second position:
sage: T = PowerSeriesRing(QQ,'t',3) sage: T.default_prec() 3 sage: T = PowerSeriesRing(QQ,3,names='t') sage: T.default_prec() 3 sage: T.variable_name() 't'
So what's the best way to preserve this functionality, but also allow creation of multivariate power series using the same (or similar) syntax as PolynomialRing
uses? e.g.
sage: R = PolynomialRing(QQ,3,'r'); R Multivariate Polynomial Ring in r0, r1, r2 over Rational Field sage: R = PolynomialRing(QQ,'r',3); R Multivariate Polynomial Ring in r0, r1, r2 over Rational Field
Note, on the other hand, that the following give errors:
sage: T = PowerSeriesRing(QQ,3,'t') Traceback (most recent call last): ... ValueError: first letter of variable name must be a letter sage: T = PowerSeriesRing(QQ,3,name='t') Traceback (most recent call last): ... TypeError: PowerSeriesRing() got multiple values for keyword argument 'name'
To me, it seems that the current behavior (assuming integer arguments are meant to be default_prec
) should be deprecated, and the allowed syntax of PowerSeriesRing
should mirror that of PolynomialRing
.
Oh, I just realized what I should try: Make PowerSeriesRing(QQ,3,'t')
and PowerSeriesRing(QQ,'t,u,v')
create multivariate power series rings, and raise deprecation warnings in the other cases (as long as the rest of sage doesn't depend too heavily on this behavior); then MPowerSeriesRing
never has to exist, PowerSeriesRing
will sort of work as expected, and eventually it can be brought properly in line with PolynomialRing
. How does that sound to other people here?
comment:39 Changed 9 years ago by
- Description modified (diff)
- Status changed from needs_work to needs_review
Ok, I've attached a new patch which carries out the plan above. After running doctests, it seems that there is a non-trivial body of code (elliptic curves, and maybe p-adics) that makes use of the syntax
T = PowerSeriesRing(QQ,'t',3)
or
T.<t> = PowerSeriesRing(QQ,3)
to construct univariate power series rings. I will mention this and a few other things in a new ticket for merging univariate and multivariate power series rings.
This new patch passes all doctests (with -long), and documentation builds cleanly, looks good.
comment:40 Changed 9 years ago by
- Work issues remove `MPowerSeriesRing`, implement through `PowerSeriesRing` deleted
comment:41 Changed 9 years ago by
- Description modified (diff)
Note: these patches apply cleanly to sage 4.6.alpha1
comment:42 Changed 9 years ago by
- Description modified (diff)
comment:43 Changed 9 years ago by
- Cc malb added
Hello,
Although we've probably missed the 4.6 release by now, maybe someone could take a look at this in time for the 4.7 release?
best, Niles
comment:44 follow-up: ↓ 45 Changed 9 years ago by
Hi, I tried to apply the patch against 4.6.alpha3 and it failed (bitrot). I can take a look once that's fixed but I think it would be best if Mike would take a look as well.
comment:45 in reply to: ↑ 44 Changed 9 years ago by
- Description modified (diff)
Replying to malb:
Hi, I tried to apply the patch against 4.6.alpha3 and it failed (bitrot). I can take a look once that's fixed but I think it would be best if Mike would take a look as well.
I finally got 4.6.alpha3 built, but I didn't have any errors applying the patches . . . can you double check that you first applied trac_1956_multi_power_series_new_4.patch and then trac_1956_uni_multi_ps_2.patch? (I did get errors when applying only the second of these) If there are still errors, could you send me a snippet of the failure message?
Thanks for staying with this!
comment:46 follow-up: ↓ 47 Changed 9 years ago by
You are right, I cannot reproduce my problem. As instructed above, I can indeed apply the patch cleanly.
Some more comments (I'm really not the best person to review this code):
Here are some performance figures:
sage: R.<t,u,v> = PowerSeriesRing(QQ); R Multivariate Power Series Ring in t, u, v over Rational Field sage: f = R.random_element(prec=20,bound=20) sage: fp = f.polynomial() sage: %timeit f^2 # power series 125 loops, best of 3: 2.89 ms per loop sage: %timeit fp^2 # polynomials 625 loops, best of 3: 82.4 µs per loop sage: R.<t,u,v> = PowerSeriesRing(GF(127)); R Multivariate Power Series Ring in t, u, v over Finite Field of size 127 sage: f = R.random_element(prec=20,bound=20) sage: fp = f.polynomial() sage: %timeit f^2 # power series 625 loops, best of 3: 1.06 ms per loop sage: %timeit fp^2 # polynomials 625 loops, best of 3: 15.3 µs per loop
I also noticed that this is rather slow:
sage: sage: R.<t,u,v> = PowerSeriesRing(GF(127)); R Multivariate Power Series Ring in t, u, v over Finite Field of size 127 sage: f = R.random_element(prec=20,bound=20) sage: p = f.polynomial() sage: fp = f.polynomial() sage: %timeit R(fp) 25 loops, best of 3: 11.1 ms per loop
Conversion to Magma didn't seem to work at all, or am I doing it wrong?
sage: magma(f) In file "/home/malb/.sage//temp/road/7093//interface//tmp7130", line 1, column 16: >> _sage_[4]:=-63*u^4 + 11*t*v^8 - 17*t^3*u^4*v^3 - 61*t^4*u^5*v^2 + 16*t^5*u^ ^ User error: Identifier 'u' has not been declared or assigned
comment:47 in reply to: ↑ 46 ; follow-ups: ↓ 48 ↓ 49 Changed 9 years ago by
Replying to malb:
Here are some performance figures: <SNIP>
Thanks; I think these are outside the scope of this ticket (although important!) Here are the numbers I get:
sage: R.<t,u,v> = PowerSeriesRing(QQ); R Multivariate Power Series Ring in t, u, v over Rational Field sage: f = R.random_element(prec=20,bound=20) sage: fp = f.polynomial() sage: %timeit f^2 # power series 125 loops, best of 3: 2.67 ms per loop sage: %timeit fp^2 # polynomials 625 loops, best of 3: 52.3 µs per loop
The arithmetic is done in a univariate power series ring over a multivariate power series ring (called the "background ring" in the code):
sage: R._bg_ps_ring() Power Series Ring in T over Multivariate Polynomial Ring in t, u, v over Rational Field sage: fb = f._bg_value sage: %timeit fb^2 125 loops, best of 3: 2.61 ms per loop
So the speed is a limitation of univariate power series rings over multivariate polynomial rings . . .
sage: R.<t,u,v> = PowerSeriesRing(GF(127)); R Multivariate Power Series Ring in t, u, v over Finite Field of size 127 sage: f = R.random_element(prec=20,bound=20) sage: fp = f.polynomial() sage: %timeit f^2 # power series 625 loops, best of 3: 1.48 ms per loop sage: %timeit fp^2 # polynomials 625 loops, best of 3: 20.8 µs per loop sage: fb = f._bg_value # univariate power series over multivariate polynomials sage: %timeit fb^2 625 loops, best of 3: 1.42 ms per loop
Conversion to Magma didn't seem to work at all, or am I doing it wrong?
I don't know anything about converting to magma, but the following are also broken:
sage: magma(fp) TypeError: Error evaluating Magma code. IN:_sage_[5]:=SageCreateWithNames(PolynomialRing(_sage_ref2,3,negdeglex),["t","u","v"]); OUT: >> _sage_[5]:=SageCreateWithNames(PolynomialRing(_sage_ref2,3,negdeglex),["t", ^ User error: Identifier 'negdeglex' has not been declared or assigned sage:magma(fb) >> _sage_[6]:=32*t^4*T^4 + 30*t^2*u*v^5*T^8 - 21*t^3*u^3*v^3*T^9 + 56*t*u^7*v^ ^ User error: Identifier 't' has not been declared or assigned
comment:48 in reply to: ↑ 47 Changed 9 years ago by
Replying to niles:
I don't know anything about converting to magma, but the following are also broken:
sage: magma(fp) TypeError: Error evaluating Magma code. IN:_sage_[5]:=SageCreateWithNames(PolynomialRing(_sage_ref2,3,negdeglex),["t","u","v"]); OUT: >> _sage_[5]:=SageCreateWithNames(PolynomialRing(_sage_ref2,3,negdeglex),["t", ^ User error: Identifier 'negdeglex' has not been declared or assigned
Magma doesn't support local orderings for multivariate polynomial rings.
comment:49 in reply to: ↑ 47 Changed 9 years ago by
Replying to niles:
Mike, would you be able to take a look at this some time?
Here are some further comparisons with polynomial arithmetic, computing fk
for various values of
k
. Obviously the power series computation has a greater and greater advantage over the polynomial computation as
k
grows, since the polynomial arithmetic computes many more (irrelevant) terms. So the time tests below aren't exactly fair -- one could compute the necessary terms by computing small polynomial powers, truncating, and repeating -- but part of the point of the power series code is to automate this for the user. One side effect is that arithmetic for small powers is slower, but for large powers the power series arithmetic is fairly good: below, I have a comparison in a simple special case:
k = 2n
and non-zero constant coefficient.
In any case, the point of this ticket is to implement the basic arithmetic. The code should be fast enough for basic use, and speed improvement can be part of a separate ticket.
Timings for
fk
v.s.
fpk
sage: BaseRingList = [ZZ,QQ,GF(11),GF(65537)] sage: for B in BaseRingList: R = PowerSeriesRing(B,4,'t'); R f = R.random_element(prec=10,bound=20) fp = f.polynomial() for k in [2,5,8,11]: print 'f**'+str(k)+': (power series)' timeit('f**k',number=5) print 'fp**'+str(k)+': (polynomial)' timeit('fp**k',number=5) print '--' ....: Multivariate Power Series Ring in t0, t1, t2, t3 over Integer Ring f**2: (power series) 5 loops, best of 3: 879 µs per loop fp**2: (polynomial) 5 loops, best of 3: 33.2 µs per loop -- f**5: (power series) 5 loops, best of 3: 4.12 ms per loop fp**5: (polynomial) 5 loops, best of 3: 2.52 ms per loop -- f**8: (power series) 5 loops, best of 3: 6.47 ms per loop fp**8: (polynomial) 5 loops, best of 3: 51.5 ms per loop -- f**11: (power series) 5 loops, best of 3: 10.8 ms per loop fp**11: (polynomial) 5 loops, best of 3: 652 ms per loop -- Multivariate Power Series Ring in t0, t1, t2, t3 over Rational Field f**2: (power series) 5 loops, best of 3: 1 ms per loop fp**2: (polynomial) 5 loops, best of 3: 44.8 µs per loop -- f**5: (power series) 5 loops, best of 3: 8.68 ms per loop fp**5: (polynomial) 5 loops, best of 3: 4.12 ms per loop -- f**8: (power series) 5 loops, best of 3: 26.4 ms per loop fp**8: (polynomial) 5 loops, best of 3: 119 ms per loop -- f**11: (power series) 5 loops, best of 3: 58.4 ms per loop fp**11: (polynomial) 5 loops, best of 3: 2.28 s per loop -- Multivariate Power Series Ring in t0, t1, t2, t3 over Finite Field of size 11 f**2: (power series) 5 loops, best of 3: 480 µs per loop fp**2: (polynomial) 5 loops, best of 3: 10.4 µs per loop -- f**5: (power series) 5 loops, best of 3: 2.49 ms per loop fp**5: (polynomial) 5 loops, best of 3: 547 µs per loop -- f**8: (power series) 5 loops, best of 3: 4.23 ms per loop fp**8: (polynomial) 5 loops, best of 3: 11 ms per loop -- f**11: (power series) 5 loops, best of 3: 8.6 ms per loop fp**11: (polynomial) 5 loops, best of 3: 190 ms per loop -- Multivariate Power Series Ring in t0, t1, t2, t3 over Finite Field of size 65537 f**2: (power series) 5 loops, best of 3: 614 µs per loop fp**2: (polynomial) 5 loops, best of 3: 18.8 µs per loop -- f**5: (power series) 5 loops, best of 3: 4.37 ms per loop fp**5: (polynomial) 5 loops, best of 3: 3.13 ms per loop -- f**8: (power series) 5 loops, best of 3: 16.5 ms per loop fp**8: (polynomial) 5 loops, best of 3: 231 ms per loop -- f**11: (power series) 5 loops, best of 3: 29.7 ms per loop fp**11: (polynomial) 5 loops, best of 3: 4.85 s per loop --
Timings for
f(2n)
v.s. iterated square-and-truncate for
fp(2n)
Here is a simple special case: for a power series f
with non-zero constant coefficient, the total-degree precision of f**k
is equal to the total-degree precision of f
for all k
. Below is a function which uses polynomial arithmetic to compute coefficients of fp**k
up to given total degree for k = 2^n
, where fp = f.polynomial()
.
sage: def tot_deg_truncate(g,d): return sum(g.monomial_coefficient(m)*m for m in [x for x in g.monomials() if x.degree() < d]) ....: sage: def pow_recursive(g,k,d): if k == 2: return tot_deg_truncate(g**2,d) else: return pow_recursive(tot_deg_truncate(g**2,d),k/2,d) ....: sage: R = PowerSeriesRing(ZZ,4,'t') sage: f = 1 + R.random_element(prec=10,bound=20) sage: f.constant_coefficient() 1 sage: fp = f.polynomial()
Check that pow_recursive
gets the right answer:
sage: f**16 - pow_recursive(fp,16,10) 0 + O(t0, t1, t2, t3)^10 sage: f**64 - pow_recursive(fp,64,10) 0 + O(t0, t1, t2, t3)^10
Compare timings; although the polynomial arithmetic is faster, they are on roughly the same order of magnitude.
sage: %timeit f**16 125 loops, best of 3: 4.8 ms per loop sage: %timeit pow_recursive(fp,16,10) 625 loops, best of 3: 1.18 ms per loop sage: %timeit f**64 125 loops, best of 3: 7.21 ms per loop sage: %timeit pow_recursive(fp,64,10) 125 loops, best of 3: 1.76 ms per loop
Increasing the precision exaggerates the difference:
sage: f = 1 + R.random_element(prec=30,bound=20) sage: f.constant_coefficient() 1 sage: fp = f.polynomial() sage: f**16 - pow_recursive(fp,16,30) 0 + O(t0, t1, t2, t3)^30 sage: f**64 - pow_recursive(fp,64,30) 0 + O(t0, t1, t2, t3)^30 sage: %timeit f**16 25 loops, best of 3: 15.9 ms per loop sage: %timeit pow_recursive(fp,16,30) 125 loops, best of 3: 1.81 ms per loop sage: %timeit f**64 25 loops, best of 3: 23.9 ms per loop sage: %timeit pow_recursive(fp,64,30) 125 loops, best of 3: 2.65 ms per loop
comment:50 Changed 9 years ago by
I just ran long doctests on 4.6.1.alpha2 and they still pass. I'll ask people for feedback on [sage-devel], if there is none within a few days, I'll just give this patch a positive review. This patch shouldn't have bit rotted so long!
comment:51 follow-up: ↓ 55 Changed 9 years ago by
- Status changed from needs_review to needs_work
Looking with my eyes... I found one bug not mentioned or fixed elsewhere:
733 except TypeError,AttributeError:
This should be
733 except (TypeError,AttributeError):
otherwise, AttributeError? gets *set* to the string that the TypeError? raises. Do check for other cases of this issue in the code.
This scares me a little:
raise TypeError("action of "+c.parent()+" on "+f.parent()+" not defined.")
the reason is that any exception such that constructing a string as part of it could take time, can potentially *MASSIVELY* slow down Sage. Just do:
raise TypeError
or
raise TypeError, "action not defined"
instead. Since otherwise, constructing that string might dominate arithmetic... (the coercion model will try to do the arithmetic, get the error, then try something else). Also, it is easy enough using "%debug" to see what all the relevant variables are when one gets an exception.
I see this also
1062 except:
I am against using naked excepts unless their is a very good reason to do so.
I also think that whenever possible all the inputs to a function should be documented. For example look at this one:
109 109 def completion(self, p, prec=20, extras=None): 110 """ 111 Return the completion of self with respect to the ideal generated 112 by the variable(s) ``p``. 113 114 INPUT: 115 116 - ``p`` -- variable or tuple of variables 117 118 EXAMPLES::
I think prec and extras should be documented. Especially important is "extras", since that is not self-explanatory at all. There is also a function def solve_linear_de(self, prec = infinity, b=None, f0=None):
with no docs on its input. Since it just does raise NotImplementedError
I think that's fine.
I don't like all this "type(obj)==type" stuff. I think isinstance(obj,type) is much better:
140 if p in self or type(p) == str and set(p).issubset(set([str(g) for g in self.gens()])): 141 p = tuple([p]) 142 elif type(p) == list or type(p) == tuple:
comment:52 follow-up: ↓ 53 Changed 9 years ago by
I should have added -- having read/skimmed through all the code, I think this is a *wonderful* contribution to Sage, and it's really excellent code. Thanks for finally pulling this off!!
comment:53 in reply to: ↑ 52 ; follow-up: ↓ 54 Changed 9 years ago by
Replying to was:
I should have added -- having read/skimmed through all the code, I think this is a *wonderful* contribution to Sage, and it's really excellent code. Thanks for finally pulling this off!!
I strongly agree -- and I'm sorry that I have not had time to contribute, not even by reviewing. I certainly expect to use this code.
Reading through the posts on this ticket, I think it is a really good example of good-natured cooperation and constructive criticism producing something excellent.
comment:54 in reply to: ↑ 53 Changed 9 years ago by
comment:55 in reply to: ↑ 51 Changed 8 years ago by
- Description modified (diff)
- Status changed from needs_work to needs_review
Thanks for the comments; these have been addressed now.
Replying to was:
733 except TypeError,AttributeError:
This should be
733 except (TypeError,AttributeError):
fixed.
This scares me a little:
raise TypeError("action of "+c.parent()+" on "+f.parent()+" not defined.")
the reason is that any exception such that constructing a string as part of it could take time, can potentially *MASSIVELY* slow down Sage. Just do:
raise TypeError
or
raise TypeError, "action not defined"
instead.
This has been fixed, and I've checked every occurrence of "raise ..." for similar problems (there are quite a few, since I've never heard of %debug, or pdb, until now -- perhaps this would be a useful addition to the developer's guide)
I see this also
1062 except:
There were two naked excepts; in both cases the code has been rewritten to avoid them.
I also think that whenever possible all the inputs to a function should be documented. For example look at this one ...
I've double-checked all functions now; the example you gave was one of the toughest because I was modifying an existing function (from the multivariate polynomial code). In any case, it's fixed now.
I don't like all this "type(obj)==type" stuff. I think isinstance(obj,type) is much better:
fixed.
comment:56 Changed 8 years ago by
Apply trac_1956_multi_power_series_new_4.patch, trac_1956_uni_multi_ps_2.patch, trac_1956_multi_ps_cleanup.patch
comment:57 Changed 8 years ago by
Hi,
I opened ticket 10480 on fast PowerSeries_poly multiplication. Using it multiplication in the multivariate power series implemented here should also be much faster, since it uses univariate power series multiplication internally.
Mario
comment:58 Changed 8 years ago by
Anyone interested in setting this to positive review?
Changed 8 years ago by
Changed 8 years ago by
Changed 8 years ago by
comment:59 follow-up: ↓ 60 Changed 8 years ago by
- Owner changed from malb to pernici
After applying trac_1956_multi_power_series_new_4.patch, trac_1956_uni_multi_ps_2.patch, trac_1956_multi_ps_cleanup.patch (main version), apply trac_10480_fast_PowerSeries_poly_multiplication2.patch from ticket #10480 and trac_1956_faster_MPowerSeries_mul.patch (version (1) in benchmarks below)
In the attached patch MPowerSeries._mul_ uses the do_mul_trunc_generic algorithm in polynomial_element.pyx.
The first benchmark has the 8 inversions posted here by Niles to compare with Magma, but with precisions multiplied by N; times are in seconds; processor:4-core: Intel Core i7 CPU 860 @ 2.80GHz on x86_64 GNU/Linux
test no. 1 2 3 4 5 6 7 8 N=1 main 2.72 0.02 0.18 10.5 0.05 1.55 2.33 2.61 (1) 0.05 0.01 0.05 0.24 0.03 0.04 0.07 0.09 N=2 main 290 1.7 0.75 1024 0.16 49.0 64.3 73.8 (1) 2.37 0.02 0.07 5.53 0.05 0.68 1.24 1.41 N=3 main 7743 67 4.4 24272 0.84 432 566 651 (1) 32 0.12 0.11 56 0.08 5.2 7.8 9.1
The next benchmark takes a random series and performs p = p*(p+1) N times; the same series is used in testing with the two versions; see attached mu.sage
N=1 2 3 4 n=2 prec=20 bound=20 main 0.002 0.005 0.005 0.005 (1) 6e-4 6e-4 6e-4 5e-4 prec=100 bound=100 main 0.084 0.81 2.6 3.3 (1) 0.003 0.005 0.006 0.006 prec=200 bound=200 main 0.52 153 3932 7337 (1) 0.03 0.31 1.1 1.2 n=4 prec=20 bound=20 main 0.003 0.003 0.003 0.003 (1) 6e-4 5e-4 5e-4 5e-4 prec=100 bound=100 main 0.087 0.19 0.19 0.19 (1) 0.002 0.002 0.002 0.002 prec=200 bound=200 main 0.59 1.97 1.91 1.66 (1) 0.008 0.008 0.008 0.008
Applying only trac_10480_fast_PowerSeries_poly_multiplication2.patch, not trac_1956_faster_MPowerSeries_mul.patch, MPowerSeries._mul_ uses univariate series multiplication, which uses the do_mul_trunc algorithm instead of the do_mul_trunc_generic algorithm; the latter is slower than the unpatched version in some test in un.sage in #10480, while the former is always faster. The latter algorithm is always faster than main in the tests above, and most of the times a few times faster than the former algorithm, so it seems a better choice for multivariate series.
With this patch _send_to_bg and _send_to_fg are faster; using the attached benchmark file bg.sage, n=number of variables, h=precision, N=100 random polynomials, bg=_send_to_bg, fg=_send_to_fg
n=2 h=20 n=2 h=100 n=4 h=20 n=4 h=100 bg fg bg fg bg fg bg fg main 0.89 0.02 42 0.06 1.18 0.02 58 0.06 (1) 0.019 0.002 0.09 0.02 0.018 0.002 0.09 0.01
comment:60 in reply to: ↑ 59 Changed 8 years ago by
Replying to pernici:
This is great -- on my machine adding your trac_1956_faster_MPowerSeries_mul.patch and trac_10480_fast_PowerSeries_poly_multiplication2.patch
from #10480 drops the total time for doctests in sage.rings.multi_power_series_ring_element
from 3.5 seconds to 2.8 seconds.
After reading the discussion at #10480, it looks like speeding up the multiplication is a somewhat subtle problem, so I propose opening a separate ticket for this issue. Then we can look into the relative advantages of do_mul_trunc_generic
and the Karatsuba algorithm, and patch MPowerSeries._mul_
to take advantage of the optimal one. I'll do this unless you think #10480 already includes this issue.
comment:61 Changed 8 years ago by
Oh, I think one of the earlier posts confused the buildbot (it scans for the word "apply"). Here's a clarification:
Apply trac_1956_multi_power_series_new_4.patch, trac_1956_uni_multi_ps_2.patch, trac_1956_multi_ps_cleanup.patch
comment:62 follow-up: ↓ 63 Changed 8 years ago by
Niles wrote:
After reading the discussion at #10480, it looks like speeding up the multiplication is a somewhat subtle problem, so I propose opening a separate ticket for this issue. Then we can look into the relative advantages of do_mul_trunc_generic and the Karatsuba algorithm, and patch MPowerSeries._mul_ to take advantage of the optimal one. I'll do this unless you think #10480 already includes this issue.
Nice idea to have a separate ticket for MPowerSeries._mul_ .
comment:63 in reply to: ↑ 62 Changed 8 years ago by
comment:64 follow-up: ↓ 65 Changed 8 years ago by
Multivariate series in one variable differ from univariate series; is it the intended behaviour?
sage: R = PowerSeriesRing(GF(127),'t');R Power Series Ring in t over Finite Field of size 127 sage: %timeit R.random_element(100) 625 loops, best of 3: 856 µs per loop sage: K = PowerSeriesRing(GF(127),1,'a'); K Multivariate Power Series Ring in a over Finite Field of size 127 sage: %timeit K.random_element(100) 5 loops, best of 3: 146 ms per loop
I noticed this since the patch trac_1956_faster_MPowerSeries_mul.patch does not work in this case, because _send_to_bg calls the degrees method
sage: %timeit K.random_element(100) AttributeError: 'sage.rings.polynomial.polynomial_zmod_flint.Polynomial_zmod_flint' object has no attribute 'degrees'
If multivariate series in one variable are meant to differ from univariate series I will modify _send_to_bg to deal with this case.
comment:65 in reply to: ↑ 64 Changed 8 years ago by
- Status changed from needs_review to needs_work
- Work issues set to multivariate series on 1 generator
Replying to pernici:
Multivariate series in one variable differ from univariate series; is it the intended behaviour?
No, they should be the same -- I'll fix this pretty soon.
comment:66 Changed 8 years ago by
Replying to pernici:
Multivariate series in one variable differ from univariate series; is it the intended behaviour?
I don't know whether it is intended, but I'd like to mention that there already is a difference between univariate polynomials and multivariate polynomials in one variable:
sage: R_uni = PolynomialRing(QQ,'x') sage: R_uni Univariate Polynomial Ring in x over Rational Field sage: R_multi = PolynomialRing(QQ,'x',1) sage: R_multi Multivariate Polynomial Ring in x over Rational Field sage: timeit('a = R_uni.random_element()') 625 loops, best of 3: 48.4 Âµs per loop sage: timeit('a = R_multi.random_element()') 625 loops, best of 3: 192 Âµs per loop sage: a = R_uni.random_element() sage: b = R_multi(a) sage: a.leading_coefficient() -27 sage: hasattr(b,'leading_coefficient') False sage: b.lc() -27 sage: hasattr(a,'lc') False
While it is clear that there is a difference in the timings for random_element
, I don't like that the names are different for methods that do essentially the same.
Things are different in the case of degrees
(which exists only for multivariate polynomials). Since the word is plural (it denotes the tuple of maximal exponents of each variable, not necessarily occuring in a single monomial), it doesn't really make sense in the univariate case. However, I do think that in that case (and similar cases) there should be a method of univariate polynomials emulating the corresponding method for multivariate polynomials with one variable.
So, I don't mind about the different timings; but I think methods should be more or less equivalent.
Changed 8 years ago by
apply on top of other patches; return univariate ring rather than multivariate ring in one variable
comment:67 follow-up: ↓ 68 Changed 8 years ago by
- Description modified (diff)
- Status changed from needs_work to needs_review
I think that the top-level constructors PolynomialRing
and PowerSeriesRing
should return univariate rings instead of multivariate rings in one variable. It is difficult to imagine a situation where someone would want access to the "multivariate in one variable" versions (e.g. they want to use some algorithm implemented for the multivariate case that is not available in the univariate case), and even more difficult to imagine that this is the preferred option of most users. Someone who does need this functionality can still import the multivariate constructor directly.
With trac_1956_one_variable_fix.patch we have:
sage: PowerSeriesRing(QQ,1,'a') Power Series Ring in a over Rational Field sage: PowerSeriesRing(QQ,1,'a') is PowerSeriesRing(QQ,'a') True
sage: from sage.rings.multi_power_series_ring import MPowerSeriesRing_generic sage: MPowerSeriesRing_generic(QQ,1,'a') Multivariate Power Series Ring in a over Rational Field
Unless someone here can think of a compelling reason to reverse this change, I think we should also open a new ticket for the polynomial case. In a related vein, we also have polynomial rings in no variables, and these are different from the base ring:
sage: PolynomialRing(ZZ,'a',0) Multivariate Polynomial Ring in no variables over Integer Ring sage: PolynomialRing(ZZ,'a',0) is ZZ False
For consistency, I've implemented the same behavior for power series rings:
sage: PowerSeriesRing(QQ,0,'a') Multivariate Power Series Ring in no variables over Rational Field
Perhaps these should both simply return the base ring -- I'm not so sure, but in any case I hope it can be handled for both polynomials and power series in a separate ticket.
comment:68 in reply to: ↑ 67 ; follow-up: ↓ 69 Changed 8 years ago by
- Status changed from needs_review to needs_work
- Work issues changed from multivariate series on 1 generator to multivariate series on 1 generator should remain different from a univariate series
Replying to niles:
I think that the top-level constructors
PolynomialRing
andPowerSeriesRing
should return univariate rings instead of multivariate rings in one variable.
It already does:
sage: R.<x> = QQ[] sage: R Univariate Polynomial Ring in x over Rational Field sage: R2.<x> = PolynomialRing(QQ) sage: R2 Univariate Polynomial Ring in x over Rational Field
Only when you explicitly request it, you'll get a multivariate ring with one variable:
sage: R3.<x> = PolynomialRing(QQ,1) sage: R3 Multivariate Polynomial Ring in x over Rational Field
It is difficult to imagine a situation where someone would want access to the "multivariate in one variable" versions
I find it very easy to imagine. I sometimes want it so, because multi- and univariate polynomials have different methods, and thus I want to make sure that my programs will always get a multivariate polynomial.
(e.g. they want to use some algorithm implemented for the multivariate case that is not available in the univariate case), and even more difficult to imagine that this is the preferred option of most users.
Agreed, and this is why the polynomial ring constructor returns a univariate ring, unless requested otherwise (in contrast to your claim).
Someone who does need this functionality can still import the multivariate constructor directly.
No, one shouldn't, it is deprecated:
sage: MPolynomialRing(QQ,1,'a') /mnt/local/king/SAGE/sage-4.6/local/bin/sage-ipython:1: DeprecationWarning: MPolynomialRing is deprecated, use PolynomialRing instead! #!/usr/bin/env python Multivariate Polynomial Ring in a over Rational Field
Instead, one should choose the arguments of the polynomial ring contructor appropriately (as I did above).
With trac_1956_one_variable_fix.patch we have:
sage: PowerSeriesRing(QQ,1,'a') Power Series Ring in a over Rational Field sage: PowerSeriesRing(QQ,1,'a') is PowerSeriesRing(QQ,'a') True
I am strongly -1 to that suggestion, because the user must have the possibility to choose the implementation by providing appropriate arguments (with the argument 1
versus without it).
sage: from sage.rings.multi_power_series_ring import MPowerSeriesRing_generic sage: MPowerSeriesRing_generic(QQ,1,'a') Multivariate Power Series Ring in a over Rational Field
Using a non-generic constructor should be deprecated, IMO.
Unless someone here can think of a compelling reason to reverse this change, I think we should also open a new ticket for the polynomial case.
I don't know if you find my argument compelling, but I am strongly against that suggestion.
In a related vein, we also have polynomial rings in no variables, and these are different from the base ring:
sage: PolynomialRing(ZZ,'a',0) Multivariate Polynomial Ring in no variables over Integer Ring sage: PolynomialRing(ZZ,'a',0) is ZZ FalseFor consistency, I've implemented the same behavior for power series rings:
sage: PowerSeriesRing(QQ,0,'a') Multivariate Power Series Ring in no variables over Rational Field
That's fine.
Perhaps these should both simply return the base ring --
No! The user has explicitly requested a polynomial ring. So, s/he must not get the base ring in return, since s/he may rely on particular methods that may only be implemented for polynomial rings, but not for general rings.
comment:69 in reply to: ↑ 68 ; follow-up: ↓ 70 Changed 8 years ago by
- Description modified (diff)
- Status changed from needs_work to needs_info
Hi Simon!
Replying to SimonKing:
Replying to niles:
It is difficult to imagine a situation where someone would want access to the "multivariate in one variable" versions
I find it very easy to imagine. I sometimes want it so, because multi- and univariate polynomials have different methods, and thus I want to make sure that my programs will always get a multivariate polynomial.
(e.g. they want to use some algorithm implemented for the multivariate case that is not available in the univariate case), and even more difficult to imagine that this is the preferred option of most users.
Agreed, and this is why the polynomial ring constructor returns a univariate ring, unless requested otherwise (in contrast to your claim).
Thanks for this -- actually I find your arguments here quite compelling. I was thinking before that arithmetic for univariate power series and polynomials are probably optimized for that case, and so would be preferable to the "multivariate in one variable" algorithms. But I neglected, as you have pointed out, that one might want to write code which, for simplicity, treats univariate and multivariate rings the same and thus depends on the methods written for multivariate rings. This would be especially likely if the code one were writing didn't depend on the optimal univariate algorithms.
Of course this is precisely not the case for pernici and others who are working on faster multiplication algorithms. Pernici, does this seem reasonable to you? How difficult will it be for you to treat the "multivariate in one variable" case?
comment:70 in reply to: ↑ 69 ; follow-up: ↓ 71 Changed 8 years ago by
Replying to niles:
Hi Simon!
Replying to SimonKing: Thanks for this -- actually I find your arguments here quite compelling. I was thinking before that arithmetic for univariate power series and polynomials are probably optimized for that case, and so would be preferable to the "multivariate in one variable" algorithms.
Sure univariate rings are usually preferable! I do think that most users could care less about the different methods of univariate and singly multivariate polynomials - they want speed!
Of course this is precisely not the case for pernici and others who are working on faster multiplication algorithms. Pernici, does this seem reasonable to you? How difficult will it be for you to treat the "multivariate in one variable" case?
I did not suggest that by default the polynomial ring (or power series ring) constructor should return a multivariate ring in one variable. I do think that
sage: R.<x> = QQ[[]]
should turn R
into a univariate ring.
What I was trying to explain was: The user must be able to override the default implementation by passing appropriate arguments to the ring constructor. Hence: PowerSeriesRing(QQ,'x')
should return the default (a univariate ring), but PowerSeriesRing(QQ,1,'x')
should return a multivariate ring.
In other words: I suggest that you just mimmick the current behaviour of th polynomial ring constructor.
Cheers, Simon
comment:71 in reply to: ↑ 70 Changed 8 years ago by
- Status changed from needs_info to needs_review
Replying to SimonKing:
In other words: I suggest that you just mimmick the current behaviour of th polynomial ring constructor.
Good--that's the behavior without the "one variable fix" patch :)
This morning I had an idea for helping this code get reviewed (note that the ticket turns 4 at the end of this month): Could we break the review process into several manageable pieces, and have different people take responsibility for different parts? For example, I'd propose something like the following:
- The underlying concept of the implementation (foreground ring & background ring, etc.)
- Documentation and tests
multi_power_series_ring.py
: the code accurately does what it claims to domulti_power_series_ring_element.py
: the code accurately does what it claims to do- Integration with the rest of sage: construction and use of
PowerSeriesRings
works correctly, and parallels behavior of polynomial rings - Completeness of the review: do these pieces give an essentially complete breakdown of the review?
comment:72 Changed 8 years ago by
I think that random_element does not give good random multivariate series; the lowest total degree is usually not zero
sage: R.<x,y> = QQ[[]] sage: def mval(p,N): ....: print '%.2f' % (sum([sum(R.random_element(p,p).exponents()[0]) for _ in range(N)])/float(N)) ....: sage: mval(20,100) 3.94 sage: mval(100,100) 8.34
while in the univariate case it is usually zero
sage: R.<x> = QQ[[]] sage: def mval(p,N): ....: print '%.2f' % (sum([R.random_element(p,p).exponents()[0] for _ in range(N)])/float(N)) ....: sage: mval(20,100) 0.02 sage: mval(100,100) 0.01
As a quick fix to get pseudo random series I changed the benchmark file mu.sage shifting by 1 the variables in the random polynomials
sage: R.<x,y> = QQ[[]] sage: for prec in range(10,50,10): ....: p = R.random_element(prec,prec) ....: p1 = p.polynomial()(x+1,y+1) + R.O(prec) ....: print 'prec=',prec ....: %timeit p^10 ....: %timeit p1^10 ....: prec= 10 625 loops, best of 3: 752 µs per loop 125 loops, best of 3: 3.77 ms per loop prec= 20 125 loops, best of 3: 3.4 ms per loop 25 loops, best of 3: 25.2 ms per loop prec= 30 125 loops, best of 3: 5.14 ms per loop 5 loops, best of 3: 151 ms per loop prec= 40 25 loops, best of 3: 13.3 ms per loop 5 loops, best of 3: 484 ms per loop
comment:73 follow-up: ↓ 74 Changed 8 years ago by
This is a follow-up of my previous mail on random_element in multivariate series.
R.random_element
calls R.__poly_ring.random_element
.
From the documentation of MPolynomialRing_generic.random_element
:
`First monomials are chosen uniformly random from the set of all
possible monomials of degree up to d
(inclusive). This means
that it is more likely that a monomial of degree d
appears than
a monomial of degree d-1
because the former class is bigger.`
This argument is relevant also for multivariate series
sage: R.<a,b,c,d> = PowerSeriesRing(QQ,4) sage: p = (1 + a + 2*b + 3*c + 4*d + R(0).O(16))^-5 sage: [len(x.exponents()) for x in p._bg_value] [1, 4, 10, 20, 35, 56, 84, 120, 165, 220, 286, 364, 455, 560, 680, 816]
There are many more terms with high degrees for the reason given in the quote.
Take a random series of comparable length
sage: p1 = R.random_element(16,4000) sage: [len(x.exponents()) for x in p1._bg_value] [0, 4, 3, 13, 22, 44, 54, 75, 114, 150, 190, 236, 307, 371, 440, 554]
The distribution is fairly similar.
However taking bound
small in R.random_element(prec,bound)
the lowest total degree is high
sage: p1 = R.random_element(16) sage: [len(x.exponents()) for x in p1._bg_value] [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1] sage: p1 = R.random_element(16,16) sage: [len(x.exponents()) for x in p1._bg_value] [0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 2, 0, 1, 3, 0, 2]
With a high number of variables the situation gets worse;
to get series starting with small lowest total degree
bound
should be extremely high
sage: R = PowerSeriesRing(QQ,100,'t') sage: p = R.random_element(10,100); sum(p.exponents()[0]) 9 sage: p = R.random_element(10,1000); sum(p.exponents()[0]) 7 sage: p = R.random_element(10,10000); sum(p.exponents()[0]) 7
In power series with finite precision prec
lower degrees are more important than higher degrees;
if this were not so it would not make sense to throw away terms
with degree higher than prec-1. So I think that lower degrees
should weight more than high degrees in the choice of random
elements.
comment:74 in reply to: ↑ 73 Changed 8 years ago by
Replying to pernici:
This is a follow-up of my previous mail on random_element in multivariate series.
Hi pernici,
Thanks for clarifying this, but I'm not sure how to proceed. I think I understand what you mean about lower degree terms being "more important" (in the same way that earlier decimal digits are "more important" in a real number). But I also think it makes a lot of sense for random_element
to be very similar between multivariate polynomials and multivariate power series -- this is the principle behind a number of implementation choices for multivariate power series. Nowhere is it promised that random_element
s are chosen with uniform distribution and, as you point out, the docstring for multivariate polynomials says explicitly that it is using a non-uniform distribution.
So now I wonder what the motivation is for changing the power series random_element
method: am I right to guess that you want to use the random elements to test your multiplication code? If so, let me suggest that it would be better to do this with a test suite, where you can write different custom algorithms for different random distributions and not have to worry about which of these should used for the random_element
method of polynomials or power series. (Having said that, I admit that I've found it difficult to get started with sage's test suite classes. I think there are some implementations for general rings and for integrals, so these might be good examples to follow if you want. )
Or maybe you have a different reason to think about random_element
?
best, Niles
comment:75 Changed 8 years ago by
Hi Niles, I agree that it has some advantages to keep the same random_element as for polynomials; at least it is well defined mathematically in terms of sets of monomials.
am I right to guess that you want to use the random elements to test your multiplication code?
Yes, my motivation for using random polynomials was for benchmarking multiplication; when I found out that random_element generates typically series without low order terms I replaced it with a generator which gave them.
Or maybe you have a different reason to think about random_element?
No.
Mario
comment:76 Changed 8 years ago by
In this ticket series with total degree truncation are considered.
How does the case of partial degree truncation work,
say with p
a polynomial
f = p + O(x^prec_x,y^prec_y)
Can one define the partial degree truncation as the total degree truncation with
prec = prec1+prec2-1
, and if f._bg_value
has valuation 0
truncate
x^prec_x and y^prec_y
?
With such a definition multiplication would work as in the total degree case, with a truncation of the result, if needed.
comment:77 follow-up: ↓ 78 Changed 8 years ago by
Please disregard my previous message :(
In this ticket series with total degree truncation are considered.
I would like to ask how the case of series with partial degree truncation works.
Consider for example the product of two series in variables x,y
with precisions prec_x,prec_y
; if they start with a constant term one can
just truncate with O(x^prec_x) , O(y^prec_y)
.
I do not see how it works if there are no constant terms; consider for example
two series with precision 3 in both variables; add a variable b
to track
the monomials which should be eliminated
sage: R.<x,y,b> = QQ[] sage: p1 = x + y + x^2 + x*y + y^2 + b*(x^3 + y^3) sage: p2 = x^2 + x*y + y^2 + b*(x^3 + y^3) sage: p3 = p1*p2 sage: p3.coefficient({b:0}) x^4 + 2*x^3*y + 3*x^2*y^2 + 2*x*y^3 + y^4 + x^3 + 2*x^2*y + 2*x*y^2 + y^3 sage: p3.coefficient({b:1}) 2*x^5 + 2*x^4*y + 2*x^3*y^2 + 2*x^2*y^3 + 2*x*y^4 + 2*y^5 + x^4 + x^3*y + x*y^3 + y^4
x^4, x^3*y, x*y^3, y^4
appear in p3.coefficient({b:1})
so they should be disregarded; the product of the series should be
3*x^2*y^2 + x^3 + 2*x^2*y + 2*x*y^2 + y^3
but I do not see how to get this simply with truncations;
neglecting 3*x^2*y^2
one would get the total degree truncation
sage: R.<x,y> = QQ[[]] sage: p1 = x + y + x^2 + x*y + y^2 + R.O(3) sage: p2 = x^2 + x*y + y^2 + R.O(3) sage: p1*p2 x^3 + 2*x^2*y + 2*x*y^2 + y^3 + O(x, y)^4
comment:78 in reply to: ↑ 77 Changed 8 years ago by
Replying to pernici:
... but I do not see how to get this simply with truncations; ...
I haven't been able to see how to do it either. There have been a number of unsuccessful tries to get multivariate series in sage, and I think the problem of separate variable precision has been one of the major barriers...
Changed 8 years ago by
comment:79 follow-ups: ↓ 81 ↓ 104 Changed 8 years ago by
attached patch to be applied after applying trac_1956_multi_power_series_new_4.patch, trac_1956_uni_multi_ps_2.patch, trac_1956_multi_ps_cleanup.patch
The univariate series has the property that its representation can be used to define it
sage: S.<t> = QQ[[]] sage: 1 + t + t^3 + O(t^5) 1 + t + t^3 + O(t^5)
In the case of multivariate series this property does not currently hold
sage: R.<x,y> = QQ[[]] sage: p = 1 + x + y^2 + R.O(5); p 1 + x + y^2 + O(x, y)^5 sage: 1 + x + y^2 + O(x, y)^5 TypeError
The attached patch is a hack to satisfy this property
sage: R.<x,y> = QQ[[]] sage: 1 + x + y^2 + O(x, y)^5 1 + x + y^2 + O(x, y)^5
comment:80 Changed 8 years ago by
- Description modified (diff)
- Summary changed from implement multivariate power series arithmetic to implement multivariate truncated power series arithmetic
comment:81 in reply to: ↑ 79 Changed 8 years ago by
Replying to pernici:
The attached patch is a hack to satisfy this property
sage: R.<x,y> = QQ[[]] sage: 1 + x + y^2 + O(x, y)^5 1 + x + y^2 + O(x, y)^5
Very nice syntax!
Btw: I allowed myself to edit the description and title of this ticket, since it took me a bit of time to figure out from the discussion which kind series it was about. That's a cool feature; thanks to all of you for working on it!
comment:82 Changed 8 years ago by
I get some strange behaviour; I want to use the reduce method for a polynomial extracted from a multivariate series
sage: S.<x,y> = QQ[[]] sage: p = (1+x+y)^3 sage: pp = p.polynomial() sage: R = S._poly_ring() sage: xv = R.gens() sage: pp.reduce([xv[0]^2,xv[1]^2]) 1
while the expected behavior is as working directly with polynomials
sage: R.<x,y> = QQ[] sage: p = (1+x+y)^3 sage: p.reduce([x^2,y^2]) 6*x*y + 3*x + 3*y + 1
After some trials I got the expected result by creating a new ring and casting the polynomial with it.
sage: S.<x,y> = QQ[[]] sage: p = (1+x+y)^3 sage: pp = p.polynomial() sage: R = PolynomialRing(QQ,S.variable_names()) sage: xv = R.gens() sage: pp = R(pp) sage: pp.reduce([xv[0]^2,xv[1]^2]) 6*x*y + 3*x + 3*y + 1
This is unsatisfactory because it should not be necessary to introduce a new MPolynomial_libsingular ring, and cast the polynomial.
Am I doing something wrong or is there a bug?
comment:83 Changed 8 years ago by
pernici wrote:
while the expected behavior is as working directly with polynomials
Sorry, there is no problem; the default ordering for multivariate series in negdeglex
,
while the default ordering of PolynomialRing
is degrevlex
, so that
reduce
acts differently in the two cases.
comment:84 Changed 8 years ago by
The method sqrt of MPowerSeries does not work because sqrt does not currently work for series on polynomial rings.
In ticket #10720 nth_root is added; applying trac_10720_power_series_nth_root_2.patch
the definition of nth_root for MPowerSeries is
sage: def nth_root(self, n): ....: return self.parent(self._bg_value.nth_root(n)) ....:
Examples:
sage: R.<x,y> = QQ[[]] sage: p = 1 + 2*x + 3*y + R.O(20) sage: nth_root(p^10,2) - p^5 0 + O(x, y)^20 sage: nth_root(p^10,5) - p^2 0 + O(x, y)^20 sage: nth_root(p^-10,2) - p^-5 0 + O(x, y)^20 sage: nth_root(p^-10,5) - p^-2 0 + O(x, y)^20
comment:85 Changed 8 years ago by
- Description modified (diff)
- Work issues multivariate series on 1 generator should remain different from a univariate series deleted
Hello all,
Here is a more refined attempt to distribute the review of this patch. If you can give any one of the following items positive review, please do so!
thanks, Niles
- Sage passes all doctests with this patch (buildbot gives this positive review)
- All code is documented and doctested thoroughly; documentation builds without error or warning
- The underlying concept of the implementation (a wrapper for certain univariate power series over multivariate polynomials) is sound
- multi_power_series_ring.py: the code accurately does what it claims to do
- multi_power_series_ring_element.py: the code accurately does what it claims to do
- Integration with the rest of sage: construction and use of PowerSeriesRings? works correctly, and parallels behavior of polynomial rings
- Performance: the multivariate power series arithmetic is fast enough to be included in Sage
- Coding: the code is free from obvious inefficiencies in error handling, memory management, etc.
- The items on this list constitute a complete review
For the buildbot:
Apply trac_1956_multi_power_series_new_4.patch, trac_1956_uni_multi_ps_2.patch, trac_1956_multi_ps_cleanup.patch
comment:86 follow-up: ↓ 87 Changed 8 years ago by
Niles, going back to the benchmarks you posted in ticket #1956 (see also bench1.sage posted here), here is a comparison with PARI
(2) latest patch in this ticket (3) PARI/GP times in ms test no. 1 2 3 4 5 6 7 8 (2) 28.2 2.3 24.8 166 8.16 21 45.2 56.7 (3) 27.8 20 19.7 94.4 16 43.2 47.2 48
PARI ranges from 8.7x slower to 1.8x faster; it is slowest in the second example; the timings for (3) fluctuate a lot, unlike those of (2); above I gave the first timings I got.
I took a few times the timings for (3) in the second example
sage: %timeit gp('(1 + T^3*(a^3 + b^3) + T^4*(c^4 + d^4) + O(T^(15)))^-20') 25 loops, best of 3: 17.3 ms per loop
I got timings from 14 ms to 23 ms, while I get regularly 2.3 ms for (2)
sage: L.<a,b,c,d> = PowerSeriesRing(QQ,4) sage: %timeit (1 + a^3 + b^3 + c^4 + d^4 + L.O(15))^-20 125 loops, best of 3: 2.3 ms per loop
even taking the best time for PARI in this example, it is still 6x slower.
I decided to make comparison with PARI/GP prompted by your post in ticket #1956 about distributing the review of that ticket.
I thought that another item there could be:
comparison with other CAS
and so I tried the comparison with PARI.
IMHO this ticket could be closed: I consider settled that generic truncated multiplication is faster than Karatsuba for multivariate series. The comparison with PARI is fairly favorable; the latest patch is from 6x faster to 2.5x slower than PARI/GP in the given benchmarks with the QQ field, and it is generic code, working with any ring.
These benchmarks should help reviewers in settling the item:
Performance: the multivariate power series arithmetic is fast enough to be included in Sage
comment:87 in reply to: ↑ 86 ; follow-up: ↓ 92 Changed 8 years ago by
Replying to pernici:
I thought that another item there could be:
comparison with other CAS
and so I tried the comparison with PARI.
Thanks for doing this. A comparison with Magma also appears in comment 8 above.
IMHO this ticket could be closed: I consider settled that generic truncated multiplication is faster than Karatsuba for multivariate series. The comparison with PARI is fairly favorable; the latest patch is from 6x faster to 2.5x slower than PARI/GP in the given benchmarks with the QQ field, and it is generic code, working with any ring.
Thanks! Really we're just waiting for someone (other than me) to click "positive review" -- this will initiate the process of merging it in the next version of Sage.
These benchmarks should help reviewers in settling the item:
Performance: the multivariate power series arithmetic is fast enough to be included in Sage
It sounds like you're willing to give that item a positive review. Are you also willing to say that "The underlying concept is sound"? If there are any other items you can say should have positive review, please do :)
comment:88 follow-up: ↓ 90 Changed 8 years ago by
niles wrote:
It sounds like you're willing to give that item a positive review. Are you also willing to say that "The underlying concept is sound"?
I think it is sound to use a background univariate series to implement total degree precision.
Concerning the 6th item (integration with the rest of sage) we should discuss what other approaches are available in Sage to deal with multivariate series.
Another approach to multivariate series is using series of series;
this can be done with Sage's univariate power series or gp
.
It is a bit cumbersome to convert a polynomial into a series of series with Sage's univariate power series
sage: R.<x1> = QQ[[]] sage: S.<x0> = R[[]] sage: p = (1 + x0 + x1)^3 sage: a = p.list() sage: q = sum((a[i] + O(x1^3))*x0^i for i in range(len(a))) + O(x0^3) sage: q^-5 1 - 15*x1 + 120*x1^2 + O(x1^3) + (-15 + 240*x1 - 2040*x1^2 + O(x1^3))*x0 + (120 - 2040*x1 + 18360*x1^2 + O(x1^3))*x0^2 + O(x0^3) sage: (q-1)^2 9*x1^2 + 18*x1^3 + O(x1^4) + (18*x1 + 54*x1^2 + O(x1^3))*x0 + (9 + 54*x1 + 90*x1^2 + O(x1^3))*x0^2 + O(x0^3)
Using gp
the notation is simpler
sage: q = gp('(1+x0+x1)^3 + O(x0^3) + O(x1^3)')
giving the same output for q^-5
and (q-1)^2
.
For long computations gp
is much faster, especially with more than two variables.
Should one integrate the multi-series multivariate series approach
with PowerSeriesRing
? It could be something like
sage: R = PowerSeriesRing(QQ, 2, 'x', "ms") sage: q = (1 + x0 + x1)^3 + O(x0^3,x1^3) sage: q 1 + 3*x1 + 3*x1^2 + O(x1^3) + (3 + 6*x1 + 3*x1^2 + O(x1^3))*x0 + (3 + 3*x1 + O(x1^3))*x0^2 + O(x0^3)
One could implement series of series with Sage's univariate power series, or better using PARI.
Clearly implementing this does not belong to this ticket.
The only thing to be decided in this ticket is if it is OK to keep the
total degree precision as the default option in PowerSeriesRing
.
comment:89 Changed 8 years ago by
In ticket #10532 the patch trac_10532_send_to_bg.patch fixes the bug in _send_to_bg in the case of multivariate series with one variable, see comments 64-70 in this ticket.
comment:90 in reply to: ↑ 88 Changed 8 years ago by
Replying to pernici:
niles wrote:
It sounds like you're willing to give that item a positive review. Are you also willing to say that "The underlying concept is sound"?
I think it is sound to use a background univariate series to implement total degree precision.
Thanks :)
Concerning the 6th item (integration with the rest of sage) we should discuss what other approaches are available in Sage to deal with multivariate series.
Another approach to multivariate series is using series of series; this can be done with Sage's univariate power series or
gp
. ... The only thing to be decided in this ticket is if it is OK to keep the total degree precision as the default option inPowerSeriesRing
.
With regard to "series of series", I believe there are some theoretical subtleties which make this tricky. The Pari User's Manual contains the following warning (1.2.4, page 22):
In the present version 2.3.3 of PARI, there are bugs in the handling of power series of power series, i.e. power series in several variables. ... This bug is difficult to correct because the mathematical problem itself contains some amount of imprecision, and it is not easy to design an intuitive generic interface for such beasts.
So I think total-degree precision should be the the default, until someone writes code to carefully handle the other possibilities (or to call gp
-- see the next point).
For long computations gp is much faster, especially with more than two variables.
Yes, a more complete interface with gp
is probably the right way to go. However I think PowerSeriesRing
, and its elements, should be native Sage objects -- this is not the case with gp
series, which makes them unsuitable for the default power series. For example:
sage: q = gp('(1+x0+x1)^3 + O(x0^3) + O(x1^3)') sage: type(q) <class 'sage.interfaces.gp.GpElement'> sage: q.parent() PARI/GP interpreter sage: q.base_ring() PARI/GP interpreter
comment:91 Changed 8 years ago by
niles wrote:
So I think total-degree precision should be the the default, until someone writes code to carefully handle the other possibilities (or to call gp -- see the next point).
Maybe total-degree precision is a better default anyway, even if there were available an excellent multi-series (I mean series of series ... N times, I do not know the standard name for these _really_ multivariate series) implementation; I have the impression that few times multi-series are used, and that univariate series (which total-degree precision multivariate series really is) or at most bivariate series are used frequently in practice.
Yes, a more complete interface with gp is probably the right way to go. However I think PowerSeriesRing??, and its elements, should be native Sage objects -- this is not the case with gp series, which makes them unsuitable for the default power series.
Right.
comment:92 in reply to: ↑ 87 Changed 8 years ago by
Replying to niles:
Replying to pernici:
I thought that another item there could be:
comparison with other CAS
and so I tried the comparison with PARI.
Thanks for doing this. A comparison with Magma also appears in comment 8 above.
I do not have Magma; according to the timings in comment 8 it seems that in the first benchmark Magma is 3 orders of magnitude slower.
I found this http://wiki.sagemath.org/sagebeatsmagma ; such a benchmark could be added there.
comment:93 Changed 8 years ago by
- Description modified (diff)
More code rot; patches should all apply to 4.6.2 now.
For the buildbot:
Apply trac_1956_multi_power_series_new_4.patch, trac_1956_uni_multi_ps_3.patch, trac_1956_multi_ps_cleanup.patch
comment:94 Changed 8 years ago by
correction to buildbot:
Apply trac_1956_multi_power_series_new_4.patch, trac_1956_uni_multi_ps_3.patch, trac_1956_multi_ps_cleanup.patch
comment:95 Changed 8 years ago by
updated trac_1956_multi_power_series_new_4.patch to include #random
flag in doctests of .random_element()
(lines 291 and 294 of this patch). Although I've tried to avoid #random
throughout, it's necessary here unless I explicitly set the random seed before the tests. However, even then, the result of the test may change from one version of Sage to the next (as it has now).
Buildbot:
Apply trac_1956_multi_power_series_new_4.patch, trac_1956_uni_multi_ps_3.patch, trac_1956_multi_ps_cleanup.patch
comment:96 follow-up: ↓ 97 Changed 8 years ago by
- Cc hlaw added
First, thanks for all your work Niles, what you've done has already saved me hours of pain!
There seems to be a problem with how coercion/evaluation works between two multivariate power series rings when the number of variables in each differs. For example:
sage: A.<a, b, c> = PowerSeriesRing(QQ) sage: B.<s, t> = PowerSeriesRing(QQ) sage: s(a, b) --------------------------------------------------------------------------- TypeError Traceback (most recent call last) /home/hlaw/Dropbox/org/<ipython console> in <module>() /home/hlaw/src/sage-src/local/lib/python2.6/site-packages/sage/rings/multi_power_series_ring_element.pyc in __call__(self, *x, **kwds) 419 for i in range(len(x)): 420 try: --> 421 xi = self.parent(x[i]) 422 except AttributeError: 423 raise AttributeError(str(x[i])+" does not coerce to parent ring.") /home/hlaw/src/sage-src/local/lib/python2.6/site-packages/sage/structure/element.so in sage.structure.element.Element.parent (sage/structure/element.c:4069)() /home/hlaw/src/sage-src/local/lib/python2.6/site-packages/sage/rings/power_series_ring.pyc in __call__(self, f, prec, check) 604 v = sage_eval(f.Eltseq()) 605 return self(v) * (self.gen(0)**f.Valuation()) --> 606 return self.__power_series_class(self, f, prec, check=check) 607 608 def construction(self): /home/hlaw/src/sage-src/local/lib/python2.6/site-packages/sage/rings/multi_power_series_ring_element.pyc in __init__(self, parent, x, prec, is_gen, check) 352 self._bg_value = parent._bg_ps_ring(x._bg_value) 353 except TypeError: --> 354 raise TypeError("Unable to coerce into background ring.") 355 356 elif xparent == parent._bg_ps_ring(): TypeError: Unable to coerce into background ring.
I would expect this to print a
. Similar permutations with s^2
, s + t
, etc. produce the same result.
Another problem, possibly related, is the following:
sage: B.<s, t> = PowerSeriesRing(QQ) sage: C.<z> = PowerSeriesRing(QQ) sage: s(z, z) # expect z 1 sage: s(z, z^2) # expect z 1 sage: s(0, z^2) # expect 0 0 sage: s(z^2, 0) # expect z^2 1 sage: (2*s^2)(z, z) # expect 2*z^2 2 sage: (2*s^2 + 3*t)(z, z) # expect 3*z + 2*z^2 5 sage: t(0, z) # expect z 1 sage: t(z, z) # expect z 1 sage: z(s + t) # expect s + t s + t sage: z(s) # expect s s sage: z(t) # expect t t
Clearly z
is being sent to 1 (one) somehow when going from QQ[[s,t]]
to QQ[[z]]
, but the maps in the other direction work fine.
Both of these problems (assuming they're distinct) disappear when one does the analogous calculations directly with multivariate polynomial rings.
Cheers, Hamish.
(Using Sage 4.6.2 +power series patches, on Ubuntu 10.10 64bit, Macbook 5.1.)
comment:97 in reply to: ↑ 96 ; follow-up: ↓ 98 Changed 8 years ago by
- Status changed from needs_review to needs_work
Replying to hlaw:
First, thanks for all your work Niles, what you've done has already saved me hours of pain!
Thanks -- glad to hear it! It has cost me many hours, so I hope the review will be finished soon (after I fix this problem :)
There seems to be a problem with how coercion/evaluation works between two multivariate power series rings when the number of variables in each differs.
Actually it never occurred to me that people would want to substitute into a power series from one ring elements from a different ring, or that such behavior should be supported, but I can see now why it might be useful. The problem is occurring in the __call__
method of MPowerSeries
, which assumes that the inputs should coerce to the parent of the power series they're being substituted into. A separate part of the problem is that elements from one power series ring are incorrectly coercing to 1 in the other power series ring:
sage: B.<s, t> = PowerSeriesRing(QQ) sage: C.<z> = PowerSeriesRing(QQ) sage: B(z) 1 sage: C(s) --------------------------------------------------------------------------- ... TypeError: Unable to coerce s (<class 'sage.rings.multi_power_series_ring_element.MPowerSeries'>) to Rational
So I will have to figure out what is going wrong here. In the meantime, if you need to actually use this, you can use ._value()
to recover the underlying multivariate polynomial, do the substitution there, and then .add_bigoh()
to get a power series again. The important idea is that the precision of the power series after substitution should be n*v
, where n
is the precision of the original power series, and v
is the minimum of the valuations of the input power series.
Something smarter might be done by mimicking the __call__
method of MPolynomial_element
, but I haven't thought through that yet -- let me know if you do :)
comment:98 in reply to: ↑ 97 ; follow-up: ↓ 101 Changed 8 years ago by
Replying to niles:
Replying to hlaw:
There seems to be a problem with how coercion/evaluation works between two multivariate power series rings when the number of variables in each differs.
I tracked down this bug -- it was unique to the case where one of the power series rings is univariate, and comes from the fact that (univariate) power series rings automatically convert any univariate power series from a given variable to another:
sage: R.<a> = PowerSeriesRing(ZZ)sage: S.<b> = PowerSeriesRing(ZZ) sage: b in R False sage: R.has_coerce_map_from(S) False sage: R(b) a
I was thinking this should be filed as a separate bug, but polynomials have the same behavior, so perhaps it is intentional. In any case, I've fixed multivariate power series to take this into account, so an error is raised when the variables don't match:
sage: B.<s, t> = PowerSeriesRing(QQ) sage: C.<z> = PowerSeriesRing(QQ) sage: B(z) ... TypeError: Cannot coerce input to polynomial ring.
Another problem, possibly related, is the following...
I also implemented a ._subs_formal
method, which is called automatically by .__call__
. This method substitutes *anything* into the power series, requiring only that multiplication and exponentiation be defined for the inputs. This gives all of the "expected" results from above. Some further examples (from the new doctests):
sage: B.<s, t> = PowerSeriesRing(QQ) sage: C.<z> = PowerSeriesRing(QQ) sage: s(z,z) z sage: f = -2/33*s*t^2 - 1/5*t^5 - s^5*t + s^2*t^4 sage: f(z,z) -2/33*z^3 - 1/5*z^5 sage: f(z,1) -1/5 - 2/33*z + z^2 - z^5 sage: RF = RealField(10) sage: f(z,RF(1)) -0.20 - 0.061*z + 1.0*z^2 - 0.00*z^3 - 0.00*z^4 - 1.0*z^5 sage: m = matrix(QQ,[[1,0,1],[0,2,1],[-1,0,0]]) sage: m [ 1 0 1] [ 0 2 1] [-1 0 0] sage: f(m,m) [ 2/33 0 1/5] [ 131/55 -1136/165 -24/11] [ -1/5 0 -23/165] sage: f(m,m) == -2/33*m^3 - 1/5*m^5 True sage: f = f.add_bigoh(10) sage: f(z,z) -2/33*z^3 - 1/5*z^5 + O(z^10) sage: f(m,m) # cannot substitute arbitrary elements when f has finite precision Traceback (most recent call last): ... AttributeError: 'sage.matrix.matrix_rational_dense.Matrix_rational_dense' object has no attribute 'add_bigoh'
Please test it out and let me know how it goes! -Niles
comment:99 Changed 8 years ago by
- Description modified (diff)
attached the patch which implements the changes described above; now updating "Apply" instructions
Buildbot:
Apply trac_1956_multi_power_series_new_4.patch, trac_1956_uni_multi_ps_3.patch, trac_1956_multi_ps_cleanup.patch, trac_1956_multi_ps_coercion.patch
comment:100 Changed 8 years ago by
apologies for the many rapid updates -- fixing some stupid typos and trying to get the buildbot to apply the right patches.
Apply trac_1956_multi_power_series_new_4.patch, trac_1956_uni_multi_ps_3.patch, trac_1956_multi_ps_cleanup.patch, trac_1956_multi_ps_coercion.patch
comment:101 in reply to: ↑ 98 ; follow-up: ↓ 102 Changed 8 years ago by
Replying to niles:
Please test it out and let me know how it goes! -Niles
Works great! The quick turnaround is appreciated!
Cheers, Hamish.
comment:102 in reply to: ↑ 101 Changed 8 years ago by
- Status changed from needs_work to needs_review
Replying to hlaw:
Works great! The quick turnaround is appreciated!
Glad to hear it! If you feel capable of reviewing some part of the patch, that would be *greatly* appreciated. Here is the current status of the various items:
- Sage passes all doctests with this patch
(positive review; buildbot)
- All code is documented and doctested thoroughly; documentation builds without error or warning
(needs review)
- The underlying concept of the implementation (a wrapper for certain univariate power series over multivariate polynomials) is sound
(positive review; Mario Pernici)
- multi_power_series_ring.py: the code accurately does what it claims to do
(needs review)
- multi_power_series_ring_element.py: the code accurately does what it claims to do
(needs review)
- Integration with the rest of sage: construction and use of
PowerSeriesRings
works correctly, and parallels behavior of polynomial rings
(needs review)
- Performance: the multivariate power series arithmetic is fast enough to be included in Sage
(positive review; Mario Pernici)
- Coding: the code is free from obvious inefficiencies in error handling, memory management, etc.
(needs review)
- The items on this list constitute a complete review
(needs review)
comment:103 Changed 8 years ago by
- Description modified (diff)
I've combined all of the partial patches into a single patch. There are still a number of items needing review.
If you've suggested improvements to this patch in the past, please verify that the issues have been resolved, and mark "positive review" for anything you can!
Much thanks, Niles
Buildbot: Apply trac_1956_combined.patch
comment:104 in reply to: ↑ 79 Changed 8 years ago by
comment:105 follow-up: ↓ 108 Changed 8 years ago by
- Status changed from needs_review to needs_work
- Sage passes all doctests with this patch (positive review; buildbot)
I think the buildbot gets confused by the "Apply" instruction in the ticket's description. Hence, it didn't run the big patch yet.
- All code is documented and doctested thoroughly; documentation builds without
error or warning (needs review)
needs work
sage -coverage devel/sage/sage/rings/multi_power_series_ring.py ---------------------------------------------------------------------- devel/sage/sage/rings/multi_power_series_ring.py SCORE devel/sage/sage/rings/multi_power_series_ring.py: 96% (30 of 31) Missing documentation: * unpickle_multi_power_series_ring_v0(base_ring, num_gens, names, order, default_prec, sparse):
These just need #indirect doctest
added:
sage -coverage devel/sage/sage/rings/multi_power_series_ring_element.py ---------------------------------------------------------------------- devel/sage/sage/rings/multi_power_series_ring_element.py SCORE devel/sage/sage/rings/multi_power_series_ring_element.py: 100% (52 of 52) Possibly wrong (function name doesn't occur in doctests): * _subs_formal(self,*x,**kwds): * _add_(left, right): * _sub_(left, right): * _mul_(left, right): * _lmul_(self, c): * _div_(self, denom_r): ----------------------------------------------------------------------
But the docs build without warnings etc.
- The underlying concept of the implementation (a wrapper for certain univariate
power series over multivariate polynomials) is sound (positive review; Mario Pernici)
- multi_power_series_ring.py: the code accurately does what it claims to do (needs_review)
Shouldn't def _element_constructor_(self,f)
use the default precision of the ring?
- multi_power_series_ring_element.py: the code accurately does what it claims to
do (needs review)
- Integration with the rest of sage: construction and use of PowerSeriesRings?
works correctly, and parallels behavior of polynomial rings (needs review)
needs work
The default_prec
parameter does not seem to work as expected (or is my expectation wrong?)
sage: P.<x,y> = PowerSeriesRing(GF(127),default_prec=20) sage: x.prec() +Infinity
- Performance: the multivariate power series arithmetic is fast enough to be
included in Sage (positive review; Mario Pernici)
- Coding: the code is free from obvious inefficiencies in error handling, memory
management, etc. (needs review)
If the performance is good (7), then I think we can just assume it's okay. We can always come back and fix stuff later which is too slow.
I have to say that I don't like the excessive use of double underscore attributes. It makes it harder to inherit from it or reach into the internals (just as you have to write _PowerSeriesRing_generic__poly_ring
) I'd just use a single underscore. I'll leave it up to you (i.e., not make my review dependent on it) whether you insist on keeping them.
- The items on this list constitute a complete review (needs review)
Fine by me.
Changed 8 years ago by
comment:106 follow-up: ↓ 107 Changed 8 years ago by
- Description modified (diff)
comment:107 in reply to: ↑ 106 Changed 8 years ago by
Replying to niles:
Patchbot: Apply trac_1956_combined_2.patch
comment:108 in reply to: ↑ 105 ; follow-ups: ↓ 109 ↓ 110 Changed 8 years ago by
- Status changed from needs_work to needs_review
Replying to malb:
Thanks for taking a look; sorry for the coverage negligence. I've added the #indirect notes
, and a docstring for unpickle
. Now sage -coverage
returns 100% for both multi_power_series
files. As a bonus, I've added a docstring for the unpickle function of univariate power series too.
I believe your other comments do not require changes -- I've elaborated below. This means we are once again ready for review :)
- multi_power_series_ring.py: the code accurately does what it claims to do (needs_review)
Shouldn't
def _element_constructor_(self,f)
use the default precision of the ring?
No, I believe not. Default precision is used only when absolutely necessary to convert an infinite precision element to a finite precision one (e.g. for inversion or reversion). But, for example, coercion of a polynomial should result in an infinite precision element, not a finite precision one. This is consistent with the behavior of univariate power series rings:
sage: P = QQ[x]; P Univariate Polynomial Ring in x over Rational Field sage: S = P.completion(P.gen()); S Power Series Ring in x over Rational Field sage: S(P.random_element()).prec() +Infinity
sage: P = QQ['x,y,z']; P Multivariate Polynomial Ring in x, y, z over Rational Field sage: S = P.completion(P.gens()); S Multivariate Power Series Ring in x, y, z over Rational Field sage: S(P.random_element()).prec() +Infinity
- Integration with the rest of sage: construction and use of
PowerSeriesRing
works correctly, and parallels behavior of polynomial rings (needs review)
needs work
The
default_prec
parameter does not seem to work as expected (or is my expectation wrong?)
sage: P.<x,y> = PowerSeriesRing(GF(127),default_prec=20) sage: x.prec() +Infinity
I believe that your expectation is wrong, as described above. Note the parallel behavior for univariate power series rings:
sage: P.<x> = PowerSeriesRing(GF(127),default_prec=20) sage: x.prec() +Infinity
- Coding: the code is free from obvious inefficiencies in error handling, memory
management, etc. (needs review)
If the performance is good (7), then I think we can just assume it's okay. We can always come back and fix stuff later which is too slow.
so...I'm reading 'positive review', right? ;) (free from *obvious* inefficiencies, etc.)
I have to say that I don't like the excessive use of double underscore attributes. It makes it harder to inherit from it or reach into the internals (just as you have to write
_PowerSeriesRing_generic__poly_ring
) I'd just use a single underscore. I'll leave it up to you (i.e., not make my review dependent on it) whether you insist on keeping them.
Excessive or not, each of the double underscore methods is necessary to replace the corresponding method of PowerSeries
or PowerSeriesRing_generic
, from which the MPowerSeries*
classes inherit. Here's a complete list of the double underscore attributes:
MPowerSeries.__init__ MPowerSeries.__reduce__ MPowerSeries.__call__ MPowerSeries.__getitem__ MPowerSeries.__invert__ MPowerSeries.__cmp__ MPowerSeries.__mod__ MPowerSeries.__lshift__ MPowerSeries.__rshift__ MPowerSeriesRing_generic.__init__ MPowerSeriesRing_generic.__reduce__ MPowerSeriesRing_generic.__cmp__
comment:109 in reply to: ↑ 108 Changed 8 years ago by
Replying to niles:
Here's a complete list of the double underscore attributes...
To clarify: that's a complete list of such attributes defined in the code I've written. You will see many more (50 or 60) if you use tab-completion to see what double underscore methods a multivariate power series or multivariate power series ring has. I have not added any new ones (from the univariate case), simply redefined some.
comment:110 in reply to: ↑ 108 ; follow-up: ↓ 111 Changed 8 years ago by
This seems to be the current status:
- Sage passes all doctests with this patch
waiting for patchbot :)
- All code is documented and doctested thoroughly; documentation builds without error or warning
positive review
- The underlying concept of the implementation (a wrapper for certain univariate power series over multivariate polynomials) is sound
positive review
- multi_power_series_ring.py: the code accurately does what it claims to do
positive review
- multi_power_series_ring_element.py: the code accurately does what it claims to do
needs review
- Integration with the rest of sage: construction and use of PowerSeriesRings? works correctly, and parallels behavior of polynomial rings
positive review
- Performance: the multivariate power series arithmetic is fast enough to be included in Sage
positive review
- Coding: the code is free from obvious inefficiencies in error handling, memory management, etc.
positive review
- The items on this list constitute a complete review
positive review
Thanks for explaining the default precision business to me! As I said, I rarely touch power series.
Excessive or not, each of the double underscore methods is necessary to replace the corresponding method of
PowerSeries
orPowerSeriesRing_generic
, from which theMPowerSeries*
classes inherit. Here's a complete list of the double underscore attributes:
Sorry, I didn't mean stuff like __init__
which of course has to be like it is. I meant stuff like __poly_ring
which could easily be _poly_ring
.
comment:111 in reply to: ↑ 110 Changed 8 years ago by
Replying to malb:
This seems to be the current status:
...
Great! Thanks :)
Thanks for explaining the default precision business to me! As I said, I rarely touch power series.
Glad my explanation made sense :)
Excessive or not, each of the double underscore methods is necessary ...
Sorry, I didn't mean stuff like
__init__
which of course has to be like it is. I meant stuff like__poly_ring
which could easily be_poly_ring
.
Ah, I agree. This came up once before (perhaps with you, or with Simon), and indeed I've eliminated the unnecessary double underscore attributes. Those listed above are the only ones remaining.
comment:112 follow-up: ↓ 114 Changed 8 years ago by
looking at trac_1956_combined_2.patch I see __ngens
, __term_order
and __poly_ring
, did you forget to upload the patch where you replaced those by _ngens
, _term_order
and _poly_ring
. Or are we talking about different things?
comment:113 Changed 8 years ago by
Shouldn't this be a value error instead of a type error? Since it's a wrong value and not type?
if len(x) != self.parent().ngens(): raise TypeError("Number of arguments does not match number of variables in parent.")
Other than that, I think the code in multi_power_series_ring_element.py
looks fine.
comment:114 in reply to: ↑ 112 ; follow-up: ↓ 115 Changed 8 years ago by
- Status changed from needs_review to needs_work
Replying to malb:
looking at trac_1956_combined_2.patch I see
__ngens
,__term_order
and__poly_ring
, did you forget to upload the patch where you replaced those by_ngens
,_term_order
and_poly_ring
. Or are we talking about different things?
Ack! No, but I was looking for double underscore attributes defined by def *
, totally forgetting about these others which are not methods. I'll upload a new version tomorrow.
comment:115 in reply to: ↑ 114 Changed 8 years ago by
...*sigh*... I wish I could update my previous comment to fix the syntax:
Replying to niles:
Replying to malb:
looking at trac_1956_combined_2.patch I see
__ngens
,__term_order
and__poly_ring
, did you forget to upload the patch where you replaced those by_ngens
,_term_order
and_poly_ring
. Or are we talking about different things?
I was looking for double underscore attributes defined by "def __*
", totally forgetting about these others which are not methods. I'll upload a new version tomorrow.
comment:116 Changed 8 years ago by
I'll make an honest attempt to look at it tomorrow even (GMT). So if by the day after tomorrow you haven't got a review, feel free to bug me about it!
comment:117 follow-up: ↓ 120 Changed 8 years ago by
Now that I've looked carefully at these, and read documents like PEP 8 style guide about double underscore attributes, I think I'm inclined to leave them as is. There are three contributing reasons: consistency with multivariate polynomial rings, consistency with univariate power series rings, and not wanting the background power series ring to be treated casually.
Here is the list of the double underscore attributes of MPowerSeriesRing_generic
:
- Those which parallel multivariate polynomials:
__term_order __ngens
- Those which parallel univariate power series:
__poly_ring __is_sparse __params
- Others, which I'd rather leave double-underscore for consistency with the rest and because I don't think the background ring should be treated casually:
__bg_power_series_ring __bg_indeterminate __base_ring
Changed 8 years ago by
comment:118 follow-up: ↓ 119 Changed 8 years ago by
- Description modified (diff)
- Status changed from needs_work to needs_review
Now I've added combined_3.patch
: It changes TypeError
to ValueError
as suggested above, and also in one other place (where again input of wrong length raises the error). The error types in the rest of the code seem good to me.
I've also changed the name of the indeterminate in the background power series ring from 'T'
to 'Tbg'
. It's conceivable that bugs could arise if a user happens to work with a multivariate power series ring, and also with a univariate power series ring having the same indeterminate name as that of the background ring of the multivariate ring.
Since I don't know exactly how these might come up, if ever, I'd like to let this be sufficient for now, and add a more robust system of protections against possible bugs in a separate ticket.
comment:119 in reply to: ↑ 118 Changed 8 years ago by
Patchbot: apply trac_1956_combined_3.patch
comment:120 in reply to: ↑ 117 Changed 8 years ago by
Replying to niles:
- Others, which I'd rather leave double-underscore for consistency with the rest and because I don't think the background ring should be treated casually:
__bg_power_series_ring __bg_indeterminate __base_ring
Why is there a __base_ring
at all? The base ring is usually taken care of by inheritance of sage.structure.parent.Parent
or other base classes. So, it should certainly be possible to work without a new attribute -- it seems to be a duplicate.
comment:121 follow-up: ↓ 122 Changed 8 years ago by
Hi, Simon good catch about the __base_ring
attribute.
@Niles: as you realised in your code (in __init__
of the ring), having double underscore attributes makes it hard to inherit because double underscore names are prefixed by the class name. If you want subclasses to use these attributes then you'll need to switch to a single underscore. The single underscore expresses that the attribute is meant to be private, i.e., not treated casually.
With double underscores there is no such thing as consistency, the point of them is to break consistency to prevent people from using them, if they are not in exactly the right class of the hierarchy. The design is broken if you have to set __poly_ring
yourself twice to make it work for the various classes in the inheritance hierarchy.
I agree that other rings do it wrong, too. But there's no need to repeat this wrong approach in a new class?
comment:122 in reply to: ↑ 121 Changed 8 years ago by
Replying to malb:
@Niles: as you realised in your code (in
__init__
of the ring), having double underscore attributes makes it hard to inherit because double underscore names are prefixed by the class name.
To be precise: Double underscore names are only problematic if they start with two underscores but do not end with two underscores.
comment:123 Changed 8 years ago by
- Description modified (diff)
Very well, I've been convinced -- thanks for being patient. The double underscore attributes mentioned above have all been replaced with single underscore versions. For __poly_ring
, I had to use _poly_ring_
, since there is a class method _poly_ring
already inherited from univariate power series.
All tests pass for the multi_power_series*
files; the rest should be fine so I'll wait for the patchbot to verify them.
comment:124 Changed 8 years ago by
Patchbot: apply trac_1956_combined_4.patch
comment:125 Changed 8 years ago by
- Status changed from needs_review to positive_review
Okay, let's push the button: positive review.
comment:126 Changed 8 years ago by
- Milestone changed from sage-4.7 to sage-4.7.1
comment:127 Changed 8 years ago by
- Status changed from positive_review to needs_work
sage -t -force_lib devel/sage/sage/rings/multi_power_series_ring.py ********************************************************************** File "/mnt/usb1/scratch/jdemeyer/merger/sage-4.7.1.alpha1/devel/sage-main/sage/rings/multi_power_series_ring.py", line 456: sage: (c,R) = M.construction(); (c,R) Expected: (CompletionFunctor, Multivariate Polynomial Ring in f0, f1, f2, f3 over Rational Field) Got: (Completion[('f0', 'f1', 'f2', 'f3')], Multivariate Polynomial Ring in f0, f1, f2, f3 over Rational Field) ********************************************************************** File "/mnt/usb1/scratch/jdemeyer/merger/sage-4.7.1.alpha1/devel/sage-main/sage/rings/multi_power_series_ring.py", line 459: sage: c Expected: CompletionFunctor Got: Completion[('f0', 'f1', 'f2', 'f3')] **********************************************************************
comment:128 Changed 8 years ago by
- Description modified (diff)
- Status changed from needs_work to needs_review
A new patch, updated according to the repr string of the completion functor. No other changes.
comment:129 Changed 8 years ago by
Patchbot: apply trac_1956_combined_5.patch
comment:130 follow-up: ↓ 131 Changed 8 years ago by
- Status changed from needs_review to needs_work
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 **********************************************************************
comment:131 in reply to: ↑ 130 Changed 8 years ago by
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?
Changed 8 years ago by
comment:132 Changed 8 years ago by
- Description modified (diff)
- Status changed from needs_work to needs_review
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
comment:133 follow-up: ↓ 134 Changed 8 years ago by
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.
comment:134 in reply to: ↑ 133 Changed 8 years ago by
comment:136 Changed 8 years ago by
- Merged in set to sage-4.7.1.alpha2
- Resolution set to fixed
- Status changed from positive_review to closed
Hi Mike,
hasn't something happened in that direction already or am I confused?
Cheers,
Michael