Opened 5 years ago
Closed 2 years ago
#19945 closed defect (fixed)
Fix Rational.__pow__
Reported by:  cheuberg  Owned by:  

Priority:  major  Milestone:  sage8.5 
Component:  basic arithmetic  Keywords:  
Cc:  dkrenn, behackl  Merged in:  
Authors:  Jeroen Demeyer  Reviewers:  Travis Scrimshaw 
Report Upstream:  N/A  Work issues:  
Branch:  5706029 (Commits)  Commit:  5706029d8199f2e7d0e66ba84f19a44f3336e04c 
Dependencies:  Stopgaps: 
Description (last modified by )
The following is unexpected:
sage: A.<n> = AsymptoticRing('QQ^n * n^QQ', ZZ); (1/2)^n  ValueError Traceback (most recent call last) <ipythoninput18784db7f7676d> in <module>() > 1 A = AsymptoticRing('QQ^n * n^QQ', ZZ, names=('n',)); (n,) = A._first_ngens(1); (Integer(1)/Integer(2))**n /usr/local/src/sageconfig/src/sage/rings/rational.pyx in sage.rings.rational.Rational.__pow__ (build/cythonized/sage/rings/rational.c:23638)() 2549 2550 if isinstance(n, Element): > 2551 return (<Element>n)._parent(self)**n 2552 try: 2553 return n.parent()(self)**n /usr/local/src/sageconfig/src/sage/structure/parent.pyx in sage.structure.parent.Parent.__call__ (build/cythonized/sage/structure/parent.c:9529)() 916 if mor is not None: 917 if no_extra_args: > 918 return mor._call_(x) 919 else: 920 return mor._call_with_args(x, args, kwds) /usr/local/src/sageconfig/src/sage/structure/coerce_maps.pyx in sage.structure.coerce_maps.DefaultConvertMap_unique._call_ (build/cythonized/sage/structure/coerce_maps.c:4978)() 153 print(type(C), C) 154 print(type(C._element_constructor), C._element_constructor) > 155 raise 156 157 cpdef Element _call_with_args(self, x, args=(), kwds={}): /usr/local/src/sageconfig/src/sage/structure/coerce_maps.pyx in sage.structure.coerce_maps.DefaultConvertMap_unique._call_ (build/cythonized/sage/structure/coerce_maps.c:4846)() 148 cdef Parent C = self._codomain 149 try: > 150 return C._element_constructor(x) 151 except Exception: 152 if print_warnings: /usr/local/src/sageconfig/local/lib/python2.7/sitepackages/sage/rings/asymptotic/asymptotic_ring.pyc in _element_constructor_(self, data, simplify, convert) 3891 return result 3892 > 3893 return self.create_summand('exact', data) 3894 3895 /usr/local/src/sageconfig/local/lib/python2.7/sitepackages/sage/rings/asymptotic/asymptotic_ring.pyc in create_summand(self, type, data, **kwds) 4335 4336 try: > 4337 return self(TM(data, **kwds), simplify=False, convert=False) 4338 except ZeroCoefficientError: 4339 return self.zero() /usr/local/src/sageconfig/src/sage/structure/parent.pyx in sage.structure.parent.Parent.__call__ (build/cythonized/sage/structure/parent.c:9529)() 916 if mor is not None: 917 if no_extra_args: > 918 return mor._call_(x) 919 else: 920 return mor._call_with_args(x, args, kwds) /usr/local/src/sageconfig/src/sage/structure/coerce_maps.pyx in sage.structure.coerce_maps.DefaultConvertMap_unique._call_ (build/cythonized/sage/structure/coerce_maps.c:4978)() 153 print(type(C), C) 154 print(type(C._element_constructor), C._element_constructor) > 155 raise 156 157 cpdef Element _call_with_args(self, x, args=(), kwds={}): /usr/local/src/sageconfig/src/sage/structure/coerce_maps.pyx in sage.structure.coerce_maps.DefaultConvertMap_unique._call_ (build/cythonized/sage/structure/coerce_maps.c:4846)() 148 cdef Parent C = self._codomain 149 try: > 150 return C._element_constructor(x) 151 except Exception: 152 if print_warnings: /usr/local/src/sageconfig/local/lib/python2.7/sitepackages/sage/rings/asymptotic/term_monoid.pyc in _element_constructor_(self, data, coefficient) 1670 except ValueError as e: 1671 raise combine_exceptions( > 1672 ValueError('%s is not in %s.' % (data, self)), e) 1673 1674 return self._create_element_(growth, coefficient) ValueError: 1/2 is not in Exact Term Monoid QQ^n * n^QQ with coefficients in Integer Ring. > *previous* ValueError: Factor 1/2 of 1/2 is neither a coefficient (in Integer Ring) nor growth (in Growth Group QQ^n * n^QQ).
Change History (36)
comment:1 Changed 5 years ago by
comment:2 in reply to: ↑ description Changed 5 years ago by
Replying to cheuberg:
sage: A.<n> = AsymptoticRing('QQ^n * n^QQ', ZZ) sage: (1/2)^n
The problem is that (1/2)^n
calls __pow__
of the rational 1/2
. This tries A(1/2)^n
, which fails since 1/2
is not in A
. This also explains why it works with coefficient ring QQ
as in comment:1.
As a workaround you can use
sage: n.rpow(1/2)
comment:3 Changed 5 years ago by
Documented at #19961.
comment:4 Changed 3 years ago by
 Component changed from asymptotic expansions to basic arithmetic
 Dependencies set to #24247
 Summary changed from Asymptotic Ring: cannot construct (1/2)^n to Fix Rational.__pow__
comment:5 Changed 3 years ago by
 Description modified (diff)
comment:6 Changed 3 years ago by
 Branch set to u/jdemeyer/ticket/19945
comment:7 Changed 3 years ago by
 Commit set to 1c9bb4f78f10157d7b02a385bbfbb1f3d313e9ae
Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
1c9bb4f  Implement Rational.__pow__ using the coercion model

comment:8 Changed 3 years ago by
 Status changed from new to needs_review
comment:9 Changed 3 years ago by
 Commit changed from 1c9bb4f78f10157d7b02a385bbfbb1f3d313e9ae to 05b066dda9ce103a4ea48a6c97596d29f085ad2e
Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
05b066d  Implement Rational.__pow__ using the coercion model

comment:10 Changed 3 years ago by
 Commit changed from 05b066dda9ce103a4ea48a6c97596d29f085ad2e to 84c056bc5aefc9c7525cc2d1c7a5758e0f33dfe9
Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
84c056b  Implement Rational.__pow__ using the coercion model

comment:11 Changed 3 years ago by
 Status changed from needs_review to needs_work
sage t long warnlong 63.9 src/sage/rings/rational.pyx ********************************************************************** File "src/sage/rings/rational.pyx", line 302, in sage.rings.rational.rational_power_parts Failed example: (1)^(1/3) Expected: (1)^(2/3) Got: 1.0*(1)^(2/3) **********************************************************************
comment:12 followup: ↓ 13 Changed 3 years ago by
Seriously? How did that happen?
comment:13 in reply to: ↑ 12 Changed 3 years ago by
Replying to jdemeyer:
Seriously? How did that happen?
No idea. I just thought I would review this ticket, merged your branch in the current develop
, and ran the test suite.
comment:14 Changed 2 years ago by
 Dependencies #24247 deleted
comment:15 Changed 2 years ago by
 Commit changed from 84c056bc5aefc9c7525cc2d1c7a5758e0f33dfe9 to b74fb2005c709b90eb2821de3b20615bea6ed936
Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
b74fb20  Implement Rational.__pow__ using the coercion model

comment:16 Changed 2 years ago by
 Commit changed from b74fb2005c709b90eb2821de3b20615bea6ed936 to b862b7eedff3d9c75ea821ce4bdec379e06a82c7
Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
b862b7e  Implement Rational.__pow__ using the coercion model

comment:17 Changed 2 years ago by
 Milestone changed from sage7.1 to sage8.5
 Status changed from needs_work to needs_review
Rebased to 8.5.beta5.
comment:18 Changed 2 years ago by
 Reviewers set to Travis Scrimshaw
In general, this LGTM. Have you run any timings to see if there is any change for Q^{Z}, Q^{Q}, and Z^{Q}? Since these are common operations and they now go through the coercion framework (if I understand everything correctly), I want to make sure they are not (significantly) affected.
comment:19 Changed 2 years ago by
I did some timings with mixed results:
QQ^ZZ
became slower by a factor 1.4 (404ns to 575ns). This is because of the overhead ofIntegerPowAction
and it's really hard to avoid this slowdown. It's also a large relative slowdown because the actual powering is very fast.
 The old code optimized the case of powering to a rational integer (an element of
QQ
with denominator 1). The new code doesn't, leading to a slowdown forQQ^QQ
with integral exponent.
QQ^QQ
is faster if the exponent is not a rational integer.
ZZ^QQ
simply delegates toQQ^QQ
(this code is not changed), so speed gains/losses are similar to theQQ^QQ
case.
I'm not entirely sure whether I should implement a special case for powering to a rational integer. Is this case common enough to warrant the extra code?
I also noticed that #24500 gives an additional slight slowdown for QQ^ZZ
(about 5%). This is because powering is always a right action, so the problem that #24500 fixes does not occur for IntegerPowAction
.
comment:20 Changed 2 years ago by
#26775 should solve most of the slowdown for QQ^ZZ
.
comment:21 Changed 2 years ago by
Never mind, my analysis at #26775 was wrong.
Still, there is something strange with benchmarks: without changing any code, the old code slowed down from 404ns to 446ns and the new code (on top of #24500) sped up from about 600ns to 515ns. So without really changing anything (but recompiling the new code), I now have a slowdown of a factor only 1.15 for QQ^ZZ
. Benchmarks are a strange thing...
comment:22 Changed 2 years ago by
The wildly differing timings must be something like a memory layout/cache issue.
Anyway, for QQ^ZZ
, I give the best timings for
sage: a = QQ(25); b = ZZ(2); timeit('a**b', repeat=20, number=100000)
Vanilla 8.5.beta5:
100000 loops, best of 20: 397 ns per loop
With #19945:
100000 loops, best of 20: 512 ns per loop
100000 loops, best of 20: 509 ns per loop
comment:23 Changed 2 years ago by
Most of the slowdown in this ticket is due to
sage: g = coercion_model.get_action; args = (QQ, ZZ, operator.pow) sage: timeit('g(*args)', repeat=20, number=100000) 100000 loops, best of 20: 133 ns per loop
comment:24 followup: ↓ 26 Changed 2 years ago by
That is what I expected. I was wondering if we wanted to special case such things to avoid going through the coercion framework. I also think it is worthwhile to special case rational powers that are integers.
comment:25 Changed 2 years ago by
#26776 should gain about 20ns.
comment:26 in reply to: ↑ 24 ; followup: ↓ 29 Changed 2 years ago by
Replying to tscrim:
I was wondering if we wanted to special case such things to avoid going through the coercion framework.
More useful (if you want to do this optimization) would be to specialcase integer powering in general in the coercion model.
comment:27 Changed 2 years ago by
sage: a = QQ(25); b = ZZ(2); timeit('a**b', repeat=20, number=100000) 100000 loops, best of 20: 493 ns per loop
comment:28 Changed 2 years ago by
With the improvements from #26776, the slowdown is "only" 20%.
comment:29 in reply to: ↑ 26 ; followup: ↓ 30 Changed 2 years ago by
Replying to jdemeyer:
Replying to tscrim:
I was wondering if we wanted to special case such things to avoid going through the coercion framework.
More useful (if you want to do this optimization) would be to specialcase integer powering in general in the coercion model.
I thought we already had some special casing for this with IntegerPowAction
?
IMO, it is really easy to accidentally construct integers that are rational numbers, so that is why I think there is some merit to special casing that. However, I don't know how often this occurs in practice. My gut tells me that exponentiation of rationals is usually not a bottleneck in most computations, so I might be being too precious here.
comment:30 in reply to: ↑ 29 Changed 2 years ago by
Replying to tscrim:
I thought we already had some special casing for this with
IntegerPowAction
?
No, we don't. It's just an ordinary action.
IMO, it is really easy to accidentally construct integers that are rational numbers, so that is why I think there is some merit to special casing that.
Sure, I'll do that.
comment:31 Changed 2 years ago by
 Status changed from needs_review to needs_work
comment:32 Changed 2 years ago by
 Commit changed from b862b7eedff3d9c75ea821ce4bdec379e06a82c7 to 5706029d8199f2e7d0e66ba84f19a44f3336e04c
Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
5706029  Implement Rational.__pow__ using the coercion model

comment:33 Changed 2 years ago by
 Status changed from needs_work to needs_review
I optimized powering by rational integers, which is now much faster than before.
comment:34 Changed 2 years ago by
Thank you and also for doing the initial timings as well. When the patchbot comes back green, you can set this to a positive review.
comment:35 Changed 2 years ago by
 Status changed from needs_review to positive_review
comment:36 Changed 2 years ago by
 Branch changed from u/jdemeyer/ticket/19945 to 5706029d8199f2e7d0e66ba84f19a44f3336e04c
 Resolution set to fixed
 Status changed from positive_review to closed
does work.