Open 8d15854a-f726-4f6b-88e7-82ec1970fbba opened 10 years ago
Branch: u/gagern/ticket/16222
These modifications are not ready for inclusion yet, they need some discussion first.
How to get the symbolic expression from expression_conversion
into AlgebraicNumber_base
.
Right now, I use direct access to an underscore-prefixed variable called _symbolic_expression
. I'm not sure whether that is acceptable. If cross-module interfaces should avoid using underscores, we'll have to look for alternatives. But I rely on the fact that the expression really describes the number in question, so this certainly isn't a variable the user should be allowed to tweak. So I don't know a better way yet.
Default state of that variable.
Right now, I have a class variable which defaults to None
and instance variables which will override this when required. Alternatives would be setting an instance variable in some (as of now non-existant) constructor, or avoding a default value altogether and using hasattr
instead. I simply wrote what I'd do in my own code, but if there is a sage convention to write things differently, I'll gladly comply.
When to compute the minpoly.
I now do the computation on demand, hoping that expressions cannot be legally modified after their creation. If they can be modified, we might want to do the minpoly computation when we do the conversion from symbolic to algebraic. Likewise if you have reasons to believe that minpoly computation should not make that conversion noticably slower than it already is.
Test case dependency.
I just noticed that my current doctest depends on ticket:14239, but that should be easy to fix. If you have a favorite example, let me know, but otherwise I'll come up with something.
The _more_precision
loop.
I have more trouble testing that _more_precision
codepath inside _exact_symbolic
. Can anyone come up with a situation where identifying the correct root would need more than default precision but where numeric minpoly computation would still succeed? If not, what shall we do? Shall we deactivate that whole loop, and simply return without taking action if we find more than one root which could match?
How to obtain the roots.
Is p.roots(self.parent(), False)
the right thing to do, or is this too high-level, should I use something more basic here?
New commits:
bb01e0b | Drive-by fix to minpoly documentation. |
491dcc4 | Faster exactification using numeric minpoly. |
Commit: 491dcc4
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.
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
Replying to @videlec:
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 non-symbolic case, or find a rigorous argument why we don't need such a proof?
Replying to @gagern:
Replying to @videlec:
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 non-symbolic 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()
.
Replying to @videlec:
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?
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))
Replying to @videlec:
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?
Replying to @gagern:
Replying to @videlec:
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.
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.
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 non-example 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.
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^16-64*x^14+1648*x^12-23552*x^10+213480*x^8-1191424*x^6+3340480*x^4-3562496*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
Replying to @videlec:
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 spin-off 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.
Replying to @videlec:
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 follow-up ticket to investigate that.
Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
4b1d01e | Drive-by 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” |
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.
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 n-th 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 |
Changed branch from u/gagern/ticket/16222 to u/vdelecroix/16222
Reviewer: Vincent Delecroix
Are there any link/dependency with #14239 or/and #17516?
Replying to @videlec:
Are there any link/dependency with #14239 or/and #17516?
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?
Replying to @gagern:
… 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?
Dependencies: #15605
I fear we'll have to wait for #15605, the ticket #18129 got duplicated to, before we can trust SR to get this right…
This is a spinoff from #14239 comment:13. There I noticed that for a large symbolic expressions
b
, the callb.minpoly()
was a lot (by several orders of magnitude) faster than the callQQbar(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.Depends on #15605
Component: number fields
Author: Martin von Gagern
Branch/Commit: u/vdelecroix/16222 @
8773468
Reviewer: Vincent Delecroix
Issue created by migration from https://trac.sagemath.org/ticket/16222