Opened 4 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: sage-7.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 vdelecroix)

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

  • Branch set to u/cheuberg/polynomials/power

comment:2 Changed 4 years ago by git

  • Commit set to cc49098e1fc111e63d344b6c2f914d0d3e5e463b

Branch pushed to git repo; I updated commit sha1. New commits:

cc49098Trac #20086: powers of constant polynomials

comment:3 Changed 4 years ago by cheuberg

  • Authors set to Clemens Heuberger
  • 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 4 years ago by git

  • Commit changed from cc49098e1fc111e63d344b6c2f914d0d3e5e463b to 25fdd0c0ceb085d68b2ddf851983a386771d164f

Branch pushed to git repo; I updated commit sha1. New commits:

25fdd0cTrac #20086: powers of constant polynomials over ZZ

comment:5 Changed 4 years ago by cheuberg

  • Summary changed from QQ[X]: allow arbitrary powers of constant polynomials to ZZ[X], QQ[X]: allow arbitrary powers of constant polynomials

comment:6 follow-up: Changed 4 years ago by 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 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 4 years ago by git

  • Commit changed from 25fdd0c0ceb085d68b2ddf851983a386771d164f to 6217beeb2520f5e8c0d1e92881c4d1c9569e9868

Branch pushed to git repo; I updated commit sha1. New commits:

6217beeTrac #20085: raise exception instead of moving to SR

comment:8 Changed 4 years ago by git

  • Commit changed from 6217beeb2520f5e8c0d1e92881c4d1c9569e9868 to 918499cf377a551a5e28ea1d6748ed50860f1afa

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

918499cTrac #20086: raise exception instead of moving to SR

comment:9 in reply to: ↑ 6 Changed 4 years ago by cheuberg

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

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 should QQ['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 follow-up ticket and to have this functionality here soon.

comment:10 Changed 4 years ago by behackl

  • Cc dkrenn added; dakrenn removed

comment:11 follow-up: Changed 4 years ago by behackl

  • 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 follow-up ticket; I'd also like to have this in as soon as possible.

I've reviewed your changes, please cross-review. If you are satisfied and if there are no other objections, I'd set this to positive_review.


New commits:

caa02c8Merge tag '7.1.beta6' into HEAD
3f44399fix segmentation fault
5934ed1fix (constant poly)^(constant poly)
e6e1c81add doctests

comment:12 Changed 4 years ago by behackl

  • Reviewers set to Benjamin Hackl

comment:13 Changed 4 years ago by git

  • Commit changed from e6e1c812ecb49be56dc4d9b4ed2e524011c391c9 to 87a22b6fc4cff82491a7b0ba35983bc5dc1299df

Branch pushed to git repo; I updated commit sha1. New commits:

87a22b6fix "referenced before assignment"

comment:14 Changed 4 years ago by cheuberg

  • Branch changed from u/behackl/polynomials/power to u/cheuberg/polynomials/power

comment:15 in reply to: ↑ 11 Changed 4 years ago by cheuberg

  • Commit changed from 87a22b6fc4cff82491a7b0ba35983bc5dc1299df to 14f7efaaf5be81405d1954473b4dc3f47100e486

Replying to behackl:

I've fixed the segmentation fault from rings/integer.pyx and implemented the case of constant^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 cross-review. 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 re-instated fix (in some other way).


New commits:

e95c9b9Trac #20086: refactor method; do not fix segmentation fault 2^R in polynomial ring
eadecbaTrac #20086: fix ReSt errors (docbuild -u -k complains)
a139868Trac #20086: handle non-constant polynomial exponents
14f7efaTrac #20086: ReSt formatting

comment:16 follow-up: Changed 4 years ago by vdelecroix

  • 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 ; follow-up: Changed 4 years ago by cheuberg

  • 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)
-> 2

Which 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 ; follow-ups: Changed 4 years ago by vdelecroix

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 exponentiation

I 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)
-> 2

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

  • Description modified (diff)
  • Status changed from needs_review to needs_work

It seems that I mis-edited 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 4 years ago by cheuberg

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

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 ticket---but then, this could be extended easily.

Benjamin

comment:22 Changed 4 years ago by git

  • Commit changed from 14f7efaaf5be81405d1954473b4dc3f47100e486 to 892109ac4f48a4a32d35671365ecceb0d1c9e0bf

Branch pushed to git repo; I updated commit sha1. New commits:

892109aTrac #20086: fix bug in previous commit (and use better doctest)

comment:23 Changed 4 years ago by cheuberg

  • 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 ; follow-up: Changed 4 years ago by cheuberg

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 sage-devel a few weeks ago. Every parent seems to have its own home-made logic about how to do coercion. Please do not try to fix that in this ticket.

comment:25 in reply to: ↑ 24 ; follow-up: Changed 4 years ago by vdelecroix

Replying to cheuberg:

[... snip ...]

There was a discussion on sage-devel a few weeks ago. Every parent seems to have its own home-made 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 a NotImplementedError.
  • disallow polynomial^polynomial even if the power is a constant polynomial. In that case, the error should be a TypeError. This is the current behavior for some polynomial rings but not all
    sage: 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 4 years ago by cheuberg

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 a NotImplementedError.
  • disallow polynomial^polynomial even if the power is a constant polynomial. In that case, the error should be a TypeError. This is the current behavior for some polynomial rings but not all
    sage: 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?

As far as I can see without implementing it and make ptestlong, this fits perfectly.

comment:27 Changed 4 years ago by git

  • Commit changed from 892109ac4f48a4a32d35671365ecceb0d1c9e0bf to 1782a3cc92f875fa3ac8a83395d204d2ca3d7f8c

Branch pushed to git repo; I updated commit sha1. New commits:

1782a3cTrac #20086: rewrite try/expect; disallow polynomial exponents

comment:28 follow-up: Changed 4 years ago by cheuberg

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 ; follow-up: Changed 4 years ago by vdelecroix

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 correct TypeError. 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 4 years ago by git

  • Commit changed from 1782a3cc92f875fa3ac8a83395d204d2ca3d7f8c to 9bda3669c14463b3061f19c959bf8a22cb467117

Branch pushed to git repo; I updated commit sha1. New commits:

9bda366Trac #20086: refactor as try/except/else

comment:31 in reply to: ↑ 29 Changed 4 years ago by cheuberg

Replying to vdelecroix:

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

done.

comment:32 Changed 4 years ago by vdelecroix

  • Authors changed from Clemens Heuberger to Clemens Heuberger, Vincent Delecroix
  • 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:

cba01ebTrac 20086: implement (polynomial)^(rational)

comment:33 Changed 4 years ago by vdelecroix

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

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

  • Commit changed from cba01eb6ebe741567cabc1e9e4a3dd0c94aa11b2 to c603c2b293819088246b34ab676f303871d5cb40

Branch pushed to git repo; I updated commit sha1. New commits:

c603c2bTrac 20086: improvements / simplification

comment:36 Changed 4 years ago by git

  • Commit changed from c603c2b293819088246b34ab676f303871d5cb40 to 42c53199edeb619d8b4867cea0b318d7d0c8ea0a

Branch pushed to git repo; I updated commit sha1. New commits:

42c5319Trac 20086: fix doctest in polynomial_quotient_ring.py

comment:37 Changed 4 years ago by behackl

  • Branch changed from u/vdelecroix/20086 to u/behackl/polynomial/rational-powers
  • 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 cross-review and set to positive_review if you are satisfied.


New commits:

4313d3fMerge tag '7.2.beta1' into polynomial/rational-powers
c8350c6improve exception messages

comment:38 Changed 4 years ago by vdelecroix

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() is True (or more generally when u.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 4 years ago by vdelecroix

  • Status changed from needs_review to needs_work

comment:40 Changed 4 years ago by behackl

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 method---and 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 non-constant 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 follow-up 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 4 years ago by git

  • Commit changed from c8350c631287dd71a9b79fb422d769a1e480cdc8 to 5bdbd3232fb25bda9de40d7a8189036eb80b0c19

Branch pushed to git repo; I updated commit sha1. New commits:

5bdbd32special treatment for zero polynomial

comment:42 follow-up: Changed 3 years ago by behackl

... any new thoughts on this?

comment:43 in reply to: ↑ 42 Changed 3 years ago by behackl

  • 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 follow-up: Changed 3 years ago by vdelecroix

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

Last edited 3 years ago by vdelecroix (previous) (diff)

comment:45 Changed 3 years ago by git

  • Commit changed from 5bdbd3232fb25bda9de40d7a8189036eb80b0c19 to 6bf385b18875572eeb4ab748eb8f25124bce6acf

Branch pushed to git repo; I updated commit sha1. New commits:

88dced3Merge tag '7.2.beta6' into polynomial/rational-powers
6bf385bhandle unit more carefully

comment:46 in reply to: ↑ 44 Changed 3 years ago by behackl

  • Status changed from needs_info to needs_review

Replying to vdelecroix:

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

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 cross-review. :-)

comment:47 Changed 3 years ago by vdelecroix

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 vdelecroix

Moreover, it is not always possible to do u = u.base_ring(u) (even for units).

comment:49 Changed 3 years ago by vdelecroix

  • Branch changed from u/behackl/polynomial/rational-powers to u/public/20086
  • Commit 6bf385b18875572eeb4ab748eb8f25124bce6acf deleted
  • Milestone changed from sage-7.1 to sage-7.2

comment:50 follow-up: Changed 3 years ago by vdelecroix

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

9bda366Trac #20086: refactor as try/except/else
cba01ebTrac 20086: implement (polynomial)^(rational)
c603c2bTrac 20086: improvements / simplification
42c5319Trac 20086: fix doctest in polynomial_quotient_ring.py
4313d3fMerge tag '7.2.beta1' into polynomial/rational-powers
c8350c6improve exception messages
5bdbd32special treatment for zero polynomial
88dced3Merge tag '7.2.beta6' into polynomial/rational-powers
6bf385bhandle unit more carefully
654e207Trac 20086: fix some issues with nth_root

comment:51 Changed 3 years ago by git

  • Commit changed from 654e207cd5b003c4bae5edbde67c5b49330646f8 to 49fe88e4a5f7058abd072f36cba075b698586841

Branch pushed to git repo; I updated commit sha1. New commits:

fb17388remove double if, minor adaptions
49fe88eadd more doctests

comment:52 in reply to: ↑ 50 Changed 3 years ago by behackl

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 follow-up: Changed 3 years ago by 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

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 git

  • Commit changed from 49fe88e4a5f7058abd072f36cba075b698586841 to 6f91df97a81fa094c04583314ad38cd5fd199cdb

Branch pushed to git repo; I updated commit sha1. New commits:

6f91df9exp --> nn

comment:55 in reply to: ↑ 53 Changed 3 years ago by behackl

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 follow-up ticket.

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

Makes sense. Done.

comment:56 follow-up: Changed 3 years ago by vdelecroix

  • Authors changed from Clemens Heuberger, Vincent Delecroix to Clemens Heuberger, Vincent Delecroix, Benjamin Hackl
  • 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 behackl

Replying to vdelecroix:

Enough to go!

Thanks for this final sprint! Now I can sleep a bit better... ;-)

comment:58 follow-up: Changed 3 years ago by nbruin

  • Status changed from positive_review to needs_info

Taking an n-th 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^n-a has a root. You will see that:

  • it actually has a decent performance over QQ (although there the algorithm should really be special-cased)
  • 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^2-b).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^2-b).roots()
[(7, 1), (-7, 1)]
sage: factor(b) #doesn't work of course
Last edited 3 years ago by nbruin (previous) (diff)

comment:59 in reply to: ↑ 58 ; follow-up: Changed 3 years ago by behackl

Replying to nbruin:

Taking an n-th 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^n-a has a root. You will see that:

  • it actually has a decent performance over QQ (although there the algorithm should really be special-cases)
  • 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 follow-up ticket; see it as a "performance enhancement" of this implementation. Also, I'm not sure how exactly the root-finding 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 ; follow-up: Changed 3 years ago by 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.

Perhaps park your code on sage.rings.polynomial.polynomial_element.Polynomial? Taking an n-th 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 behackl

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 n-th 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.

Yes, that is true---however, 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 slow---this 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?

Last edited 3 years ago by behackl (previous) (diff)

comment:62 follow-up: Changed 3 years ago by vdelecroix

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 behackl

  • 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 follow-up: Changed 3 years ago by behackl

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?, respectively---so implementing it in a superclass is not possible.

Ideas?

comment:65 in reply to: ↑ 64 Changed 3 years ago by behackl

Replying to behackl:

Ideas?

  • Leave the code in categories.unique_factorization_domain, raise a NotImplementedError ...
    • ... 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 git

  • Commit changed from 6f91df97a81fa094c04583314ad38cd5fd199cdb to 7940c9f126743aadb02ef8325dd474f9986592ce

Branch pushed to git repo; I updated commit sha1. New commits:

7940c9fduplicate code for univariate/multivariate polynomial rings

comment:67 Changed 3 years ago by behackl

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

Last edited 3 years ago by behackl (previous) (diff)

comment:68 Changed 3 years ago by git

  • Commit changed from 7940c9f126743aadb02ef8325dd474f9986592ce to 0fc07442120a84047511625a8fe4bbd1cef5bb60

Branch pushed to git repo; I updated commit sha1. New commits:

0fc0744add remark regarding duplicated code

comment:69 Changed 3 years ago by git

  • Commit changed from 0fc07442120a84047511625a8fe4bbd1cef5bb60 to 8cef2c7bce3fc7379f95d6c7c3bd6fb21e64c0ab

Branch pushed to git repo; I updated commit sha1. New commits:

8cef2c7fix a doctest in polynomial_quotient_ring.py

comment:70 Changed 3 years ago by behackl

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 vdelecroix

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 vdelecroix

And if the degree is not divisible by n you already know that there is no n-th root (in multivariate case, degree="sum of degrees of monomials"). This is very cheap to test.

comment:73 follow-up: Changed 3 years ago by git

  • Commit changed from 8cef2c7bce3fc7379f95d6c7c3bd6fb21e64c0ab to 753462d07465a87797609856bec781797da9a6b1

Branch pushed to git repo; I updated commit sha1. New commits:

753462dmore refactoring after moving code to rings.polynomial

comment:74 in reply to: ↑ 73 Changed 3 years ago by behackl

Replying to git:

Branch pushed to git repo; I updated commit sha1. New commits:

753462dmore refactoring after moving code to rings.polynomial

Simplified the code and added the additional check for self.degree() % n.

comment:75 follow-up: Changed 3 years ago by 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))

Concerning further optimization: testing degree is of course a worthwhile cheap trick. A follow-up ticket should probably use square-free factorization; something along the lines:

  • if the characteristic p divides n then first check that the polynomial only has p-th powers of the variables in it. Take p-th root (i.e., replace variables and take p-th 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 divides f, take n-th root of f/(g^n), multiply by g and return that.

In positive characteristic we may first need to take all p-th powers out.

Last edited 3 years ago by nbruin (previous) (diff)

comment:76 Changed 3 years ago by git

  • Commit changed from 753462d07465a87797609856bec781797da9a6b1 to e56adf4325379d05441efb4dd9487c2cc7904ddc

Branch pushed to git repo; I updated commit sha1. New commits:

e56adf4fix ordinals

comment:77 in reply to: ↑ 75 ; follow-up: Changed 3 years ago by 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. :-)

Concerning further optimization: testing degree is of course a worthwhile cheap trick. A follow-up ticket should probably use square-free factorization; something along the lines:

  • if the characteristic p divides n then first check that the polynomial only has p-th powers of the variables in it. Take p-th root (i.e., replace variables and take p-th 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 divides f, take n-th root of f/(g^n), multiply by g and return that.

In positive characteristic we may first need to take all p-th powers out.

Of course, there is always potential to improve the performance. However, like you said, this should happen in a follow-up ticket.

comment:78 in reply to: ↑ 77 Changed 3 years ago by vdelecroix

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 git

  • Commit changed from e56adf4325379d05441efb4dd9487c2cc7904ddc to e8adfb10e935659d9701aa69b888571b569e222a

Branch pushed to git repo; I updated commit sha1. New commits:

e8adfb1remove ordinals, improve error messages

comment:80 follow-up: Changed 3 years ago by vdelecroix

  • 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 ; follow-up: Changed 3 years ago by behackl

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 git

  • Commit changed from e8adfb10e935659d9701aa69b888571b569e222a to 2c46e087a34666657dee8cb2f168a08d8afad488

Branch pushed to git repo; I updated commit sha1. New commits:

2c46e08return one from original parent, not base ring

comment:83 in reply to: ↑ 81 Changed 3 years ago by behackl

  • Status changed from needs_work to needs_review

comment:84 follow-up: Changed 3 years ago by vdelecroix

Patchbot reports doctest failure due to the change in the error message.

comment:85 Changed 3 years ago by vdelecroix

  • Description modified (diff)
  • Status changed from needs_review to needs_work

I have an implementation for a much faster implementation of n-th root at #20571.

Last edited 3 years ago by vdelecroix (previous) (diff)

comment:86 Changed 3 years ago by vdelecroix

  • Description modified (diff)

comment:87 Changed 3 years ago by vdelecroix

  • Description modified (diff)

comment:88 in reply to: ↑ 84 Changed 3 years ago by behackl

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
Last edited 3 years ago by behackl (previous) (diff)

comment:89 follow-up: Changed 3 years ago by 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.

comment:90 Changed 3 years ago by git

  • Commit changed from 2c46e087a34666657dee8cb2f168a08d8afad488 to 5fc1027e43d601a791568b22f12ec70957c47519

Branch pushed to git repo; I updated commit sha1. New commits:

5fc1027fix doctests

comment:91 in reply to: ↑ 89 Changed 3 years ago by behackl

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 behackl

  • Status changed from needs_work to needs_review

comment:93 Changed 3 years ago by vdelecroix

  • Status changed from needs_review to positive_review

Let's move to #20571.

comment:94 Changed 3 years ago by vbraun

  • Branch changed from public/20086 to 5fc1027e43d601a791568b22f12ec70957c47519
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.