Sage: Ticket #11941: Solve and assumptions too aggressive with cube root of negative numbers
https://trac.sagemath.org/ticket/11941
<p>
<a class="closed ticket" href="https://trac.sagemath.org/ticket/6515" title="defect: assume doesn't interact well with solve (closed: fixed)">#6515</a> did a great job helping us start to catch some assumptions when we do solving.
</p>
<p>
However, <a class="ext-link" href="http://ask.sagemath.org/question/824/real-solution-of-x38-0"><span class="icon"></span>this ask.sagemath.org post</a> catches a case where it's too aggressive, because Sage says that <code>(-1)^(1/3)</code> is not real.
</p>
<pre class="wiki">sage: solve(x^3+1==0,x)
[x == 1/2*I*(-1)^(1/3)*sqrt(3) - 1/2*(-1)^(1/3), x == -1/2*I*(-1)^(1/3)*sqrt(3) - 1/2*(-1)^(1/3), x == (-1)^(1/3)]
sage: assume(x,'real')
sage: solve(x^3+1==0,x)
[]
</pre><p>
What's weird about this is that the Maxima in Sage should just return <code>x==-1</code>.
</p>
<pre class="wiki">(%i2) display2d:false;
(%o2) false
(%i3) solve(x^3+1=0,x);
(%o3) [x = -(sqrt(3)*%i-1)/2,x = (sqrt(3)*%i+1)/2,x = -1]
</pre><p>
Not sure what's going on with that.
</p>
en-usSagehttps://trac.sagemath.org/chrome/site/logo_sagemath_trac.png
https://trac.sagemath.org/ticket/11941
Trac 1.1.6jdemeyerThu, 03 Nov 2011 16:14:43 GMTmilestone deleted
https://trac.sagemath.org/ticket/11941#comment:1
https://trac.sagemath.org/ticket/11941#comment:1
<ul>
<li><strong>milestone</strong>
<em>sage-4.7.3</em> deleted
</li>
</ul>
<p>
Milestone sage-4.7.3 deleted
</p>
TicketmboratkoThu, 17 Nov 2011 06:08:04 GMTmilestone set
https://trac.sagemath.org/ticket/11941#comment:2
https://trac.sagemath.org/ticket/11941#comment:2
<ul>
<li><strong>milestone</strong>
set to <em>sage-4.8</em>
</li>
</ul>
<p>
The fact that Maxima normally returns -1 and Sage returns (-1)<sup>(1/3) is a bit odd, as you mentioned. At a more basic level, Sage doesn't seem to think that (-1)</sup>(1/3) is in RR:
</p>
<pre class="wiki">sage: (-1)^(1/3) in RR
False
sage: (2)^(1/3) in RR
True
</pre><p>
So if we fix that problem, then at least it would return (-1)<sup>(1/3). I also suspect that it would properly simplify to -1 at that point as well, based on the following example:
</sup></p>
<pre class="wiki">sage: solve(x^3-8==0,x)
[x == I*sqrt(3) - 1, x == -I*sqrt(3) - 1, x == 2]
sage: solve(x^3+8==0,x)
[x == I*(-1)^(1/3)*sqrt(3) - (-1)^(1/3), x == -I*(-1)^(1/3)*sqrt(3) - (-1)^(1/3), x == 2*(-1)^(1/3)]
</pre><p>
(Note that the (1/3) exponent appears everywhere next to the -1, as if some rule specifies that sage should not simplify it out.)
</p>
<p>
I am new, however, and I am not sure where next to look.
</p>
TicketkcrismanThu, 17 Nov 2011 14:18:27 GMT
https://trac.sagemath.org/ticket/11941#comment:3
https://trac.sagemath.org/ticket/11941#comment:3
<p>
Well, in general we do <em>not</em> want to do this. It's been discussed ad nauseam <a class="ext-link" href="http://groups.google.com/group/sage-support/browse_thread/thread/66d45649ca8c1cc3?pli=1"><span class="icon"></span>many</a> <a class="ext-link" href="http://ask.sagemath.org/question/260/fractional-power-to-negative-number"><span class="icon"></span>times</a>, and the sense is that:
</p>
<ul><li>Similar programs don't necessarily do this
</li><li><code>(-1)^(1/3)</code> is not really <code>-1</code> but a primitive complex root of <code>-1</code>.
</li></ul><p>
This ticket is about the fact that Maxima returns three solutions to the equation, but when we do the <code>assume(x,'real')</code> they <em>all</em> mysteriously vanish!
</p>
TicketburcinThu, 17 Nov 2011 16:05:10 GMT
https://trac.sagemath.org/ticket/11941#comment:4
https://trac.sagemath.org/ticket/11941#comment:4
<p>
Here's what maple does:
</p>
<pre class="wiki">> solve( x^3+1= 0, x);
1/2 1/2
-1, 1/2 - 1/2 I 3 , 1/2 + 1/2 I 3
</pre><p>
I wonder why maxima is returning <code>(-1)^(1/3)</code>. Maybe we should ask the Maxima developers.
</p>
TicketkcrismanThu, 17 Nov 2011 16:12:25 GMT
https://trac.sagemath.org/ticket/11941#comment:5
https://trac.sagemath.org/ticket/11941#comment:5
<blockquote class="citation">
<p>
I wonder why maxima is returning <code>(-1)^(1/3)</code>. Maybe we should ask the Maxima developers.
</p>
</blockquote>
<p>
No, no! See this example from the description:
</p>
<pre class="wiki">(%i2) display2d:false;
(%o2) false
(%i3) solve(x^3+1=0,x);
(%o3) [x = -(sqrt(3)*%i-1)/2,x = (sqrt(3)*%i+1)/2,x = -1]
</pre><p>
We are somehow getting the <code>(-1)^(1/3)</code> by doing something nonstandard in Maxima, apparently. But their vanilla thing is just right.
</p>
TicketmboratkoThu, 17 Nov 2011 16:15:08 GMT
https://trac.sagemath.org/ticket/11941#comment:6
https://trac.sagemath.org/ticket/11941#comment:6
<p>
It seems that sage sets domain: complex (I was made aware of this by burcin in IRC). You do get this result as follows:
</p>
<p>
{{{(%i8) domain:complex;
</p>
<p>
(%o8) complex
(%i9) solve(x<sup>3+1=0,x);
</sup></p>
<p>
(%o9) [x = ((-1)<sup>(1/3)*sqrt(3)*%i-(-1)</sup>(1/3))/2,
</p>
<blockquote>
<p>
x = -((-1)<sup>(1/3)*sqrt(3)*%i+(-1)</sup>(1/3))/2,x = (-1)<sup>(1/3)]
</sup></p>
</blockquote>
<p>
}}}
</p>
TicketmboratkoThu, 17 Nov 2011 16:16:50 GMT
https://trac.sagemath.org/ticket/11941#comment:7
https://trac.sagemath.org/ticket/11941#comment:7
<p>
Ugh, all my comments have screwed up formatting. I wish I could edit them (can I?). I should have been:
</p>
<pre class="wiki">
(%i8) domain:complex;
(%o8) complex
(%i9) solve(x^3+1=0,x);
(%o9) [x = ((-1)^(1/3)*sqrt(3)*%i-(-1)^(1/3))/2,
x = -((-1)^(1/3)*sqrt(3)*%i+(-1)^(1/3))/2,x = (-1)^(1/3)]
</pre>
TicketkcrismanThu, 17 Nov 2011 16:19:53 GMT
https://trac.sagemath.org/ticket/11941#comment:8
https://trac.sagemath.org/ticket/11941#comment:8
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/11941#comment:6" title="Comment 6">mboratko</a>:
</p>
<blockquote class="citation">
<p>
It seems that sage sets domain: complex (I was made aware of this by burcin in IRC). You do get this result as follows:
</p>
</blockquote>
<p>
Yes, we do, but I didn't bother checking that. Good work.
</p>
<p>
So of course now the question becomes what the "right" thing to do is? I don't think we want to set and unset <code>domain:real/complex</code> in Maxima every time we use <code>solve</code>, because presumably this would break other things. Or? At any rate we definitely need to keep <code>domain:complex</code> in general, if I recall correctly other problems that occur without it.
</p>
<pre class="wiki">(%i1) (-1)^(1/3);
(%o1) - 1
(%i2) domain:complex;
(%o2) complex
(%i3) (-1)^(1/3);
1/3
(%o3) (- 1)
</pre><p>
Typically we would want the latter answer, e.g in
</p>
<pre class="wiki">sage: a = (-1)^(1/3)
sage: a.simplify()
(-1)^(1/3)
</pre>
TicketmboratkoThu, 17 Nov 2011 16:53:47 GMT
https://trac.sagemath.org/ticket/11941#comment:9
https://trac.sagemath.org/ticket/11941#comment:9
<p>
I've attached a fairly limited workaround, but it does give (semi?) desirable behavior:
</p>
<pre class="wiki">sage: assume(x, 'real')
sage: solve(x^3+1==0,x)
[x == (-1)^(1/3)]
</pre><p>
If it is then up to the user to interpret the result, this seems ok. If they want to programmatically use the result later, it's really no good. I suppose if they were aware of this and wanted to fix the result to be real then they could manually employ the same method as in my patch to the result.
</p>
<p>
I guess, as a more general question, do we want the results of solve to always return a result with domain: real, and if so can we make this change for just this function? If there is a use case where this is undesirable, then I think the only option is make the assumptions file not only check the returned value, but also modify it (so that if assume is real, then the value is actually replaced with the real values). Perhaps another option is to allow the user to specify the desired domain of results from solve.
</p>
TicketmboratkoThu, 17 Nov 2011 17:16:55 GMTattachment set
https://trac.sagemath.org/ticket/11941
https://trac.sagemath.org/ticket/11941
<ul>
<li><strong>attachment</strong>
set to <em>trac_11941.patch</em>
</li>
</ul>
<p>
limited workaround for assumption and solve
</p>
TicketjdemeyerTue, 13 Aug 2013 15:35:53 GMTmilestone changed
https://trac.sagemath.org/ticket/11941#comment:10
https://trac.sagemath.org/ticket/11941#comment:10
<ul>
<li><strong>milestone</strong>
changed from <em>sage-5.11</em> to <em>sage-5.12</em>
</li>
</ul>
Ticketvbraun_spamThu, 30 Jan 2014 21:20:52 GMTmilestone changed
https://trac.sagemath.org/ticket/11941#comment:11
https://trac.sagemath.org/ticket/11941#comment:11
<ul>
<li><strong>milestone</strong>
changed from <em>sage-6.1</em> to <em>sage-6.2</em>
</li>
</ul>
Ticketvbraun_spamTue, 06 May 2014 15:20:58 GMTmilestone changed
https://trac.sagemath.org/ticket/11941#comment:12
https://trac.sagemath.org/ticket/11941#comment:12
<ul>
<li><strong>milestone</strong>
changed from <em>sage-6.2</em> to <em>sage-6.3</em>
</li>
</ul>
Ticketvbraun_spamSun, 10 Aug 2014 16:51:03 GMTmilestone changed
https://trac.sagemath.org/ticket/11941#comment:13
https://trac.sagemath.org/ticket/11941#comment:13
<ul>
<li><strong>milestone</strong>
changed from <em>sage-6.3</em> to <em>sage-6.4</em>
</li>
</ul>
TicketrwsSat, 31 Jan 2015 09:05:26 GMTstatus changed
https://trac.sagemath.org/ticket/11941#comment:14
https://trac.sagemath.org/ticket/11941#comment:14
<ul>
<li><strong>status</strong>
changed from <em>new</em> to <em>needs_review</em>
</li>
</ul>
<p>
You need to set "needs review" if you want someone to look at your code. I'll do that now.
</p>
TicketrwsWed, 11 Mar 2015 17:22:31 GMTstatus, upstream changed; work_issues set
https://trac.sagemath.org/ticket/11941#comment:15
https://trac.sagemath.org/ticket/11941#comment:15
<ul>
<li><strong>status</strong>
changed from <em>needs_review</em> to <em>needs_work</em>
</li>
<li><strong>work_issues</strong>
set to <em>report upstream</em>
</li>
<li><strong>upstream</strong>
changed from <em>N/A</em> to <em>Not yet reported upstream; Will do shortly.</em>
</li>
</ul>
<p>
As Karl-Dieter says in <a class="ticket" href="https://trac.sagemath.org/ticket/11941#comment:8" title="Comment 8">comment:8</a> it is not desirable to switch back and forth with Maxima domain commands, so I don't think your patch is the right way to solve it, esp. since you don't get the right result, either. Presumably the failure of assumption should be reported upstream.
</p>
TicketrwsMon, 05 Dec 2016 06:51:53 GMTmilestone changed; cc set
https://trac.sagemath.org/ticket/11941#comment:16
https://trac.sagemath.org/ticket/11941#comment:16
<ul>
<li><strong>cc</strong>
<em>pelegm</em> added
</li>
<li><strong>milestone</strong>
changed from <em>sage-6.4</em> to <em>sage-7.6</em>
</li>
</ul>
<p>
As found in <a class="closed ticket" href="https://trac.sagemath.org/ticket/22017" title="defect: Unreadable real solution for a very simple equation (closed: invalid)">#22017</a> SymPy gets it right. If SymPy is better than Maxima with symbolic polynomial roots then I think we should switch to SymPy for the special case.
</p>
TicketpelegmMon, 05 Dec 2016 07:19:21 GMT
https://trac.sagemath.org/ticket/11941#comment:17
https://trac.sagemath.org/ticket/11941#comment:17
<p>
Note that sympy also incorporates assumptions into the solver:
</p>
<div class="wiki-code"><div class="code"><pre>In <span class="p">[</span><span class="mi">2</span><span class="p">]:</span> x <span class="o">=</span> var<span class="p">(</span><span class="s">'x'</span><span class="p">,</span> real<span class="o">=</span><span class="bp">True</span><span class="p">)</span>
In <span class="p">[</span><span class="mi">3</span><span class="p">]:</span> solve<span class="p">(</span>Eq<span class="p">(</span>x<span class="o">**</span><span class="mi">3</span><span class="p">,</span> <span class="o">-</span><span class="mi">8</span><span class="p">))</span>
Out<span class="p">[</span><span class="mi">3</span><span class="p">]:</span> <span class="p">[</span><span class="o">-</span><span class="mi">2</span><span class="p">]</span>
</pre></div></div>
TicketrwsMon, 05 Dec 2016 07:39:51 GMT
https://trac.sagemath.org/ticket/11941#comment:18
https://trac.sagemath.org/ticket/11941#comment:18
<p>
Switching to SymPy would also depend on something like <a class="closed ticket" href="https://trac.sagemath.org/ticket/22024" title="defect: symbolic placeholder for complex root (closed: fixed)">#22024</a>.
</p>
TicketrwsThu, 02 Nov 2017 07:09:23 GMT
https://trac.sagemath.org/ticket/11941#comment:19
https://trac.sagemath.org/ticket/11941#comment:19
<p>
With <a class="closed ticket" href="https://trac.sagemath.org/ticket/22024" title="defect: symbolic placeholder for complex root (closed: fixed)">#22024</a> we have now:
</p>
<pre class="wiki">sage: solve(x^3+1==0,x,algorithm='sympy')
[x == -1, x == -1/2*I*sqrt(3) + 1/2, x == 1/2*I*sqrt(3) + 1/2]
sage: solve(x^3+1==0,x,algorithm='sympy', domain='real')
[x == -1]
</pre><p>
which, as I understand the notorious discussion, isn't right either because in the complex domain Maxima's results are not reproduced?
</p>
TicketcharpentTue, 14 Nov 2017 15:40:04 GMT
https://trac.sagemath.org/ticket/11941#comment:20
https://trac.sagemath.org/ticket/11941#comment:20
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/11941#comment:19" title="Comment 19">rws</a>:
</p>
<blockquote class="citation">
<p>
With <a class="closed ticket" href="https://trac.sagemath.org/ticket/22024" title="defect: symbolic placeholder for complex root (closed: fixed)">#22024</a> we have now:
</p>
<pre class="wiki">sage: solve(x^3+1==0,x,algorithm='sympy')
[x == -1, x == -1/2*I*sqrt(3) + 1/2, x == 1/2*I*sqrt(3) + 1/2]
sage: solve(x^3+1==0,x,algorithm='sympy', domain='real')
[x == -1]
</pre><p>
which, as I understand the notorious discussion, isn't right either because in the complex domain Maxima's results are not reproduced?
</p>
</blockquote>
<p>
Hmmm... things seems to have changed on Maxima's front. Compare :
</p>
<pre class="wiki">charpent@p-202-021:~$ sage -maxima
;;; Loading #P"/usr/local/sage-8/local/lib/ecl/sb-bsd-sockets.fas"
;;; Loading #P"/usr/local/sage-8/local/lib/ecl/sockets.fas"
;;; Loading #P"/usr/local/sage-8/local/lib/ecl/defsystem.fas"
;;; Loading #P"/usr/local/sage-8/local/lib/ecl/cmp.fas"
Maxima 5.39.0 http://maxima.sourceforge.net
using Lisp ECL 16.1.2
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) forget();
(%o1) []
(%i2) solve(x^3+1=0,x);
sqrt(3) %i - 1 sqrt(3) %i + 1
(%o2) [x = - --------------, x = --------------, x = - 1]
2 2
(%i3) declare(x,real);
(%o3) done
(%i4) solve(x^3+1=0,x);
sqrt(3) %i - 1 sqrt(3) %i + 1
(%o4) [x = - --------------, x = --------------, x = - 1]
2 2
(%i5) quit();
</pre><p>
[ Note : this is "our" Maxima ; but Maxima 5.40.0 as packaged in Debian and Cocalc's version both give the same answers... ]
</p>
<p>
and :
</p>
<pre class="wiki">charpent@p-202-021:~$ sage
┌────────────────────────────────────────────────────────────────────┐
│ SageMath version 8.1.rc0, Release Date: 2017-11-08 │
│ Type "notebook()" for the browser-based notebook interface. │
│ Type "help()" for help. │
└────────────────────────────────────────────────────────────────────┘
┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓
┃ Warning: this is a prerelease version, and it may be unstable. ┃
┗━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┛
sage: forget();
sage: solve(x^3+1==0,x)
[x == 1/2*I*sqrt(3)*(-1)^(1/3) - 1/2*(-1)^(1/3), x == -1/2*I*sqrt(3)*(-1)^(1/3) - 1/2*(-1)^(1/3), x == (-1)^(1/3)]
sage: assume(x,"real")
sage: solve(x^3+1==0,x)
[]
sage: quit
Exiting Sage (CPU time 0m1.63s, Wall time 1m2.97s).
</pre><p>
Maxima's second answer may be disputable (it doesn't account for the declaration of <code>x</code> as real), but Sage's is indisputably wrong, <em>wrong</em>, <strong>wrong</strong>.
</p>
<p>
I'm painfully tempted to file a new ticket and flag it as critical. Advice ?
</p>
TicketkcrismanTue, 14 Nov 2017 15:45:01 GMT
https://trac.sagemath.org/ticket/11941#comment:21
https://trac.sagemath.org/ticket/11941#comment:21
<p>
As pointed out <a class="ticket" href="https://trac.sagemath.org/ticket/11941#comment:6" title="Comment 6">earlier</a>, this is due to <code>domain:complex</code>:
</p>
<pre class="wiki">Maxima 5.39.0 http://maxima.sourceforge.net
using Lisp ECL 16.1.2
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) domain:complex;
(%o1) complex
(%i2) solve(x^3+1=0,x);
1/3 1/3
(- 1) sqrt(3) %i - (- 1)
(%o2) [x = ------------------------------,
2
1/3 1/3
(- 1) sqrt(3) %i + (- 1) 1/3
x = - ------------------------------, x = (- 1) ]
2
</pre>
TicketkcrismanTue, 14 Nov 2017 15:46:25 GMT
https://trac.sagemath.org/ticket/11941#comment:22
https://trac.sagemath.org/ticket/11941#comment:22
<p>
And if you can figure out how to deal with this - in Maxima or elsewhere - please do! I don't think this was ever reported upstream.
</p>
TicketcharpentTue, 14 Nov 2017 16:24:58 GMT
https://trac.sagemath.org/ticket/11941#comment:23
https://trac.sagemath.org/ticket/11941#comment:23
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/11941#comment:21" title="Comment 21">kcrisman</a>:
</p>
<blockquote class="citation">
<p>
As pointed out <a class="ticket" href="https://trac.sagemath.org/ticket/11941#comment:6" title="Comment 6">earlier</a>, this is due to <code>domain:complex</code>:
</p>
<pre class="wiki">Maxima 5.39.0 http://maxima.sourceforge.net
using Lisp ECL 16.1.2
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) domain:complex;
(%o1) complex
(%i2) solve(x^3+1=0,x);
1/3 1/3
(- 1) sqrt(3) %i - (- 1)
(%o2) [x = ------------------------------,
2
1/3 1/3
(- 1) sqrt(3) %i + (- 1) 1/3
x = - ------------------------------, x = (- 1) ]
2
</pre></blockquote>
<p>
OK. I agree that this is only <em>disputable</em>.
</p>
<p>
My beef is with Sage's <em>second</em> answer, which tells us that Sage is unable to find a real root to <code>x^3+1==O</code>. Maxima returns a list of three <em>candidate</em> answers, among whom two turn out to be unacceptable. Sage turns out <em>no</em> candidate.
</p>
TicketkcrismanTue, 14 Nov 2017 18:16:50 GMT
https://trac.sagemath.org/ticket/11941#comment:24
https://trac.sagemath.org/ticket/11941#comment:24
<p>
You can definitely feel free to fix it or report upstream though! It is ugly to say the least.
</p>
<p>
The second answer I can't quite explain. Typically <code>solve</code> has nothing to do with our assumptions - least of all the <code>declare</code> syntax Maxima uses for what we do with things like assuming real or integer. But we do have some minimal checking (see below).
</p>
<pre class="wiki">(%i3) declare(x,real);
(%o3) done
(%i4) solve(x^3+1=0,x);
1/3 1/3
(- 1) sqrt(3) %i - (- 1)
(%o4) [x = ------------------------------,
2
1/3 1/3
(- 1) sqrt(3) %i + (- 1) 1/3
x = - ------------------------------, x = (- 1) ]
2
</pre><p>
However, note that without <code>domain:complex</code> we get
</p>
<pre class="wiki">(%i1) declare(x,real);
(%o1) done
(%i2) solve(x^3+1=0,x);
sqrt(3) %i - 1 sqrt(3) %i + 1
(%o2) [x = - --------------, x = --------------, x = - 1]
2 2
</pre><p>
as perhaps noted above.
</p>
<p>
<a class="ext-link" href="https://github.com/sagemath/sage/blob/master/src/sage/symbolic/expression.pyx#L11436"><span class="icon"></span>Here is the relevant code</a> for how Sage checks for assumptions with solving.
</p>
<pre class="wiki"> # make sure all the assumptions are satisfied
from sage.symbolic.assumptions import assumptions
to_check = assumptions()
if to_check:
for ix, soln in reversed(list(enumerate(X))):
if soln.lhs().is_symbol():
if any([a.contradicts(soln) for a in to_check]):
del X[ix]
if multiplicities:
del ret_multiplicities[ix]
continue
</pre><p>
Apparently something is going wrong here with <code>x == -1</code>, but I'm not sure what.
</p>
TicketcharpentTue, 14 Nov 2017 19:36:15 GMT
https://trac.sagemath.org/ticket/11941#comment:25
https://trac.sagemath.org/ticket/11941#comment:25
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/11941#comment:24" title="Comment 24">kcrisman</a>:
</p>
<blockquote class="citation">
<p>
You can definitely feel free to fix it or report upstream though! It is ugly to say the least.
</p>
</blockquote>
<p>
Inded. But the point of <em>not</em> replacing (-1)^(1/n) by -1 is well taken : the first expression may, after all, be <em>any</em> <em>n</em>th root of -1 (see an example below).
</p>
<p>
So I'm not sure it's a bug. But yes, it's ugly as hell... (more tolerable in \LaTeX...).
</p>
<blockquote class="citation">
<p>
The second answer I can't quite explain. Typically <code>solve</code> has nothing to do with our assumptions - least of all the <code>declare</code> syntax Maxima uses for what we do with things like assuming real or integer. But we do have some minimal checking (see below).
</p>
<pre class="wiki">(%i3) declare(x,real);
(%o3) done
(%i4) solve(x^3+1=0,x);
1/3 1/3
(- 1) sqrt(3) %i - (- 1)
(%o4) [x = ------------------------------,
2
1/3 1/3
(- 1) sqrt(3) %i + (- 1) 1/3
x = - ------------------------------, x = (- 1) ]
2
</pre><p>
However, note that without <code>domain:complex</code> we get
</p>
<pre class="wiki">(%i1) declare(x,real);
(%o1) done
(%i2) solve(x^3+1=0,x);
sqrt(3) %i - 1 sqrt(3) %i + 1
(%o2) [x = - --------------, x = --------------, x = - 1]
2 2
</pre><p>
as perhaps noted above.
</p>
<p>
<a class="ext-link" href="https://github.com/sagemath/sage/blob/master/src/sage/symbolic/expression.pyx#L11436"><span class="icon"></span>Here is the relevant code</a> for how Sage checks for assumptions with solving.
</p>
<pre class="wiki"> # make sure all the assumptions are satisfied
from sage.symbolic.assumptions import assumptions
to_check = assumptions()
if to_check:
for ix, soln in reversed(list(enumerate(X))):
if soln.lhs().is_symbol():
if any([a.contradicts(soln) for a in to_check]):
del X[ix]
if multiplicities:
del ret_multiplicities[ix]
continue
</pre><p>
Apparently something is going wrong here with <code>x == -1</code>, but I'm not sure what.
</p>
</blockquote>
<p>
I think that our code for testing that an expression is real is too weak. After all,
</p>
<pre class="wiki">sage: ((-1)^(1/3)).is_real()
False
</pre><p>
A workaround is to force the evaluation of each root "in Sage terms", as demonstrates the following crock :
</p>
<pre class="wiki">sage: Sols=solve(x^3+1==0,x);Sols
[x == 1/2*I*sqrt(3)*(-1)^(1/3) - 1/2*(-1)^(1/3), x == -1/2*I*sqrt(3)*(-1)^(1/3) - 1/2*(-1)^(1/3), x == (-1)^(1/3)]
sage: [t.rhs().is_real() for t in Sols]
[False, False, False]
</pre><p>
None of these roots is <em>known as</em> real (in direct contradiction of d'Alembert's theorem, no less...). Try to force a re-evaluation of these expressions :
</p>
<pre class="wiki">sage: [t.rhs().real_part()+I*t.rhs().imag_part() for t in Sols]
[-1, -1/2*I*sqrt(3) + 1/2, 1/2*I*sqrt(3) + 1/2]
</pre><blockquote>
<p>
(Note that this implies that, in that specific case, <code>(-1)^(1/3)</code> is <code>(I*sqrt(3)+1)/2</code>...)
</p>
</blockquote>
<p>
Now, these re-evaluated quantities <em>can</em> be effectively tested for "reality" :
</p>
<pre class="wiki">sage: [(t.rhs().real_part()+I*t.rhs().imag_part()).is_real() for t in Sols]
[True, False, False]
</pre><p>
I do not know what code uses the <code>.contradict()</code> method for the assertion <code>x is real</code>, but it may fall in the same trap.
</p>
<p>
The problem is now to know <em>what code</em> is to be fixed : assumptions ? Or more general algebraic code ? Is this problem specific to Maxima-generated expressions, or more general ? How to force re-evaluation (<code>real_part()</code> and <code>imag_part()</code> may be highly nontrivial, or even impossible for some expressions) ?
</p>
<p>
Advice more than welcome...
</p>
TicketcharpentWed, 15 Nov 2017 20:38:15 GMT
https://trac.sagemath.org/ticket/11941#comment:26
https://trac.sagemath.org/ticket/11941#comment:26
<p>
One more data point : Maxima seems to be able to solve the specific test which is problematic for the current Sage assumption code. Consider :
</p>
<pre class="wiki">sage: ## A few abbreviations, I'm lazy
sage: def mrhs(x):return(maxima_lib.rhs(x))
sage: def mreal(x):return(maxima_lib.featurep(x,"real"))
sage: def msolve(e,v):return(maxima_lib.solve(*[t._maxima_lib_() for t in [e,v]]
....: ))
sage: assumptions()
[]
sage: maxima_lib.facts()
[kind(sinh,one_to_one),kind(log,one_to_one),kind(tanh,one_to_one),kind(log,increasing)]
sage: [mreal(mrhs(t)) for t in msolve(x^3+1==0,x)]
[true, false, false]
</pre><p>
Questions :
</p>
<ul><li>Should we use this ? If so
<ul><li>where ?
</li><li>with what generality ?
</li></ul></li></ul><p>
Advice necessary...
</p>
TicketrwsThu, 16 Nov 2017 06:52:17 GMT
https://trac.sagemath.org/ticket/11941#comment:27
https://trac.sagemath.org/ticket/11941#comment:27
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/11941#comment:25" title="Comment 25">charpent</a>:
</p>
<blockquote class="citation">
<p>
I think that our code for testing that an expression is real is too weak. After all,
</p>
<pre class="wiki">sage: ((-1)^(1/3)).is_real()
False
</pre></blockquote>
<p>
The False from is_real (and any of these functions) just means "I don't know" in absence of a Python tri-state logic. It may be possible to return Unknown here by implementing <code>is_complex</code>.
</p>
TicketrwsThu, 16 Nov 2017 06:56:19 GMT
https://trac.sagemath.org/ticket/11941#comment:28
https://trac.sagemath.org/ticket/11941#comment:28
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/11941#comment:26" title="Comment 26">charpent</a>:
</p>
<blockquote class="citation">
<pre class="wiki">sage: [mreal(mrhs(t)) for t in msolve(x^3+1==0,x)]
[true, false, false]
</pre></blockquote>
<p>
Is that root really real? In what domain is Maxima at that point?
</p>
TicketcharpentThu, 16 Nov 2017 07:52:25 GMT
https://trac.sagemath.org/ticket/11941#comment:29
https://trac.sagemath.org/ticket/11941#comment:29
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/11941#comment:28" title="Comment 28">rws</a>:
</p>
<blockquote class="citation">
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/11941#comment:26" title="Comment 26">charpent</a>:
</p>
<blockquote class="citation">
<pre class="wiki">sage: [mreal(mrhs(t)) for t in msolve(x^3+1==0,x)]
[true, false, false]
</pre></blockquote>
<p>
Is that root really real?
</p>
</blockquote>
<p>
That's what Maxima says. Whatever it does is probably more credible than what we do (see below).
</p>
<blockquote class="citation">
<p>
In what domain is Maxima at that point?
</p>
</blockquote>
<p>
Look for yourself :
</p>
<pre class="wiki">sage: from sage.interfaces.maxima_lib import maxima_lib as ml
sage: def mrhs(x):return(ml.rhs(x))
sage: def msolve(e,v):return(ml.solve(*[t._maxima_lib_() for t in [e,v]]))
sage: def mreal(x):return(ml.featurep(x,"real"))
sage: ml.ev("domain")
complex
sage: sol=msolve(x^3+1==0,x);sol
[_SAGE_VAR_x=((-1)^(1/3)*sqrt(3)*%i-(-1)^(1/3))/2,_SAGE_VAR_x=-((-1)^(1/3)*sqrt(3)*%i+(-1)^(1/3))/2,_SAGE_VAR_x=(-1)^(1/3)]
sage: ml.ev("domain")
complex
sage: [mreal(mrhs(t)) for t in sol]
[true, false, false]
sage: [mrhs(t).sage().n() for t in sol]
[-1.00000000000000 + 1.11022302462516e-16*I,
0.500000000000000 - 0.866025403784439*I,
0.500000000000000 + 0.866025403784439*I]
</pre><p>
The numerical values <em>tend to confirm</em> that the first root is real.
</p>
<p>
And <strong>that</strong> is our problem : the <code>.contradicts</code> code (in <code>$SAGE_ROOT/src/sage/symbolic/assumptions.py</code>) is piss-poor : it coerces the value to be tested to <code>CC</code> and tests if the resulting coercion belongs to <code>RR</code>. Aaaaarghhh...
</p>
<p>
That said, I stumbled on another problem. Our current code doesn't pass <em>declarations</em> to Maxima correctly :
</p>
<pre class="wiki">age: assume(z,"integer")
sage: assume(z>0)
sage: assumptions()
[z is integer, z > 0]
sage: ml.facts()
[kind(sinh,one_to_one),kind(log,one_to_one),kind(tanh,one_to_one),kind(log,increasing),_SAGE_VAR_z>0]
sage: maxima_calculus("facts()")
[kind(sinh,one_to_one),kind(log,one_to_one),kind(tanh,one_to_one),kind(log,increasing),_SAGE_VAR_z>0]
</pre><p>
I haven't (yet) checked trac to see if this is known. If not, that's a nice ticket to file.
</p>
<p>
I don't (yet) have a solution.
</p>
Ticket