symbolic radical expression for algebraic number

[gagern]

It would be nice to have a function which converts algebraic numbers to symbolic expressions using radicals, at least for those cases where this is possible (i.e. the cases where the minimal polynomial of the algebraic number has a degree less than 5). For the quadratic case, I have a code snippet to illustrate what I'm asking for:

def AA2SR(x):
    x = AA(x)
    p = x.minpoly()
    if p.degree() < 2:
        return SR(QQ(x))
    if p.degree() > 2:
        raise TypeError("Cannot handle degrees > 2")
    c, b, a = p.coeffs()
    y = (-b+sqrt(b*b-4*a*c))/(2*a)
    if x == AA(y):
        return y
    y = (-b-sqrt(b*b-4*a*c))/(2*a)
    assert x == AA(y)
    return y

def QQbar2SR(x):
    x = QQbar(x)
    return AA2SR(x.real()) + I*AA2SR(x.imag())

These functions could then be applied to vectors or matrices over algebraic real or complex numbers in order to obtain an exact description of the relevant values using radicals.

This request is a spin-off from my question on Ask Sage.

Follow-up tickets: #17516, #17517

Depends on #17495
Depends on #16964

Component: number fields

Author: Martin von Gagern, Jeroen Demeyer

Branch/Commit: 4626286

Reviewer: Marc Mezzarobba, Jeroen Demeyer, Vincent Delecroix

Issue created by migration from https://trac.sagemath.org/ticket/14239

[gagern]

Attachment: trac_14239_algebraic_to_radicals.patch.gz

[gagern]

comment:1

I wrapped my code into two methods radical_expression for the AlgebraicNumber and AlgebraicReal classes.

Only second degree polynomials supported so far. I guess this should cover the most useful use cases: third and fourth degree are certainly possible, but don't help very much with the readability of the result in most cases. Square roots, on the other hand, are everywhere, so being able to represent these is a huge win imho.

The code I wrote does rely little on the internal workings of algebraic numbers and their descriptions, mostly because I haven't dug into that code yet. I guess the check which solution of the quadratic equation to choose might benefit from a more direct comparison of separating intervals, but at the moment I see no urgent performance issue with most use cases of this function.

[mezzarobba]

Reviewer: Marc Mezzarobba

[mezzarobba]

comment:3

I think the feature you are looking for sort of exists. But it is available for elements of number fields, not algebraic numbers. And it is pretty inconvenient to use with algebraic numbers because alg.as_number_field_element() returns alg as an element of an abstract number field (plus an embedding of that number field into QQbar or AA), while conversion to symbolic expressions expects a number field with a registered embedding into CC or RR.

Yet, with the example from your question on AskSage, one can do:

sage: a = N[0][0]
sage: nf, val, nf2alg = a.as_number_field_element()
sage: approx = CC(nf2alg(nf.gen()))
sage: embedded_nf = NumberField(nf.polynomial(), 'a', embedding=approx)
sage: nf2nf = sage.rings.number_field.number_field_morphisms.NumberFieldEmbedding(nf, embedded_nf, embedded_nf.gen())
sage: SR(nf2nf(val))
1/5*sqrt(5)

So I disagree with the approach you take in your patch: IMO, we should make something like the above example work automatically and share the actual conversion code between algebraic numbers and number fields instead of implementing essentially the same feature twice.

[edit: add missing word]

[gagern]

comment:5

Replying to @mezzarobba:

Yet, with the example from your question on AskSage, one can do:

I just gave that code a try, and on sage 5.12 it failed in the assignment of embedded_nf with

    embedded_nf = NumberField(nf.polynomial(), 'a', embedding=approx)
  File "parent.pyx", line 761, in sage.structure.parent.Parent.__getattr__ (sage/structure/parent.c:6823)
  File "misc.pyx", line 251, in sage.structure.misc.getattr_from_other_class (sage/structure/misc.c:1606)
AttributeError: 'RationalField_with_category' object has no attribute 'polynomial'

Probably due to an outdated sage, so feel free to ignore this if that's the case. Can't update this system just now.

[mezzarobba]

comment:6

Replying to @gagern:

Probably due to an outdated sage

Yes, I think so. (I tried again on develop from a few days ago, and it works.)

[gagern]

Branch: u/gagern/ticket/14239

[gagern]

Commit: 1ddc68c

[gagern]

New commits:

1ddc68c

Implement method radical_expression for QQbar and AA elements.

[sagetrac-git]

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

30f7109

Check radical expression using QQbar even for AA elements.

[sagetrac-git]

Changed commit from 1ddc68c to 30f7109

[sagetrac-git]

Changed commit from 30f7109 to 039f4bf

[sagetrac-git]

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

039f4bf

Fix radical expression for rational numbers.

[gagern]

comment:11

I notice that in real world applications, the check for self == QQbar(res) takes excessively long. Example:

QQ[x](16*x^6 - 392*x^4 + 133*x^2 + 900).roots(AA, False)[-1].radical_expression()

It gets stuck in several layers of exactify, but eventually the keyboard interrupt ends up here:

  File "sage/rings/qqbar.py", line 6798, in exactify
    red_elt, red_back, red_pol = do_polred(newpol_sage_y)
  File "sage/rings/qqbar.py", line 1645, in do_polred
    new_poly, elt_back = poly._pari_().polredbest(flag=1)
  File "gen.pyx", line 8048, in sage.libs.pari.gen.gen.polredbest (sage/libs/pari/gen.c:42279)
  File "c_lib.pyx", line 73, in sage.ext.c_lib.sig_raise_exception (sage/ext/c_lib.c:872)

On the one hand, I wonder whether there is some better way to detect an inexact conversion. I haven't looked enough at the internal structure of a symbolic expression to know how to handle this.

On the other hand, I wonder why this comparison is taking so long. Or even whether it is simply taking excessively long, or perhaps even got stuck completely. Perhaps this is an indication of some more fundamental problem. If you want to invastigate this, you can do so without checking out this branch. Simply using the example above, the relevant operation is the following:

a = QQ[x](16*x^6 - 392*x^4 + 133*x^2 + 900).roots(AA, False)[-1]
b = -1/1296*((9*I*sqrt(3)*(2/9*I*sqrt(443)*sqrt(3) + 53/27)^(1/3) + 9*(2/9*I*sqrt(443)*sqrt(3) + 53/27)^(1/3) - 37*I*sqrt(3)/(2/9*I*sqrt(443)*sqrt(3) + 53/27)^(1/3) + 37/(2/9*I*sqrt(443)*sqrt(3) + 53/27)^(1/3) + 18)*(9*I*sqrt(3)*(2/9*I*sqrt(443)*sqrt(3) + 53/27)^(1/3) + 9*(2/9*I*sqrt(443)*sqrt(3) + 53/27)^(1/3) - 37*I*sqrt(3)/(2/9*I*sqrt(443)*sqrt(3) + 53/27)^(1/3) + 37/(2/9*I*sqrt(443)*sqrt(3) + 53/27)^(1/3) - 12) - 3456)*sqrt(-9/2*I*sqrt(3)*(2/9*I*sqrt(443)*sqrt(3) + 53/27)^(1/3) - 9/2*(2/9*I*sqrt(443)*sqrt(3) + 53/27)^(1/3) + 37/2*I*sqrt(3)/(2/9*I*sqrt(443)*sqrt(3) + 53/27)^(1/3) - 37/2/(2/9*I*sqrt(443)*sqrt(3) + 53/27)^(1/3) + 6)
a == b

Even as it stands, having the as_radical_expression method is better than not having it. So this comment here should be no reason not to merge the branch, in my opinion. On the other hand, I would very much welcome input to make things work faster for complicated expressions like those described above.

Barring any good ideas, I might eventually end up adding an optional argument exact with a default value of True. Setting that to False would skip this check, and therefore possibly return an approximate symbolic ring element.

[nbruin]

comment:12

The problem is, an expression like this:

(2/9*I*sqrt(443)*sqrt(3) + 53/27)^(1/3)

has (at least) 3 possible values. You need to specify branch cuts to single out a particular value in, say, CC. Algebraically, sqrt(443), sqrt(3), and I have the same problem. So when you first construct this element you are making essentially the ring

QQ[x,y,z,w]/(x^2-443,y^2-3,z^2+1,w^3-(2/9*z*x*y + 53/27))

which is of rather high degree.

The other elements you add to this will just add to the degree, since the next time a sqrt(443) is encountered, it is not algebraically clear that this is supposed to be the same as the one designated by x. The complex embedding that comes with AA or Qbar specifies some of the ambiguity (e.g., the branch cut of sqrt(443) will ensure the "positive" root always), but once you start taking roots of complex numbers, the cuts might not even do what they algebraically are supposed to (e.g., -(z)^(1/3) != (-z)^(1/3) ).

While in the expression for b that you give, it is exactly the same expression that you are taking the cube root of repeatedly, one can only see this after careful inspection. For a human it's clear you're meaning the "same" cube root every time, but the computer has no immediate way of finding that out. A priori, all the choices made by sqrt... and (...)^(1/3) are considered as algebraically unrelated and only with laborious computations using the embedding in CC does one rediscover those relations. If instead you do

u=Qbar(sqrt(-3))
v=Qbar(sqrt(443))*u
w=Qbar( (2/9*v+53/27)^(1/3))

and then construct b from that, you should already get a little faster results.

The problems you are running into are fairly well-understood and real. Symbolic notation as you use it fails to express the full algebraic relations once you start repeating related expressions.

EDIT: apologies. The fact that the value didn't match was due to a missubstitution on my part. I've corrected the result below. The value b computed does match what it's supposed to be:

sage: u=QQbar(sqrt(-3))
sage: v=QQbar(sqrt(443))*u
sage: w=(2/9*v + 53/27)^(1/3)
sage: b = -1/1296*((9*u*w + 9*w - 37*u/w + 37/w + 18)*(9*u*w + 9*w - 37*u/w + 37/w - 12) - 3456)*sqrt(-9/2*u*w - 9/2*w + 37/2*u/w - 37/2/w + 6)
sage: b
4.904821984561304? + 0.?e-16*I
sage: 16*b.minpoly() #still takes a long time
16*x^6 - 392*x^4 + 133*x^2 + 900

The fundamental problem remains, though: the symbolic version doesn't quite carry the same information. In addition, some of the slowness might be due to ``QQbar` itself, since:

sage: t=(-9/2*u*w - 9/2*w + 37/2*u/w - 37/2/w + 6)
sage: t.minpoly() # this is fast
x^3 - 18*x^2 - 891*x + 2916
sage: r=sqrt(t)
sage: r.minpoly() # this is quite slow
x^6 - 18*x^4 - 891*x^2 + 2916
sage: b = -1/1296*((9*u*w + 9*w - 37*u/w + 37/w + 18)*(9*u*w + 9*w - 37*u/w + 37/w - 12) - 3456)*r
sage: b.minpoly() #this is now pretty quick
x^6 - 49/2*x^4 + 133/16*x^2 + 225/4

so QQbar is slow exactifying an element r for which it has all relevant algebraic information (it's the square root of t, which has a known minimal polynomial)

Computer algebra systems such as maple have their own symbolic ways of dealing with the ambiguity introduced by defining elements as "roots" of polynomials:

> lprint(solve(x^5+4*x+1));
RootOf(_Z^5+4*_Z+1,index = 1), RootOf(_Z^5+4*_Z+1,index = 2), RootOf(_Z^5+4*_Z
+1,index = 3), RootOf(_Z^5+4*_Z+1,index = 4), RootOf(_Z^5+4*_Z+1,index = 5)

which allows it to simplify RootOf(_Z^5+4*_Z+1,index = 1)-RootOf(_Z^5+4*_Z+1,index = 1) to 0 but leave RootOf(_Z^5+4*_Z+1,index = 1)-RootOf(_Z^5+4*_Z+1,index = 2) alone. This doesn't capture the algebraic dependencies between the roots, but at least allows identification of each root individually.

[gagern]

comment:13

Thanks for explaining a bit of what's going on here to cause such delays. Alternate ways to compute things more efficiently are all very well, but won't be of use unless I can make this automatic somehow. After all, on the high level all I do is compare two numbers, without wanting to think about implementation details.

While investigating this problem, I noticed that for my original length expression b the call b.minpoly() works really fast, whereas QQbar(b).minpoly() takes ages. Makes me wonder whether the exactification should be somehow based on the minimal polynomial of the symbolic expression. Or at least leverage some of the machinery which makes that so fast. That should speed up the kind of comparison I'm using. Or is that computation somehow inexact?

Let me try to think this through. Naively I'd rewrite sage.symbolic.expression_conversions.algebraic to first compute the minpoly of the expression in question. I fear that this might take long, which is a problem in those cases where we don't need the exact representation. So it would probably be better to only do this on demand. Which raises several related questions.

• Can we store a reference to a symbolic expression without having to worry that it will get modified? Are symbolic expressions read-only after construction?
• Can we afford the cost of storing whole symbolic expressions in addition to their AA or QQbar converted forms? Or should we try to reconstruct the symbolic expression from the descriptor DAG when needed?
• How do we link the symbolic expression with the algebraic number? We could invent a descriptor which is layered around the one returned from the current conversion approach, but whose exactify method will compute the minimal polynomial of the whole expression instead of recursing.
• On the other hand, we might need access to the current interval, so an alternative would be storing such information in AlgebraicNumber_base and modify its exactify method to take the additional information into account.
• Are there cases where a symbolic minpoly would take significantly more time than the corresponding recursive exactify? If so, can we somehow detect these cases before we start the symbolic minpoly computation?
• We need to find roots for that polynomial, which I'd so via a simple p.roots(self.parent(), False) unless someone suggests otherwise.
• Wen need to compare them to the approximation returned by the nested descriptor. Which is probably the main reason to have this in AlgebraicNumber_base instead of some ANDescr descendant.
• We might have to increase precision until the interval contains no more than a single solution. Simple with AlgebraicNumber_base._more_precision but difficult or wasteful for ANDescr.

Does this approach make sense?

[sagetrac-git]

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

cb09ec6

Fix formatting of doctest block.

[sagetrac-git]

Changed commit from 8486d6c to cb09ec6

[sagetrac-git]

Changed commit from cb09ec6 to 4ccb707

[sagetrac-git]

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

472c754	Implement symbolic (radical) roots of polynomials.
989674d	Merge branch ticket/17517 into ticket/14239.
4ccb707	Delegate computation of symbolic roots to the roots method of the polynomial.

[sagetrac-git]

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

890cf9b	Increase precision when converting to radical expression.
ee3823a	Avoid code duplication between qqbar and nfe.
52512cc	Simplify NumberFieldElement._symbolic_
229b1fb	Use refine_embedding to convert to AA/QQbar
981942b	Add an example in degree 10 with degree 2 subfield
3423450	Fall back to convertig generator to symbolic.
73df99f	Drop __symbolic member from NumberFieldElement.
da4aa6e	Indicate temporary nature of generator-based conversion.
a5a9245	Fix formatting of doctest block.
4626286	Delegate computation of symbolic roots to the roots method of the polynomial.

[sagetrac-git]

Changed commit from 4ccb707 to 4626286

[gagern]

comment:97

Rebased since we had a conflict with e216b6f from #16908.

Am I right that we are just waiting for Jeroen or someone else to give this a positive review? Or is there still something missing which I should address in order to let this proceed?

[rwst]

comment:98

Passes make ptestlong. Can't comment on the code, sorry.

[videlec]

comment:99

What is the point of

-        Test :trac:`14895`::
-
-            sage: K.<sqrt2> = QuadraticField(2)
+            sage: K.<sqrt2> = QuadraticField(2) # :trac:`14895`

I really do not like the attribute symbolic. Why this result needs to be cached? You can possibly cache in the parent a symbolic version of the generator.

The fact that SR(a) when a is a number field element returns possibly a numerical approximation is a very bad regression. You switch from something exact to something approximative.

[videlec]

comment:100

What is the point of looking for an approximation

k = ( K._n()*CC(K.gen()).log() / CC(two_pi_i) ).real().round() # n ln z / (2 pi i)

[videlec]

comment:101

Instead of embedding.im_gens()[0] you can use the clearer embedding.gen_image().

[gagern]

Changed author from Martin von Gagern to Martin von Gagern, Jeroen Demeyer

[gagern]

comment:102

Replying to @videlec:

What is the point of

-        Test :trac:`14895`::
-
-            sage: K.<sqrt2> = QuadraticField(2)
+            sage: K.<sqrt2> = QuadraticField(2) # :trac:`14895`

That one is from 4626286 by Jeroen, but you have the patch reversed. I assume that :trac:`…` work in plain text but not in doctest comments, so this made a lot of sense to me.

I really do not like the attribute symbolic. Why this result needs to be cached? You can possibly cache in the parent a symbolic version of the generator.

Could it be you are seeing the wrong patch direction here as well? We are dropping the __symbolic member and its caching functionality, since we don't expect to need it except in rare cases as a temporary solution.

The fact that SR(a) when a is a number field element returns possibly a numerical approximation is a very bad regression. You switch from something exact to something approximative.

Yes, you are definitely seeing the patch in the wrong direction. A good thing that you pasted the patch lines above, otherwise I would have ben thoroughly confused.

---

Replying to @videlec:

What is the point of looking for an approximation

k = ( K._n()*CC(K.gen()).log() / CC(two_pi_i) ).real().round() # n ln z / (2 pi i)

I must confess that I was a bit unhappy with that myself when I read it. But the code just changed place and location, it wasn't really changed. And I figure that if you ever have a cyclotomic field of such a great order that you are affected by rounding errors here, its dimension would be high enough that you'd not get any real work done with it in any case. So I decided not to complain or modify this, since it's not directly related to the issue at hand and I feard that discussing one more item might prevent stuff from getting accepted.

Are you sure you want this handled with this ticket here? If so I'll have to work out how to find the right k using exact arithmetic. I guess the integers are all there somewhere, but I'll have to look closer. It's been a while.

---

Replying to @videlec:

Instead of embedding.im_gens()[0] you can use the clearer embedding.gen_image().

Thanks, that should be a quick and obvious improvement. Will push a commit to that effect once I know whether to address the approximation issue above as well.

[videlec]

comment:103

Replying to @gagern:

Yes, you are definitely seeing the patch in the wrong direction. A good thing that you pasted the patch lines above, otherwise I would have ben thoroughly confused.

Oups... sorry

Replying to @videlec:

What is the point of looking for an approximation

k = ( K._n()*CC(K.gen()).log() / CC(two_pi_i) ).real().round() # n ln z / (2 pi i)

Are you sure you want this handled with this ticket here? If so I'll have to work out how to find the right k using exact arithmetic. I guess the integers are all there somewhere, but I'll have to look closer. It's been a while.

It should be easy, NumberField_cyclotomic does have a method .zeta_order(). But then, the exact root of unity depends on the embedding... and I do not see how to avoid going to QQbar (or using in an ugly way interval approximations). It is really a pity that we use embedding into RLF/CLF by default and not AA/QQbar (see also #18103, #18104).

Vincent

[gagern]

comment:104

Replying to @videlec:

It should be easy, NumberField_cyclotomic does have a method .zeta_order(). But then, the exact root of unity depends on the embedding... and I do not see how to avoid going to QQbar (or using in an ugly way interval approximations).

It is really a pity that we use embedding into RLF/CLF by default and not AA/QQbar (see also #18103, #18104).

Yes, using qqbar for the generator would certainly make things better. How about we leave the approximation in place for now, and file a separate ticket to get rid of that once the generator comes from qqbar?

[gagern]

comment:105

Replying to @videlec:

Instead of embedding.im_gens()[0] you can use the clearer embedding.gen_image().

That won't work: at that point (since the refine_embedding line a bit above that), the object in question is a morphism.NumberFieldHomomorphism_im_gens and not a number_field_morphisms.NumberFieldEmbedding any more, so it doesn't have a gen_image method. I have no idea whether that's how it should be, but I'm fairly certain that if this should change, some other ticket should take care of that.

So I'm back to asking for review, even though I changed nothing.

[videlec]

comment:106

Replying to @gagern:

Replying to @videlec:

Instead of embedding.im_gens()[0] you can use the clearer embedding.gen_image().

That won't work: at that point (since the refine_embedding line a bit above that), the object in question is a morphism.NumberFieldHomomorphism_im_gens and not a number_field_morphisms.NumberFieldEmbedding any more, so it doesn't have a gen_image method. I have no idea whether that's how it should be, but I'm fairly certain that if this should change, some other ticket should take care of that.

So I'm back to asking for review, even though I changed nothing.

That's good to me!

Vincent

[videlec]

Changed reviewer from Marc Mezzarobba, Jeroen Demeyer to Marc Mezzarobba, Jeroen Demeyer, Vincent Delecroix

[vbraun]

Changed branch from u/gagern/ticket/14239 to 4626286