Opened 6 years ago
Last modified 5 years ago
#16222 needs_info enhancement
Faster exactification using numeric minpoly
Reported by:  gagern  Owned by:  

Priority:  major  Milestone:  sage6.4 
Component:  number fields  Keywords:  
Cc:  Merged in:  
Authors:  Martin von Gagern  Reviewers:  Vincent Delecroix 
Report Upstream:  N/A  Work issues:  
Branch:  u/vdelecroix/16222 (Commits)  Commit:  8773468739544ed9f12b836c7f28a06be37ffd94 
Dependencies:  #15605  Stopgaps: 
Description
This is a spinoff from comment:13:ticket:14239. There I noticed that for a large symbolic expressions b
, the call b.minpoly()
was a lot (by several orders of magnitude) faster than the call QQbar(b).minpoly()
. We should try to make this speed gain available to those algebraic numbers which were constructed from symbolic expressions.
I now know that the speed gain is almost certainly due to the numeric algorithm in calculus.minpoly
. So what we should in my opinion do is use that numeric algorithm to obtain a minimal polynomial of the whole expression, and if that succeeds, base subsequent exactifications on that polynomial instead of the nested tree of algebraic number descriptors.
Change History (26)
comment:1 Changed 6 years ago by
 Branch set to u/gagern/ticket/16222
 Created changed from 04/23/14 19:33:11 to 04/23/14 19:33:11
 Modified changed from 04/23/14 19:33:11 to 04/23/14 19:33:11
comment:2 Changed 6 years ago by
 Commit set to 491dcc4884823cab627f318ea63f0e2504f9f7c2
comment:3 Changed 5 years ago by
 Commit changed from 491dcc4884823cab627f318ea63f0e2504f9f7c2 to 8597313c855ce7a413cfcabead498059d2cfcbf6
comment:4 Changed 5 years ago by
 Status changed from new to needs_review
Rebased to current develop, and changed test case to avoid dependency on #14239. I'm now requesting review since that's what I need even if I consider it likely that the reviewer will not agree with all the choices I made as outlined above.
comment:5 Changed 5 years ago by
 Milestone changed from sage6.3 to sage6.4
comment:6 followup: ↓ 7 Changed 5 years ago by
Hello,
Did you try using the magic algdep
from pari? It perform (very quickly) some LLL to find the polynomial with smallest coefficient for which x
is almost a root. With your example
sage: a = (443/96*I*sqrt(443)*sqrt(3) + 833939/1728)^(1/3) sage: b = sqrt(144*a + 9205/a + 1176)/12 sage: c = QQbar(b) sage: gp.algdep(c.n(prec=150).real(), 6) 16*x^6  392*x^4 + 133*x^2 + 900
The advantage is that it would potentially speed up non symbolic cases.
Vincent
comment:7 in reply to: ↑ 6 ; followup: ↓ 8 Changed 5 years ago by
Replying to vdelecroix:
Did you try using the magic
algdep
from pari?
My code builds on that, since expression.minpoly calls calculus.minpoly which in turn calls algdep. But it does come with a verification step, ensuring that the found polynomial is in fact the one we need. I'm not sure I'd trust the output from algdep
directly, particularly not unless I had some very good ideas of what degree to expect.
The advantage is that it would potentially speed up non symbolic cases.
Well, one alternative would be to directly feed the descriptor DAG into calculus.minpoly
. Computing a numerical_approx
from these would be easy enough. The hard part would be the symbolic verification step, g(ex).simplify_trig().canonicalize_radical() == 0
. I guess that's done in Maxima. Do you have an idea how we might either implement this for the nonsymbolic case, or find a rigorous argument why we don't need such a proof?
comment:8 in reply to: ↑ 7 ; followup: ↓ 9 Changed 5 years ago by
Replying to gagern:
Replying to vdelecroix:
Did you try using the magic
algdep
from pari?My code builds on that, since expression.minpoly calls calculus.minpoly which in turn calls algdep. But it does come with a verification step, ensuring that the found polynomial is in fact the one we need. I'm not sure I'd trust the output from
algdep
directly, particularly not unless I had some very good ideas of what degree to expect.
Of course we can not believe the output of algdep
. But at least we can bound the degree. And once a polynomial is given we can check irreducibility. So I guess we can produce a candidate quite reliably.
The advantage is that it would potentially speed up non symbolic cases.
Well, one alternative would be to directly feed the descriptor DAG into
calculus.minpoly
. Computing anumerical_approx
from these would be easy enough. The hard part would be the symbolic verification step,g(ex).simplify_trig().canonicalize_radical() == 0
. I guess that's done in Maxima. Do you have an idea how we might either implement this for the nonsymbolic case, or find a rigorous argument why we don't need such a proof?
Right, checking that it is actually an annihilator polynomial is harder... I do not see anything different from calling __nonzero__
which indeed calls .exactify()
.
comment:9 in reply to: ↑ 8 Changed 5 years ago by
Replying to vdelecroix:
So I guess we can produce a candidate quite reliably.
I'm unsure about the relation between “candidate” and “reliable”. The former implies something that needs to be verified, the latter like I could rely on it without verification. But I guess what you mean is that we could construct a candidate which would very likely be the correct one. But “very likely” isn't enough in my opinion.
Right, checking that it is actually an annihilator polynomial is harder... I do not see anything different from calling
__nonzero__
which indeed calls.exactify()
.
In theory, we could go the reverse route: take the descriptor dag and try building a symbolic expression from this which in turn can be fed to Maxima for verification. But I'm far from certain that this is a good idea for those cases which didn't originate in a symbolic expression in the first place. Which brings us back to the request as it currently stands, right?
comment:10 followup: ↓ 11 Changed 5 years ago by
I just notice that your ticket is based on a very old version of Sage (6.3.beta5)... it would be nice to upgrade to the current 6.5.rc0.
More about the ticket, I am worried because self._symbolic_expression
is completely ignored in the following cases
sage: a = QQbar(sqrt(2)) + 1 sage: b = QQbar(sqrt(2)) + QQbar(sqrt(3))
comment:11 in reply to: ↑ 10 ; followup: ↓ 12 Changed 5 years ago by
Replying to vdelecroix:
More about the ticket, I am worried because
self._symbolic_expression
is completely ignored in the following cases […]
So what do you propose? That we always try to build a symbolic expression from the descriptor dag?
comment:12 in reply to: ↑ 11 Changed 5 years ago by
Replying to gagern:
Replying to vdelecroix:
More about the ticket, I am worried because
self._symbolic_expression
is completely ignored in the following cases […]So what do you propose? That we always try to build a symbolic expression from the descriptor dag?
No. But at least when performing arithmetic operations on QQbar
or AA
elements the operation are also done on the corresponding symbolic expressions.
comment:13 Changed 5 years ago by
 Commit changed from 8597313c855ce7a413cfcabead498059d2cfcbf6 to 998534a33a502a86518f6c541f99cac1ffacdaa5
comment:14 followup: ↓ 15 Changed 5 years ago by
 Status changed from needs_review to needs_work
I find that it's hard to draw the line. Your example of QQbar(sqrt(2)) + 1
shows that we want symbolic expressions for minpoly construction even if some arguments were not explicitely constructed from symbolic expressions. But where does it stop? In my demo implementation, I decided that caching the symbolic expression after every elementary operation would duplicate a lot of information. So I decided to recursively construct such symbolic expressions. The result is pretty much a construction of symbolic expressions from the descriptor dag, even though I didn't think of it that way when I started coding. It can generate symbolic expression for a lot of things.
And I'm using it unconditionally, which leads to a lot of doctests failing. Some of them for better:
line 3653, rt3b.as_number_field_element(): Expected: … defining polynomial y^4  4*y^2 + 1, a^2 + 2 … Got: … defining polynomial y^2  3, a … line 3755, rt2b._exact_value(): Expected: a^3  3*a where a^4  4*a^2 + 1 = 0 and a in 1.931851652578137? Got: a where a^2  2 = 0 and a in 1.414213562373095? line 7171, y.is_complex(): Expected: True Got: False
But some of them for worse:
line 360, sage_input(n, verify=True): Expected: AA(2) Got: v = sqrt(AA(2)) v*v line 7881, sage_input(n): Expected: QQbar(sqrt(AA(2))) + 1 Got: v = sqrt(QQbar(3)) QQbar(sqrt(AA(2))) + v/v
I don't have any performance comparisons yet, but I guess that there might well be cases where the current approach was fast enough, but the whole detour via symbolic expressions, PARI's algdep and Maxima for verification might have a severe impact.
So my key concern right now is, when do we want to use this approach, and when not to? If multithreading in Python were a more natural concept, I'd be so bold as to suggest that we start both approaches, and the one which terminates first wins, aborting the other. This might make results somewhat unpredictable, though. And given the state of threading support in CPython, I doubt this is realistic in any case.
Here are a few more realistic ways how we could make that case distinctions in a way that might make sense. I'm not completely happy with any of them.
 Only do the algdep approach if at least one number in the descriptor dag originated from a symbolic expression
 Limit the reconstruction of symbolic expressions to fewer operations, mainly elementary arithmetic operations
 Give control to the user: make exactification using algdep not part of the normal exactification procedure, but instead the result of an explicit method invocation
 Try to detect exactifications which are likely very slow using the standard approach, e.g. indicated by the expected size of the minpoly for some union field, and try algdep only on those
comment:15 in reply to: ↑ 14 ; followup: ↓ 18 Changed 5 years ago by
Replying to gagern:
But some of them for worse:
line 360, sage_input(n, verify=True): Expected: AA(2) Got: v = sqrt(AA(2)) v*v line 7881, sage_input(n): Expected: QQbar(sqrt(AA(2))) + 1 Got: v = sqrt(QQbar(3)) QQbar(sqrt(AA(2))) + v/v
Not really for worse. The previous doctests assumed that at some point in the code .exactify()
would have been called for all nodes of the description tree... Right now, the symbolic code ignore the exactification of the nodes.
EDIT: I removed the nonexample I previously wrote here.
Here are a few more realistic ways how we could make that case distinctions in a way that might make sense. I'm not completely happy with any of them.
 Only do the algdep approach if at least one number in the descriptor dag originated from a symbolic expression
 Limit the reconstruction of symbolic expressions to fewer operations, mainly elementary arithmetic operations
 Give control to the user: make exactification using algdep not part of the normal exactification procedure, but instead the result of an explicit method invocation
 Try to detect exactifications which are likely very slow using the standard approach, e.g. indicated by the expected size of the minpoly for some union field, and try algdep only on those
Option 3 is definitely needed.
I am against option 1 as we do not want to keep the origin of the objects. By the way, do you have an idea to have quick symbolic representation for
sage: QQbar(sqrt2).sqrt() + QQbar(3).sqrt()
I.e. roots of x^n  a
have a simple symbolic counterpart.
For options 2 and 4 there is a right balance to find.
comment:16 followup: ↓ 17 Changed 5 years ago by
Not exactly related to this ticket: I just found out some maxima magic to get the polynomial of symbolic expressions involving square root and such (p. 1041 of the maxima manual):
(%i1) load(to_poly_solve)$ (%i2) first(elim_allbut(first(to_poly(eq,[1,x])),[x])); 4 2 (%o2) [x  24 x + 4] (%i3) eq: x = sqrt(sqrt(3) + sqrt(5) + 1) + sqrt(7); (%o3) x = sqrt(7) + sqrt(sqrt(5) + sqrt(3) + 1) (%i3) first(elim_allbut(first(to_poly(eq, [1,x])), [x])); 16 14 12 10 8 6 (%o3) [x  64 x + 1648 x  23552 x + 213480 x  1191424 x 4 2 + 3340480 x  3562496 x + 524176]
It is way much much faster than what we have in Sage for minpoly:
sage: maxima.load("to_poly_solve") "/opt/sage_flatsurf/local/share/maxima/5.35.1/share/to_poly_solve/to_poly_solve.mac" sage: maxima("eq: x = sqrt(sqrt(3) + sqrt(5) + 1) + sqrt(7);") x=sqrt(7)+sqrt(sqrt(5)+sqrt(3)+1) sage: maxima("first(elim_allbut(first(to_poly(eq, [1,x])), [x]));") [x^1664*x^14+1648*x^1223552*x^10+213480*x^81191424*x^6+3340480*x^43562496*x^2+524176] sage: %timeit maxima("first(elim_allbut(first(to_poly(eq, [1,x])), [x]));") 10 loops, best of 3: 18 ms per loop sage: a = sqrt(7) + sqrt(sqrt(5) + sqrt(3) + 1) sage: %timeit a.minpoly() 1 loops, best of 3: 449 ms per loop
comment:17 in reply to: ↑ 16 Changed 5 years ago by
Replying to vdelecroix:
Not exactly related to this ticket: I just found out some maxima magic to get the polynomial of symbolic expressions involving square root and such
Nice! Do you want to open a spinoff ticket, asking for a way to exploit this functionality in Sage? As you said, it's not exactly on topic here, but it would be a shame if that idea got lost once we have this issue here wrapped up.
comment:18 in reply to: ↑ 15 Changed 5 years ago by
Replying to vdelecroix:
Option 3 is definitely needed.
For options 2 and 4 there is a right balance to find.
I guess once we have a method which can be called by hand, we can use that to experiment a bit, in our daily work, and can ask others to do the same. That might give us some intuition about when to use which. So I'd keep this out of this ticket for now, and instead post a followup ticket to investigate that.
comment:19 Changed 5 years ago by
 Commit changed from 998534a33a502a86518f6c541f99cac1ffacdaa5 to 96942f09ea8a9cddf2b2dc8fece352935ddfe3b2
Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
4b1d01e  Driveby fix to minpoly documentation.

f065ca4  Faster exactification using numeric minpoly.

8cc7bc7  Build symbolic expression from descriptor graph of algebraic number.

96942f0  Trac #16222: Add argument “symbolic” to method “exactify”

comment:20 Changed 5 years ago by
 Status changed from needs_work to needs_review
Rebased to 6.6.rc2 and added public interface z.exactify(symbolic=True)
.
I considered “numeric” instead of “symbolic”, but I at least tend to associate “numeric” with “inexact”, so I felt that term might be more confusing. Also considered “algoiritm='symbolic'” but I think in the long run we might have various alternatives which could be tried one after the other, so a flag to enable or disable each sounds more useful. The same holds for a new method as opposed to an argument to exactify
. So that's why I chose the interface the way I did.
Now that my new code is not called by default, it shouldn't have any adverse effects on other parts of the code, which might make this easier to review. I hope this can be merged, to enable some evaluation in preparation for #18122 which I filed for automatically choosing between algorithms. I fear that might be hard to get right.
comment:21 Changed 5 years ago by
 Branch changed from u/gagern/ticket/16222 to u/vdelecroix/16222
 Commit changed from 96942f09ea8a9cddf2b2dc8fece352935ddfe3b2 to 8773468739544ed9f12b836c7f28a06be37ffd94
Hello,
It's a pity that we got None
in
sage: AA(2).sqrt()._descr.quick_symbolic()
I added a new commit to handle that.
On the other hand, did you know that the symbolic ring is very broken to deal with nth root
sage: a = (1)^(1/3) sage: a.n() 0.500000000000000 + 0.866025403784439*I sage: a**2 1
As bool(I^(2/3) == (1)^(1/3))
is True
we might get into troubles with subtle examples. I am very worried by that.
Vincent
New commits:
8773468  Trac 16222: support quick symbolic for roots of x^n a

comment:22 Changed 5 years ago by
 Reviewers set to Vincent Delecroix
 Status changed from needs_review to needs_info
comment:23 followup: ↓ 24 Changed 5 years ago by
comment:24 in reply to: ↑ 23 ; followup: ↓ 25 Changed 5 years ago by
Replying to vdelecroix:
No, since in general they are not what I'd consider fast. To obtain an exact symbolic number using #14239, we need to have a minimal polynomial, while this one here is aiming to obtain that polynomial more quickly. #17516 would extend #14239 to give results in a larger number of cases, but won't affect this here either. If at all, there is some overlap with #17886, since if that works out relly well, it might make this one here superfluous. I don't think so, though, and since we are only offering this here as an alternative for now, we can use #18122 to evaluate those alternatives once we get there.
… did you know that the symbolic ring is very broken …
No, I did not, but I share your deep concerns. Have you filed a ticket for this yet?
comment:25 in reply to: ↑ 24 Changed 5 years ago by
comment:26 Changed 5 years ago by
 Dependencies set to #15605
These modifications are not ready for inclusion yet, they need some discussion first.
expression_conversion
intoAlgebraicNumber_base
._more_precision
loop.New commits:
Driveby fix to minpoly documentation.
Faster exactification using numeric minpoly.