Opened 5 years ago

Closed 23 months ago

# Fix Rational.__pow__

Reported by: Owned by: cheuberg major sage-8.5 basic arithmetic dkrenn, behackl Jeroen Demeyer Travis Scrimshaw N/A 5706029 (Commits) 5706029d8199f2e7d0e66ba84f19a44f3336e04c

The following is unexpected:

```sage: A.<n> = AsymptoticRing('QQ^n * n^QQ', ZZ); (1/2)^n
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-18-784db7f7676d> 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/sage-config/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/sage-config/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/sage-config/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/sage-config/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/sage-config/local/lib/python2.7/site-packages/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/sage-config/local/lib/python2.7/site-packages/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/sage-config/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/sage-config/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/sage-config/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/sage-config/local/lib/python2.7/site-packages/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).
```

### comment:1 Changed 5 years ago by cheuberg

```sage: A.<n> = AsymptoticRing('QQ^n*n^QQ', QQ)
sage: (1/2)^n
(1/2)^n
```

does work.

### comment:2 in reply to: ↑ description Changed 5 years ago by dkrenn

```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 dkrenn

Documented at #19961.

### comment:4 Changed 3 years ago by jdemeyer

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

• Description modified (diff)

### comment:6 Changed 3 years ago by jdemeyer

• Branch set to u/jdemeyer/ticket/19945

### comment:7 Changed 3 years ago by git

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

• Authors set to Jeroen Demeyer
• Status changed from new to needs_review

### comment:9 Changed 3 years ago by git

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

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

• Status changed from needs_review to needs_work
```sage -t --long --warn-long 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 follow-up: ↓ 13 Changed 3 years ago by jdemeyer

Seriously? How did that happen?

### comment:13 in reply to: ↑ 12 Changed 3 years ago by mmezzarobba

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 jdemeyer

• Dependencies #24247 deleted

### comment:15 Changed 2 years ago by git

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

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

• Milestone changed from sage-7.1 to sage-8.5
• Status changed from needs_work to needs_review

Rebased to 8.5.beta5.

### comment:18 Changed 2 years ago by tscrim

• Reviewers set to Travis Scrimshaw

In general, this LGTM. Have you run any timings to see if there is any change for QZ, QQ, and ZQ? 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 jdemeyer

I did some timings with mixed results:

1. `QQ^ZZ` became slower by a factor 1.4 (404ns to 575ns). This is because of the overhead of `IntegerPowAction` and it's really hard to avoid this slow-down. It's also a large relative slow-down because the actual powering is very fast.
1. 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 slow-down for `QQ^QQ` with integral exponent.
1. `QQ^QQ` is faster if the exponent is not a rational integer.
1. `ZZ^QQ` simply delegates to `QQ^QQ` (this code is not changed), so speed gains/losses are similar to the `QQ^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 slow-down 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 jdemeyer

#26775 should solve most of the slow-down for `QQ^ZZ`.

### comment:21 Changed 2 years ago by jdemeyer

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 slow-down of a factor only 1.15 for `QQ^ZZ`. Benchmarks are a strange thing...

### comment:22 Changed 2 years ago by jdemeyer

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

With #24500 and #19945:

```100000 loops, best of 20: 509 ns per loop
```
Last edited 2 years ago by jdemeyer (previous) (diff)

### comment:23 Changed 2 years ago by jdemeyer

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
```
Last edited 2 years ago by jdemeyer (previous) (diff)

### comment:24 follow-up: ↓ 26 Changed 2 years ago by tscrim

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:26 in reply to: ↑ 24 ; follow-up: ↓ 29 Changed 23 months ago by jdemeyer

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 special-case integer powering in general in the coercion model.

### comment:27 Changed 23 months ago by jdemeyer

With #26776 + #19945:

```sage: a = QQ(25); b = ZZ(2); timeit('a**b', repeat=20, number=100000)
100000 loops, best of 20: 479 ns per loop
```
Last edited 23 months ago by jdemeyer (previous) (diff)

### comment:28 Changed 23 months ago by jdemeyer

With the improvements from #26776, the slow-down is "only" 20%.

### comment:29 in reply to: ↑ 26 ; follow-up: ↓ 30 Changed 23 months ago by 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 special-case 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 23 months ago by jdemeyer

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

• Status changed from needs_review to needs_work

### comment:32 Changed 23 months ago by git

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

• Status changed from needs_work to needs_review

I optimized powering by rational integers, which is now much faster than before.

### comment:34 Changed 23 months ago by tscrim

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

• Status changed from needs_review to positive_review

### comment:36 Changed 23 months ago by vbraun

• Branch changed from u/jdemeyer/ticket/19945 to 5706029d8199f2e7d0e66ba84f19a44f3336e04c
• Resolution set to fixed
• Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.