Opened 4 years ago

Closed 12 months ago

#19945 closed defect (fixed)

Fix Rational.__pow__

Reported by: cheuberg Owned by:
Priority: major Milestone: sage-8.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 jdemeyer)

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).

Change History (36)

comment:1 Changed 4 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 4 years ago by dkrenn

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 4 years ago by dkrenn

Documented at #19961.

comment:4 Changed 2 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 2 years ago by jdemeyer

  • Description modified (diff)

comment:6 Changed 2 years ago by jdemeyer

  • Branch set to u/jdemeyer/ticket/19945

comment:7 Changed 2 years ago by git

  • Commit set to 1c9bb4f78f10157d7b02a385bbfbb1f3d313e9ae

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

1c9bb4fImplement Rational.__pow__ using the coercion model

comment:8 Changed 2 years ago by jdemeyer

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

comment:9 Changed 2 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:

05b066dImplement Rational.__pow__ using the coercion model

comment:10 Changed 2 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:

84c056bImplement Rational.__pow__ using the coercion model

comment:11 Changed 22 months 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: Changed 22 months ago by jdemeyer

Seriously? How did that happen?

comment:13 in reply to: ↑ 12 Changed 22 months ago by mmezzarobba

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

  • Dependencies #24247 deleted

comment:15 Changed 13 months 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:

b74fb20Implement Rational.__pow__ using the coercion model

comment:16 Changed 13 months 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:

b862b7eImplement Rational.__pow__ using the coercion model

comment:17 Changed 13 months 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 13 months 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 13 months 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 13 months ago by jdemeyer

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

comment:21 Changed 13 months 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 13 months 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 13 months ago by jdemeyer (previous) (diff)

comment:23 Changed 13 months 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 13 months ago by jdemeyer (previous) (diff)

comment:24 follow-up: Changed 13 months 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:25 Changed 13 months ago by jdemeyer

#26776 should gain about 20ns.

comment:26 in reply to: ↑ 24 ; follow-up: Changed 13 months ago by 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 special-case integer powering in general in the coercion model.

comment:27 Changed 13 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 13 months ago by jdemeyer (previous) (diff)

comment:28 Changed 13 months ago by jdemeyer

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

comment:29 in reply to: ↑ 26 ; follow-up: Changed 13 months ago by tscrim

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

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

  • Status changed from needs_review to needs_work

comment:32 Changed 13 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:

5706029Implement Rational.__pow__ using the coercion model

comment:33 Changed 13 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 13 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 13 months ago by jdemeyer

  • Status changed from needs_review to positive_review

comment:36 Changed 12 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.