Opened 3 years ago
Closed 3 years ago
#20086 closed defect (fixed)
rational powers in ZZ[X] and QQ[X]
Reported by:  cheuberg  Owned by:  

Priority:  major  Milestone:  sage7.2 
Component:  basic arithmetic  Keywords:  
Cc:  behackl, dkrenn  Merged in:  
Authors:  Clemens Heuberger, Vincent Delecroix, Benjamin Hackl  Reviewers:  Benjamin Hackl, Vincent Delecroix 
Report Upstream:  N/A  Work issues:  
Branch:  5fc1027 (Commits)  Commit:  5fc1027e43d601a791568b22f12ec70957c47519 
Dependencies:  Stopgaps: 
Description (last modified by )
Until now,
sage: R.<x> = ZZ[] sage: R(1)^(1/2) Traceback (most recent call last): ... TypeError: rational is not an integer
because only integer exponents were allowed for polynomials.
Implement arbitrary powers of constant polynomials by handing over to the rational field.
This was originally observed in the asymptotic ring:
sage: P.<R> = QQ[] sage: A.<Z> = AsymptoticRing('T^QQ', P) sage: sqrt(Z) Traceback (most recent call last): ... ArithmeticError: Cannot take T to the exponent 1/2 in Exact Term Monoid T^QQ with coefficients in Univariate Polynomial Ring in R over Rational Field since its coefficient 1 cannot be taken to this exponent. > *previous* TypeError: rational is not an integer
follow up: #20571
Change History (94)
comment:1 Changed 3 years ago by
 Branch set to u/cheuberg/polynomials/power
comment:2 Changed 3 years ago by
 Commit set to cc49098e1fc111e63d344b6c2f914d0d3e5e463b
comment:3 Changed 3 years ago by
 Component changed from asymptotic expansions to basic arithmetic
 Description modified (diff)
 Status changed from new to needs_review
 Summary changed from Asymptotic ring: fix exponentiation to QQ[X]: allow arbitrary powers of constant polynomials
comment:4 Changed 3 years ago by
 Commit changed from cc49098e1fc111e63d344b6c2f914d0d3e5e463b to 25fdd0c0ceb085d68b2ddf851983a386771d164f
Branch pushed to git repo; I updated commit sha1. New commits:
25fdd0c  Trac #20086: powers of constant polynomials over ZZ

comment:5 Changed 3 years ago by
 Summary changed from QQ[X]: allow arbitrary powers of constant polynomials to ZZ[X], QQ[X]: allow arbitrary powers of constant polynomials
comment:6 followup: ↓ 9 Changed 3 years ago by
The fact that for a in QQ, the value of a^(1/2)
has its parent depending on the actual value of a is a compromise for novice use (in calculus etc.) Normally, one would raise an error if a^(1/2)
does not lie in QQ. Certainly, promoting to SR is a very bad choice for anything else than novice use.
I see you take the effort of trying to put the result back in the original parent. Don't you think it's better to force that? i.e., raise an error if it doesn't work, rather than give a result back in SR?
The problem is that if someone is computing with polynomials, it's almost certainly not desired to end up in SR (where things like QQ['x']['x']
get squashed), and if an error isn't raised, it's very hard to detect that it happened.
Also, if you're implementing the (partial) map a:>a^(n/m)
, why not extend it properly? When
(QQ['x'](4))^(1/2)
works, why should QQ['x'](x^2)^(1/2)
fail?
comment:7 Changed 3 years ago by
 Commit changed from 25fdd0c0ceb085d68b2ddf851983a386771d164f to 6217beeb2520f5e8c0d1e92881c4d1c9569e9868
Branch pushed to git repo; I updated commit sha1. New commits:
6217bee  Trac #20085: raise exception instead of moving to SR

comment:8 Changed 3 years ago by
 Commit changed from 6217beeb2520f5e8c0d1e92881c4d1c9569e9868 to 918499cf377a551a5e28ea1d6748ed50860f1afa
Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
918499c  Trac #20086: raise exception instead of moving to SR

comment:9 in reply to: ↑ 6 Changed 3 years ago by
Replying to nbruin:
The fact that for a in QQ, the value of
a^(1/2)
has its parent depending on the actual value of a is a compromise for novice use (in calculus etc.) Normally, one would raise an error ifa^(1/2)
does not lie in QQ. Certainly, promoting to SR is a very bad choice for anything else than novice use.I see you take the effort of trying to put the result back in the original parent. Don't you think it's better to force that? i.e., raise an error if it doesn't work, rather than give a result back in SR?
The problem is that if someone is computing with polynomials, it's almost certainly not desired to end up in SR (where things like
QQ['x']['x']
get squashed), and if an error isn't raised, it's very hard to detect that it happened.
very valid points, thank you. I now raise an exception.
Also, if you're implementing the (partial) map
a:>a^(n/m)
, why not extend it properly? When(QQ['x'](4))^(1/2)
works, why shouldQQ['x'](x^2)^(1/2)
fail?
I needed a fix a bug which did bite me somewhere else, see initial description. Of course, it would be nice to have that, too; but I'd prefer to leave that to a followup ticket and to have this functionality here soon.
comment:10 Changed 3 years ago by
 Cc dkrenn added; dakrenn removed
comment:11 followup: ↓ 15 Changed 3 years ago by
 Branch changed from u/cheuberg/polynomials/power to u/behackl/polynomials/power
 Commit changed from 918499cf377a551a5e28ea1d6748ed50860f1afa to e6e1c812ecb49be56dc4d9b4ed2e524011c391c9
I've fixed the segmentation fault from rings/integer.pyx
and implemented the case of constant^constant
.
Also, +1 for extending this functionality on a followup ticket; I'd also like to have this in as soon as possible.
I've reviewed your changes, please crossreview. If you are satisfied and if there are no other objections, I'd set this to positive_review
.
New commits:
caa02c8  Merge tag '7.1.beta6' into HEAD

3f44399  fix segmentation fault

5934ed1  fix (constant poly)^(constant poly)

e6e1c81  add doctests

comment:12 Changed 3 years ago by
 Reviewers set to Benjamin Hackl
comment:13 Changed 3 years ago by
 Commit changed from e6e1c812ecb49be56dc4d9b4ed2e524011c391c9 to 87a22b6fc4cff82491a7b0ba35983bc5dc1299df
Branch pushed to git repo; I updated commit sha1. New commits:
87a22b6  fix "referenced before assignment"

comment:14 Changed 3 years ago by
 Branch changed from u/behackl/polynomials/power to u/cheuberg/polynomials/power
comment:15 in reply to: ↑ 11 Changed 3 years ago by
 Commit changed from 87a22b6fc4cff82491a7b0ba35983bc5dc1299df to 14f7efaaf5be81405d1954473b4dc3f47100e486
Replying to behackl:
I've fixed the segmentation fault from
rings/integer.pyx
and implemented the case ofconstant^constant
.
That's weird, but as integer tries to convert the base to the base of the exponent and expects the polynomial ring to handle it, we have to cope with it.
I've reviewed your changes, please crossreview. If you are satisfied and if there are no other objections, I'd set this to
positive_review
.
I refactored the code: that way it seems to be more suitable for future extension and handles polynomial exponents gracefully.
In a first commit, I removed handling of the segmentation fault because I thought it should be handeled by the base ring; later on, I understood that it has to be handeled here, so reinstated fix (in some other way).
New commits:
e95c9b9  Trac #20086: refactor method; do not fix segmentation fault 2^R in polynomial ring

eadecba  Trac #20086: fix ReSt errors (docbuild u k complains)

a139868  Trac #20086: handle nonconstant polynomial exponents

14f7efa  Trac #20086: ReSt formatting

comment:16 followup: ↓ 17 Changed 3 years ago by
 Status changed from needs_review to needs_info
Why doing something that complicated for dealing with strange exponents? Why not simply
if input is not an integer: convert to an integer do the exponentiation
Why are you special casing degree 0 polynomial? If I understand correctly these examples won't behave similarly
sage: R.<x> = QQ[] sage: ((x+1)^2)^(1/2) > TypeError sage: ((R(2))^2)^(1/2) > 2
Which is weird.
comment:17 in reply to: ↑ 16 ; followup: ↓ 18 Changed 3 years ago by
 Status changed from needs_info to needs_review
Replying to vdelecroix:
Why doing something that complicated for dealing with strange exponents? Why not simply
if input is not an integer: convert to an integer do the exponentiation
I do not understand this comment. Which input? Base or Exponent? Where should this code be?
Why are you special casing degree 0 polynomial? If I understand correctly these examples won't behave similarly
sage: R.<x> = QQ[] sage: ((x+1)^2)^(1/2) > TypeError sage: ((R(2))^2)^(1/2) > 2Which is weird.
I need a quick fix for R(1)^(1/2)
. If somebody has time to implement ((x+1)^2)^(1/2)
very soon, I'd be happy. I do not have time soon. However, I want to have the code associated with a recently submitted paper in 7.1.
Basically, this fix here simply branches to existing code. Computing ((x+1)^2)^(1/2)
needs new mathematical code (involving square free decomposition).
Therefore, I propose to include this now on the basis that while this is not a perfect and definitive solution, it is better than the previous behaviour.
comment:18 in reply to: ↑ 17 ; followups: ↓ 20 ↓ 24 Changed 3 years ago by
Replying to cheuberg:
Replying to vdelecroix:
Why doing something that complicated for dealing with strange exponents? Why not simply
if input is not an integer: convert to an integer do the exponentiationI do not understand this comment. Which input? Base or Exponent? Where should this code be?
Let me be more precise then
def __pow__(self, exp): cdef long n try: n = exp # old code for integer exponent except TypeError: n = QQ.coerce(exp) # new code for rational exponent
What I do not understand is why are you testing if the exponent is a polynomial...
The following is currently a TypeError
sage: 1^(ZZ['x'].one()) Traceback (most recent call last): ... TypeError: 'sage.rings.polynomial.polynomial_integer_dense_flint.Polynomial_integer_dense_flint' object cannot be interpreted as an index
The same should happen when the integer 1
is replaced by the polynomial 1
.
Why are you special casing degree 0 polynomial? If I understand correctly these examples won't behave similarly
sage: R.<x> = QQ[] sage: ((x+1)^2)^(1/2) > TypeError sage: ((R(2))^2)^(1/2) > 2Which is weird.
I need a quick fix for
R(1)^(1/2)
. If somebody has time to implement((x+1)^2)^(1/2)
very soon, I'd be happy. I do not have time soon. However, I want to have the code associated with a recently submitted paper in 7.1.
Are you sure there was a bug? In Sage the integer 1 is *not* the 1 from ZZ[x] (though they are equal through coercion). Some softwares behave differently to that respect (e.g. GAP) where there is only one 1 which is an integer and not anything else. In Sage (but not in GAP) it is already the case that operations change with respect to the parents even if the objects are equal
sage: Zmod(10)(3) == 3 sage: Zmod(10)(5) == 5 sage: log(Zmod(10)(3)) 3 sage: log(Zmod(10)(5)) Traceback (most recent call last): ... ZeroDivisionError: Inverse does not exist.
Basically, this fix here simply branches to existing code. Computing
((x+1)^2)^(1/2)
needs new mathematical code (involving square free decomposition).
Indeed. Your modifications obfuscate the code for only a very trivial case and a discutable behavior of powering with polynomials.
Therefore, I propose to include this now on the basis that while this is not a perfect and definitive solution, it is better than the previous behaviour.
Not sure it is better. Sage used to consider operations based on parents. Powers are of course a special type of operation and with some respect might be treated appart. But "(polynomial)^{(polynomial)" is not well defined. And "(polynomial)}(rational)" is well defined in some situations (and to that respect, your code improves the current situation a little).
Moreover, the current "power promotion" for ZZ is very bad
sage: 3^(1/2) sqrt(3) sage: parent(_) Symbolic Ring sage: 4^(1/2) 2 sage: parent(_) Rational Field
Making it available for constant polynomial is not that much of an improvement.
comment:19 Changed 3 years ago by
 Description modified (diff)
 Status changed from needs_review to needs_work
It seems that I misedited the description of the problem at some stage. Is now fixed.
Apart from that, the current code apparently leads to problems.
comment:20 in reply to: ↑ 18 Changed 3 years ago by
Replying to vdelecroix:
Moreover, the current "power promotion" for ZZ is very bad
sage: 3^(1/2) sqrt(3) sage: parent(_) Symbolic Ring sage: 4^(1/2) 2 sage: parent(_) Rational Field
see nbruin's comment 6.
comment:21 Changed 3 years ago by
Hi Vincent!
Just to clarify the previous need for checking whether the exponent is a polynomial: in our previous approach, we were hoping for the coefficient ring to carry out the exponentiation, i.e. if we had
sage: P.<z> = QQ[] sage: P(1/4)^(1/2) 1/2
then we would compute this by letting QQ
do the exponentiation. However, in the doctest over at rings/integer.pyx
we have something like
sage: P.<t> = QQ[] sage: 2^t
With the old implementation, this lead to an infinite loop (because coercion would put the base always back to P again...), and thus I added the check for the polynomial in the exponent.
In any case, I'm all for improving the code and first trying to deal with integers, and then with rationals afterwards. Nevertheless, I'd still just handle constant polynomials in the base on this ticketbut then, this could be extended easily.
Benjamin
comment:22 Changed 3 years ago by
 Commit changed from 14f7efaaf5be81405d1954473b4dc3f47100e486 to 892109ac4f48a4a32d35671365ecceb0d1c9e0bf
Branch pushed to git repo; I updated commit sha1. New commits:
892109a  Trac #20086: fix bug in previous commit (and use better doctest)

comment:23 Changed 3 years ago by
 Status changed from needs_work to needs_review
Problem found by patchbot was a stupid mistake, hopefully fixed now.
comment:24 in reply to: ↑ 18 ; followup: ↓ 25 Changed 3 years ago by
Replying to vdelecroix:
Let me be more precise then
def __pow__(self, exp): cdef long n try: n = exp # old code for integer exponent except TypeError: n = QQ.coerce(exp) # new code for rational exponent
ok, if you prefer it like that, we can probably do that.
I need a quick fix for
R(1)^(1/2)
. If somebody has time to implement((x+1)^2)^(1/2)
very soon, I'd be happy. I do not have time soon. However, I want to have the code associated with a recently submitted paper in 7.1.Are you sure there was a bug?
Do you think that disallowing R(1)^(1/2)
is the desired behaviour?
In Sage the integer 1 is *not* the 1 from ZZ[x] (though they are equal through coercion). Some softwares behave differently to that respect (e.g. GAP) where there is only one 1 which is an integer and not anything else. In Sage (but not in GAP) it is already the case that operations change with respect to the parents even if the objects are equal
sage: Zmod(10)(3) == 3 sage: Zmod(10)(5) == 5 sage: log(Zmod(10)(3)) 3 sage: log(Zmod(10)(5)) Traceback (most recent call last): ... ZeroDivisionError: Inverse does not exist.
I have no idea how this is related to this problem, sorry.
Basically, this fix here simply branches to existing code. Computing
((x+1)^2)^(1/2)
needs new mathematical code (involving square free decomposition).Indeed. Your modifications obfuscate the code for only a very trivial case and a discutable behavior of powering with polynomials.
It might seem trivial to you, but it did cost me an hour while writing a paper because basically, asymptotic rings using polynomial rings as coefficient rings could not compute square roots, and, sorry, I need that.
Therefore, I propose to include this now on the basis that while this is not a perfect and definitive solution, it is better than the previous behaviour.
Not sure it is better. Sage used to consider operations based on parents. Powers are of course a special type of operation and with some respect might be treated appart. But "(polynomial)^{(polynomial)" is not well defined. And "(polynomial)}(rational)" is well defined in some situations (and to that respect, your code improves the current situation a little).
There was a discussion on sagedevel a few weeks ago. Every parent seems to have its own homemade logic about how to do coercion. Please do not try to fix that in this ticket.
comment:25 in reply to: ↑ 24 ; followup: ↓ 26 Changed 3 years ago by
Replying to cheuberg:
[... snip ...]
There was a discussion on sagedevel a few weeks ago. Every parent seems to have its own homemade logic about how to do coercion. Please do not try to fix that in this ticket.
Please do not try to complicate it in this ticket either. A possible solution could be:
 allow
polynomial^rational
as much as we can (i.e. for some constant polynomials for the sake of this ticket). For more general polynomial we can for now raise aNotImplementedError
.  disallow
polynomial^polynomial
even if the power is a constant polynomial. In that case, the error should be aTypeError
. This is the current behavior for some polynomial rings but not allsage: ZZ['x'].gen() ** ZZ['x'].one() > TypeError sage: QQ['x'].gen() ** ZZ['x'].one() > TypeError sage: RR['x'].gen() ** ZZ['x'].one() x
What do you think? Would that solution fit your needs?
Vincent
comment:26 in reply to: ↑ 25 Changed 3 years ago by
Replying to vdelecroix:
A possible solution could be:
 allow
polynomial^rational
as much as we can (i.e. for some constant polynomials for the sake of this ticket). For more general polynomial we can for now raise aNotImplementedError
. disallow
polynomial^polynomial
even if the power is a constant polynomial. In that case, the error should be aTypeError
. This is the current behavior for some polynomial rings but not allsage: ZZ['x'].gen() ** ZZ['x'].one() > TypeError sage: QQ['x'].gen() ** ZZ['x'].one() > TypeError sage: RR['x'].gen() ** ZZ['x'].one() xWhat do you think? Would that solution fit your needs?
As far as I can see without implementing it and make ptestlong
, this fits perfectly.
comment:27 Changed 3 years ago by
 Commit changed from 892109ac4f48a4a32d35671365ecceb0d1c9e0bf to 1782a3cc92f875fa3ac8a83395d204d2ca3d7f8c
Branch pushed to git repo; I updated commit sha1. New commits:
1782a3c  Trac #20086: rewrite try/expect; disallow polynomial exponents

comment:28 followup: ↓ 29 Changed 3 years ago by
This is now refactored as a try/expect block; polynomial exponents are now disallowed.
I am not completely comfortable with the long block wrapped within try
/expect
because I am not sure that we catch the correct TypeError
. What do you think about that?
comment:29 in reply to: ↑ 28 ; followup: ↓ 31 Changed 3 years ago by
Replying to cheuberg:
This is now refactored as a try/expect block; polynomial exponents are now disallowed.
Thanks.
I am not completely comfortable with the long block wrapped within
try
/expect
because I am not sure that we catch the correctTypeError
. What do you think about that?
Right. It is better to avoid long try/except
. One possibility is
try: nn = pyobject_to_long(exp) except TypeError: # rational code else: # integer code
comment:30 Changed 3 years ago by
 Commit changed from 1782a3cc92f875fa3ac8a83395d204d2ca3d7f8c to 9bda3669c14463b3061f19c959bf8a22cb467117
Branch pushed to git repo; I updated commit sha1. New commits:
9bda366  Trac #20086: refactor as try/except/else

comment:31 in reply to: ↑ 29 Changed 3 years ago by
Replying to vdelecroix:
Right. It is better to avoid long
try/except
. One possibility istry: nn = pyobject_to_long(exp) except TypeError: # rational code else: # integer code
done.
comment:32 Changed 3 years ago by
 Branch changed from u/cheuberg/polynomials/power to u/vdelecroix/20086
 Commit changed from 9bda3669c14463b3061f19c959bf8a22cb467117 to cba01eb6ebe741567cabc1e9e4a3dd0c94aa11b2
 Summary changed from ZZ[X], QQ[X]: allow arbitrary powers of constant polynomials to rational powers in ZZ[X] and QQ[X]
I implemented a generic nth_root
method and used it...
New commits:
cba01eb  Trac 20086: implement (polynomial)^(rational)

comment:33 Changed 3 years ago by
Though, since num
and den
are relatively prime it should be cheaper to do self.nth_root(den) ** num
instead of (self ** num).nth_root(den)
... I will think about it.
Anyway, it would be good to see whether it works for other polynomial rings. For example QQbar['x']
or when the base ring is a number field or a finite field...
comment:34 Changed 3 years ago by
When the polynomial is of degree 0 it would also be faster to actually also relies on nth_root
:
sage: %timeit 27.nth_root(3) 1000000 loops, best of 3: 388 ns per loop sage: %timeit 27 ** (1/3) 100000 loops, best of 3: 7.51 µs per loop
As nth_root
should be guaranteed to return an element of the same ring there will also be much less coercion involved.
comment:35 Changed 3 years ago by
 Commit changed from cba01eb6ebe741567cabc1e9e4a3dd0c94aa11b2 to c603c2b293819088246b34ab676f303871d5cb40
Branch pushed to git repo; I updated commit sha1. New commits:
c603c2b  Trac 20086: improvements / simplification

comment:36 Changed 3 years ago by
 Commit changed from c603c2b293819088246b34ab676f303871d5cb40 to 42c53199edeb619d8b4867cea0b318d7d0c8ea0a
Branch pushed to git repo; I updated commit sha1. New commits:
42c5319  Trac 20086: fix doctest in polynomial_quotient_ring.py

comment:37 Changed 3 years ago by
 Branch changed from u/vdelecroix/20086 to u/behackl/polynomial/rationalpowers
 Commit changed from 42c53199edeb619d8b4867cea0b318d7d0c8ea0a to c8350c631287dd71a9b79fb422d769a1e480cdc8
Hi! I've reviewed this refactored implementation (thanks for implementing it!), and everything looks good to me.
I've added a small reviewer commit; please crossreview and set to positive_review
if you are satisfied.
New commits:
4313d3f  Merge tag '7.2.beta1' into polynomial/rationalpowers

c8350c6  improve exception messages

comment:38 Changed 3 years ago by
Nope. This is actually working because of a bug see #20214 ;(. What I wrote in the generic nth_root
should always be an infinite loop (because of u.nth_root(n)
).
One possibility:
 make it work when
u.is_one()
isTrue
(or more generally whenu.multiplicative_order()
is finite)  raise an error if
u
is not in the above case and there are some non trivial factors
The aim of the second item is that in children classes (like polynomials) you might want to do
def nth_root(self, n): if self.degree() <= 0: # factorize unit using base ring method XXX else: return super(XXX).nth_root(self, n)
comment:39 Changed 3 years ago by
 Status changed from needs_review to needs_work
comment:40 Changed 3 years ago by
Well, wouldn't it be more natural to let the unit be an element of the base ring? (As far as I see the parent of the unit of polynomials over QQ always is from QQ, for example.)
Of course, we can also add checks like if u
is one etc., but with this ticket only polynomial rings over the rationals and the integers use this methodand both of them have implemented a nth_root
method; this is why I still think that this would be good to go and why special treatment isn't needed. And the inconsistency mentioned in #20214 only introduces a recursion of depth 2, which would be resolved if the unit
method would behave for integers like for rationals.
It isn't even necessary to separately implement nth_root
for polynomials such that the case of constant polynomials is handled by the base ring: over ZZ
, the overall coefficient is decomposed w.r.t. PFD and handled like a nonconstant polynomial (which is what I would have implemented in the base ring as well). For QQ, the overall factor is in the unit and handled separately in the nth_root
method.
Letting other polynomial rings profit from this procedure should be realized in a followup ticket, IMHO.
However, what should be added before this ships is a special treatment of the zero polynomial. I'll push this in a minute.
comment:41 Changed 3 years ago by
 Commit changed from c8350c631287dd71a9b79fb422d769a1e480cdc8 to 5bdbd3232fb25bda9de40d7a8189036eb80b0c19
Branch pushed to git repo; I updated commit sha1. New commits:
5bdbd32  special treatment for zero polynomial

comment:42 followup: ↓ 43 Changed 3 years ago by
... any new thoughts on this?
comment:43 in reply to: ↑ 42 Changed 3 years ago by
 Status changed from needs_work to needs_info
Replying to behackl:
... any new thoughts on this?
I would very much like to have this in 7.2
, because I/we need it for a paper we have written. Thus I'd like to discuss this once again:
IIRC the possibly infinite recursion could be avoided altogether if we would convert the unit f.unit()
to the corresponding base ring. This wouldn't require any additional fix for the units, and the code affects only ZZ[]
and QQ[]
anyway, so I'm rather sure that it works.
Any objections to this approach?
comment:44 followup: ↓ 46 Changed 3 years ago by
It is fine for me if you modify the generic nth_root
to
u = f.unit() if u.is_one(): # do nothing elif self.parent() != self.base_ring(): # try to factorize the unit in the base ring else: # raise a NotImplementedError
EDIT: small modif in the code self != self.base_ring()
> self.parent() != self.base_ring()
comment:45 Changed 3 years ago by
 Commit changed from 5bdbd3232fb25bda9de40d7a8189036eb80b0c19 to 6bf385b18875572eeb4ab748eb8f25124bce6acf
comment:46 in reply to: ↑ 44 Changed 3 years ago by
 Status changed from needs_info to needs_review
Replying to vdelecroix:
It is fine for me if you modify the generic
nth_root
tou = f.unit() if u.is_one(): # do nothing elif self.parent() != self.base_ring(): # try to factorize the unit in the base ring else: # raise a NotImplementedErrorEDIT: small modif in the code
self != self.base_ring()
>self.parent() != self.base_ring()
Thanks for the suggestion! I've pushed some changes such that the unit is handled with more care. Apart from that, I've reviewed the documentation changes and doctests already in the past, this is still fine for me. Please crossreview. :)
comment:47 Changed 3 years ago by
With
if u.parent() != u.base_ring(): u = u.base_ring(u) try: ans = u.nth_root(n) except AttributeError: raise NotImplementedError(...)
You have a risk of infinite recursion. You should use nth_root
of u
only if u.parent()
is not self.parent()
anymore.
comment:48 Changed 3 years ago by
Moreover, it is not always possible to do u = u.base_ring(u)
(even for units).
comment:49 Changed 3 years ago by
 Branch changed from u/behackl/polynomial/rationalpowers to u/public/20086
 Commit 6bf385b18875572eeb4ab748eb8f25124bce6acf deleted
 Milestone changed from sage7.1 to sage7.2
comment:50 followup: ↓ 52 Changed 3 years ago by
 Branch changed from u/public/20086 to public/20086
 Commit set to 654e207cd5b003c4bae5edbde67c5b49330646f8
I added a commit to fix the issues I mentioned. As a reviewer, I think that the code lacks example. It currently only tests ZZ[x]
but is intended to work for any unique factorization domain...
Last 10 new commits:
9bda366  Trac #20086: refactor as try/except/else

cba01eb  Trac 20086: implement (polynomial)^(rational)

c603c2b  Trac 20086: improvements / simplification

42c5319  Trac 20086: fix doctest in polynomial_quotient_ring.py

4313d3f  Merge tag '7.2.beta1' into polynomial/rationalpowers

c8350c6  improve exception messages

5bdbd32  special treatment for zero polynomial

88dced3  Merge tag '7.2.beta6' into polynomial/rationalpowers

6bf385b  handle unit more carefully

654e207  Trac 20086: fix some issues with nth_root

comment:51 Changed 3 years ago by
 Commit changed from 654e207cd5b003c4bae5edbde67c5b49330646f8 to 49fe88e4a5f7058abd072f36cba075b698586841
comment:52 in reply to: ↑ 50 Changed 3 years ago by
Replying to vdelecroix:
I added a commit to fix the issues I mentioned. As a reviewer, I think that the code lacks example. It currently only tests
ZZ[x]
but is intended to work for any unique factorization domain...
I wasn't too happy with the repeated if
and raise
statements, so I tried to clean that up a bit. Also, I've added doctests for QQ[x]
, multivariate polynomial rings, and the number field in x^2  2
. Does anything more exotic come to your mind that you would like to have tested as well?
comment:53 followup: ↓ 55 Changed 3 years ago by
Would be nice to have a non polynomial examples. But currently, order in number fields does not know whether they are principal ideal domain or unique factorization domain
sage: R = ZZ[I] sage: R in PrincipalIdealDomains() False sage: R in UniqueFactorizationDomains() False
In __pow__
for integer polynomials, in the case the exponent is an integer you should only use nn
and not exp
(i.e. you should replace if exp == 0
and if exp < 0
by if nn == 0
and nn < 0
).
comment:54 Changed 3 years ago by
 Commit changed from 49fe88e4a5f7058abd072f36cba075b698586841 to 6f91df97a81fa094c04583314ad38cd5fd199cdb
Branch pushed to git repo; I updated commit sha1. New commits:
6f91df9  exp > nn

comment:55 in reply to: ↑ 53 Changed 3 years ago by
Replying to vdelecroix:
Would be nice to have a non polynomial examples. But currently, order in number fields does not know whether they are principal ideal domain or unique factorization domain
sage: R = ZZ[I] sage: R in PrincipalIdealDomains() False sage: R in UniqueFactorizationDomains() False
Yes, the gaussian integers were also my first idea for such an example, but getting this to work should rather be a followup ticket.
In
__pow__
for integer polynomials, in the case the exponent is an integer you should only usenn
and notexp
(i.e. you should replaceif exp == 0
andif exp < 0
byif nn == 0
andnn < 0
).
Makes sense. Done.
comment:56 followup: ↓ 57 Changed 3 years ago by
 Reviewers changed from Benjamin Hackl to Benjamin Hackl, Vincent Delecroix
 Status changed from needs_review to positive_review
Enough to go!
comment:57 in reply to: ↑ 56 Changed 3 years ago by
Replying to vdelecroix:
Enough to go!
Thanks for this final sprint! Now I can sleep a bit better... ;)
comment:58 followup: ↓ 59 Changed 3 years ago by
 Status changed from positive_review to needs_info
Taking an nth root of an element in a ring R by computing a factorization is an insane way to go about it. A much better generic strategy is to hope that the univariate polynomial ring over R has a root finding algorithm and see if the polynomial x^na
has a root. You will see that:
 it actually has a decent performance over QQ (although there the algorithm should really be specialcased)
 it will work over most fields, including the ones that are not constructed as fraction fields of rings with a factorization algorithm.
 you don't have to mess around with the unit part that a factorization algorithm probably won't recognize.
Illustration:
sage: a=next_prime(10^10)*next_prime(10^11) sage: b=a^2 sage: %timeit (R.0^2b).roots() 1000 loops, best of 3: 428 µs per loop sage: %timeit b.factor() 100 loops, best of 3: 3.68 ms per loop
sage: k.<r>=NumberField(x^2+5) sage: b=k(7^2) sage: (k['x'].0^2b).roots() [(7, 1), (7, 1)] sage: factor(b) #doesn't work of course
comment:59 in reply to: ↑ 58 ; followup: ↓ 60 Changed 3 years ago by
Replying to nbruin:
Taking an nth root of an element in a ring R by computing a factorization is an insane way to go about it. A much better generic strategy is to hope that the univariate polynomial ring over R has a root finding algorithm and see if the polynomial
x^na
has a root. You will see that:
 it actually has a decent performance over QQ (although there the algorithm should really be specialcases)
 it will work over most fields, including the ones that are not constructed as fraction fields of rings with a factorization algorithm.
 you don't have to mess around with the unit part that a factorization algorithm probably won't recognize.
While I like the general idea of this approach, I'm for discussing it on a followup ticket; see it as a "performance enhancement" of this implementation. Also, I'm not sure how exactly the rootfinding algorithm for these univariate polynomial rings can be properly motivated, for me neither roots()
nor any_root()
did the job:
sage: R.<x> = QQ[] sage: P.<X> = R[] sage: a = X^2  (x^2 + 2*x + 1) sage: a.any_root(R) Traceback (most recent call last): ... TypeError: Unable to coerce Principal ideal (1) of Univariate Polynomial Ring in x over Rational Field (<class 'sage.rings.polynomial.ideal.Ideal_1poly_field'>) to Rational
Do you strongly object against setting this back to positive_review
?
comment:60 in reply to: ↑ 59 ; followup: ↓ 61 Changed 3 years ago by
Replying to behackl:
Do you strongly object against setting this back to
positive_review
?
Yes, because it is parking code in the wrong spot.
Perhaps park your code on sage.rings.polynomial.polynomial_element.Polynomial
? Taking an nth root of a polynomial via factorization isn't quite as bad as in general.
Since NumberFieldElement
etc. already provides an nth_root
method, it may be sufficient to just provide the method on Polynomial.
comment:61 in reply to: ↑ 60 Changed 3 years ago by
Replying to nbruin:
Replying to behackl:
Do you strongly object against setting this back to
positive_review
?Yes, because it is parking code in the wrong spot.
I was hoping for a comment from Vincent, as he started with the generic solution on this branch.
In any case: I do not quite understand why this would be in the wrong place. Unique factorization domains should provide a factor
method, and providing a generic nth_root
method for them doesn't strike me as bad.
It might be true that this ticket only changes the behavior of polynomial rings, but I don't see a downside of providing a generic solution. But maybe I'm missing something? :)
Perhaps park your code on
sage.rings.polynomial.polynomial_element.Polynomial
? Taking an nth root of a polynomial via factorization isn't quite as bad as in general.Since
NumberFieldElement
etc. already provides annth_root
method, it may be sufficient to just provide the method on Polynomial.
Yes, that is truehowever, I think that the code provided on this ticket is generic enough to work for all unique factorization domains. Or doesn't it? Even if it might be not performant, I think that a generic approach is certainly allowed to be slowthis problem can be handled by special case implementations.
I just don't see a benefit from moving the code back to polynomials only. Could you elaborate more why this generic solution should be degraded to a special case?
comment:62 followup: ↓ 63 Changed 3 years ago by
The code provided is far to be working on any UFD as you have to factor the unit! But we know for sure that it will work for polynomial rings whose base ring provides a nth_root
method. It would make sense to move the code to generic polynomial. As there is not a lot of examples of UFD in Sage it is hard to say that this code is useful for general UFD.
comment:63 in reply to: ↑ 62 Changed 3 years ago by
 Status changed from needs_info to needs_work
Replying to vdelecroix:
The code provided is far to be working on any UFD as you have to factor the unit!
I might still be thinking too much of polynomial rings when thinking of UFDs. ;)
But we know for sure that it will work for polynomial rings whose base ring provides a
nth_root
method. It would make sense to move the code to generic polynomial.
The argument regarding the unit and the fact that there are not that much (exotic) UFDs implemented in Sage convice me; thanks for the clarification! :) I'll move the code.
comment:64 followup: ↓ 65 Changed 3 years ago by
Moving the code to sage.rings.polynomial.polynomial_element.Polynomial
only provides nth_root
for univariate polynomial rings, as it seems. The doctests for multivariate rings fail with
AttributeError: 'sage.rings.polynomial.multi_polynomial_libsingular.MPolynomial_libsingular' object has no attribute 'nth_root'
I'm not a particularly big fried of copying the code to sage.rings.polynomial.multi_polynomial_element.MPolynomial_element
as well. The univariate and multivariate polynomial ring elements inherit from CommutativeAlgebraelement? and CommutativeRingElement?, respectivelyso implementing it in a superclass is not possible.
Ideas?
comment:65 in reply to: ↑ 64 Changed 3 years ago by
Replying to behackl:
Ideas?
 Leave the code in
categories.unique_factorization_domain
, raise aNotImplementedError
... ... if something goes wrong (current behavior), or
 ... when the calling parent is not a polynomial ring (that seems relatively restrictive to me)
 Move it to either univariate or multivariate polynomials, call the same method from the other one (I'm thinking of something like
nth_root = sage.rings.polynomial.polynomial_element.nth_root
or so...)  Move it to both univaraite and multivariate polynomial rings
I'll try to do the second one and push the branch again, let me know if you feel that there is a better solution.
comment:66 Changed 3 years ago by
 Commit changed from 6f91df97a81fa094c04583314ad38cd5fd199cdb to 7940c9f126743aadb02ef8325dd474f9986592ce
Branch pushed to git repo; I updated commit sha1. New commits:
7940c9f  duplicate code for univariate/multivariate polynomial rings

comment:67 Changed 3 years ago by
 Status changed from needs_work to needs_review
My attempt to implement the second strategy failed spectacularly with
TypeError: descriptor 'nth_root' for 'sage.rings.polynomial.polynomial_element.Polynomial' objects doesn't apply to 'sage.rings.polynomial.multi_polynomial_libsingular.MPolynomial_libsingular' object
so I've duplicated the code and adapted the doctests a bit.
I'm not particularly happy with this solution, but I don't see a better one*. Suggestions are very welcome. Back to needs_review
again...
*EDIT: that doesn't live in categories.unique_factorization_domains
. There seems to be nothing sufficiently in common between the univariate and multivariate polynomial rings.
comment:68 Changed 3 years ago by
 Commit changed from 7940c9f126743aadb02ef8325dd474f9986592ce to 0fc07442120a84047511625a8fe4bbd1cef5bb60
Branch pushed to git repo; I updated commit sha1. New commits:
0fc0744  add remark regarding duplicated code

comment:69 Changed 3 years ago by
 Commit changed from 0fc07442120a84047511625a8fe4bbd1cef5bb60 to 8cef2c7bce3fc7379f95d6c7c3bd6fb21e64c0ab
Branch pushed to git repo; I updated commit sha1. New commits:
8cef2c7  fix a doctest in polynomial_quotient_ring.py

comment:70 Changed 3 years ago by
I've added comments that there is a duplicated version of the code. Also, I've fixed the failing doctest.
Even though we are too late for 7.2
, I'd still greatly appreciate if someone could review my last few changes. The code itself is already reviewed.
comment:71 Changed 3 years ago by
Since now the code is in the polynomial ring you can simplify it a lot. Just remove all of
u = f.unit() if u.is_one(): ... else: ....
and do
u = self.base_ring()(f.unit()) try: u.nth_root(n) except AttributeError: raise NotImplementedError
comment:72 Changed 3 years ago by
And if the degree is not divisible by n
you already know that there is no nth root (in multivariate case, degree="sum of degrees of monomials"). This is very cheap to test.
comment:73 followup: ↓ 74 Changed 3 years ago by
 Commit changed from 8cef2c7bce3fc7379f95d6c7c3bd6fb21e64c0ab to 753462d07465a87797609856bec781797da9a6b1
Branch pushed to git repo; I updated commit sha1. New commits:
753462d  more refactoring after moving code to rings.polynomial

comment:74 in reply to: ↑ 73 Changed 3 years ago by
comment:75 followup: ↓ 77 Changed 3 years ago by
Perhaps don't bother getting a grammatically correct ordinal. Note that the current code doesn't get it correct anyway: It's 21st, 22nd, 23rd.
Better to formulate the error message in a way that doesn't depend on correct ordinal spelling, e.g.
ValueError("(%s)^(1/%s) does not lie in ring"%(f,n))
Concerning further optimization: testing degree is of course a worthwhile cheap trick. A followup ticket should probably use squarefree factorization; something along the lines:
 if the characteristic p divides n then first check that the polynomial only has pth powers of the variables in it. Take pth root (i.e., replace variables and take pth root of coefficients) and take (n/p)th root of resulting polynomial
 In characteristic 0 then take g=f/GCD(f,f.derivative()), check that
g^n
dividesf
, take nth root off/(g^n)
, multiply by g and return that.
In positive characteristic we may first need to take all pth powers out.
comment:76 Changed 3 years ago by
 Commit changed from 753462d07465a87797609856bec781797da9a6b1 to e56adf4325379d05441efb4dd9487c2cc7904ddc
Branch pushed to git repo; I updated commit sha1. New commits:
e56adf4  fix ordinals

comment:77 in reply to: ↑ 75 ; followup: ↓ 78 Changed 3 years ago by
Replying to nbruin:
Perhaps don't bother getting a grammatically correct ordinal. Note that the current code doesn't get it correct anyway: It's 21st, 22nd, 23rd.
Better to formulate the error message in a way that doesn't depend on correct ordinal spelling, e.g.
ValueError("(%s)^(1/%s) does not lie in ring"%(f,n))
I do like the version with explicit ordinals in the errors more, so I fixed them. :)
Concerning further optimization: testing degree is of course a worthwhile cheap trick. A followup ticket should probably use squarefree factorization; something along the lines:
 if the characteristic p divides n then first check that the polynomial only has pth powers of the variables in it. Take pth root (i.e., replace variables and take pth root of coefficients) and take (n/p)th root of resulting polynomial
 In characteristic 0 then take g=f/GCD(f,f.derivative()), check that
g^n
dividesf
, take nth root off/(g^n)
, multiply by g and return that.In positive characteristic we may first need to take all pth powers out.
Of course, there is always potential to improve the performance. However, like you said, this should happen in a followup ticket.
comment:78 in reply to: ↑ 77 Changed 3 years ago by
Replying to behackl:
Replying to nbruin:
Perhaps don't bother getting a grammatically correct ordinal. Note that the current code doesn't get it correct anyway: It's 21st, 22nd, 23rd.
Better to formulate the error message in a way that doesn't depend on correct ordinal spelling, e.g.
ValueError("(%s)^(1/%s) does not lie in ring"%(f,n))
I do like the version with explicit ordinals in the errors more, so I fixed them. :)
Having the following code executed each time the function run is useless
if 10 <= n % 100 < 20: postfix = 'th' else: postfix = {1:'st', 2:'nd', 3:'rd'}.get(n % 10, 'th')
Moreover, having these four lines potentially anywhere in Sage because you like is not the way to go. If you really think it is better this way, then implement a function ordinal_str
that would really be doctested.
In the error message proposed by Nils there was the explicit mention of the parent which is a good thing.
comment:79 Changed 3 years ago by
 Commit changed from e56adf4325379d05441efb4dd9487c2cc7904ddc to e8adfb10e935659d9701aa69b888571b569e222a
Branch pushed to git repo; I updated commit sha1. New commits:
e8adfb1  remove ordinals, improve error messages

comment:80 followup: ↓ 81 Changed 3 years ago by
 Status changed from needs_review to needs_work
sage: R.<x> = ZZ[] sage: parent(R.one().nth_root(3)) Integer Ring
comment:81 in reply to: ↑ 80 ; followup: ↓ 83 Changed 3 years ago by
Replying to vdelecroix:
sage: R.<x> = ZZ[] sage: parent(R.one().nth_root(3)) Integer Ring
Fixed, will push in a moment.
comment:82 Changed 3 years ago by
 Commit changed from e8adfb10e935659d9701aa69b888571b569e222a to 2c46e087a34666657dee8cb2f168a08d8afad488
Branch pushed to git repo; I updated commit sha1. New commits:
2c46e08  return one from original parent, not base ring

comment:83 in reply to: ↑ 81 Changed 3 years ago by
 Status changed from needs_work to needs_review
comment:84 followup: ↓ 88 Changed 3 years ago by
Patchbot reports doctest failure due to the change in the error message.
comment:85 Changed 3 years ago by
 Description modified (diff)
 Status changed from needs_review to needs_work
I have an implementation for a much faster implementation of nth root at #20571.
comment:86 Changed 3 years ago by
 Description modified (diff)
comment:87 Changed 3 years ago by
 Description modified (diff)
comment:88 in reply to: ↑ 84 Changed 3 years ago by
Replying to vdelecroix:
Patchbot reports doctest failure due to the change in the error message.
Indeed, and this can be fixed straightforward. However, there is (once again ...) a slight inconvenience: there are superfluous parentheses in, e.g.
sage: P.<x> = ZZ[] sage: P(4).nth_root(3) Traceback (most recent call last): ... ValueError: (4)^(1/3) does not lie in ...
or
sage: x.nth_root(3) Traceback (most recent call last): ... ValueError: (x)^(1/3) does not lie in ...
Should we live with that? Fixing it would require querying something like self.is_constant()
and self in self.parent().gens()
or so.
EDIT: self.is_constant()
would not work:
sage: P.<x> = ZZ[] sage: Q.<y> = P[] sage: Q(x+1).is_constant() True
comment:89 followup: ↓ 91 Changed 3 years ago by
I would think most people would prioritize #20571 over such typographical issues, so I wouldn't sweat too hard over parentheses at this point. Another issue: while it's nice to have informative error messages, error strings that cost a lot to be constructed only to be caught by an "except" statement can mean a performance penalty. This is more an issue in python than in many other languages because there's a culture in python (and perhaps even more so in some parts of sage) to use exceptions for regular program flow control, not just for error conditions.
comment:90 Changed 3 years ago by
 Commit changed from 2c46e087a34666657dee8cb2f168a08d8afad488 to 5fc1027e43d601a791568b22f12ec70957c47519
Branch pushed to git repo; I updated commit sha1. New commits:
5fc1027  fix doctests

comment:91 in reply to: ↑ 89 Changed 3 years ago by
Replying to nbruin:
I would think most people would prioritize #20571 over such typographical issues, so I wouldn't sweat too hard over parentheses at this point. Another issue: while it's nice to have informative error messages, error strings that cost a lot to be constructed only to be caught by an "except" statement can mean a performance penalty. This is more an issue in python than in many other languages because there's a culture in python (and perhaps even more so in some parts of sage) to use exceptions for regular program flow control, not just for error conditions.
Yes, that's my feeling too. While it is somehow unfortunate, I can very well live with it.
I've fixed the doctests and tested the entire src/sage/rings/polynomial
directory; let's see what the patchbot thinks.
comment:92 Changed 3 years ago by
 Status changed from needs_work to needs_review
comment:93 Changed 3 years ago by
 Status changed from needs_review to positive_review
Let's move to #20571.
comment:94 Changed 3 years ago by
 Branch changed from public/20086 to 5fc1027e43d601a791568b22f12ec70957c47519
 Resolution set to fixed
 Status changed from positive_review to closed
Branch pushed to git repo; I updated commit sha1. New commits:
Trac #20086: powers of constant polynomials