solve should not convert floating point to rationals when solving

[videlec]

The following (obtained on 8.2.beta6) is very bad

sage: solve(1.0 * x^2 - 1.5*x + 2.0, x)
[x == -1/4*I*sqrt(23) + 3/4, x == 1/4*I*sqrt(23) + 3/4]

The input is an equation with floating point numbers (likely to be approximations). The answer obtained are exact symbolic expressions... it becomes crazy with

sage: solve(1.13157771r * x^2 - 1.2241351312r*x + 2.0000401231r, x)
[x == -1/1516050631807733202922*I*sqrt(3389945193423588592231412963046087595757907) + 190736613975697/352629861450338,
 x == 1/1516050631807733202922*I*sqrt(3389945193423588592231412963046087595757907) + 190736613975697/352629861450338]

CC: @sagetrac-tmonteil @rwst

Component: basic arithmetic

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

[videlec]

Description changed:

--- 
+++ 
@@ -4,4 +4,11 @@
 sage: solve(1.0 * x^2 - 1.5*x + 2.0, x)
 [x == -1/4*I*sqrt(23) + 3/4, x == 1/4*I*sqrt(23) + 3/4]
 ```
-The input is an equation with floating point numbers (likely to be approximations). The answer obtained are exact symbolic expressions.
+The input is an equation with floating point numbers (likely to be approximations). The answer obtained are exact symbolic expressions... it becomes crazy with
+
+```
+sage: solve(1.13157771r * x^2 - 1.2241351312r*x + 2.0000401231r, x)
+[x == -1/1516050631807733202922*I*sqrt(3389945193423588592231412963046087595757907) + 190736613975697/352629861450338,
+ x == 1/1516050631807733202922*I*sqrt(3389945193423588592231412963046087595757907) + 190736613975697/352629861450338]
+```
+

[rwst]

comment:2

We are following Maxima's lead here:

(%i3) solve(1.0 * x^2 - 1.5*x + 2.0, x);

rat: replaced 2.0 by 2/1 = 2.0

rat: replaced -1.5 by -3/2 = -1.5

rat: replaced 1.0 by 1/1 = 1.0
                        sqrt(23) %i - 3      sqrt(23) %i + 3
(%o3)            [x = - ---------------, x = ---------------]
                               4                    4

Note that our documentation states: Whenever possible, answers will be symbolic, but with systems of equations, at times approximations will be given by Maxima, due to the underlying algorithm.

The workaround would be to use SymPy:

sage: solve(1.0 * x^2 - 1.5*x + 2.0, x, algorithm='sympy')
[x == (0.750000000000000 - 1.19895788082818*I),
 x == (0.750000000000000 + 1.19895788082818*I)]
sage: solve(1.13157771r * x^2 - 1.2241351312r*x + 2.0000401231r, x, algorithm='s
....: ympy')
[x == (0.540897509902347 - 1.21445837087478*I),
 x == (0.540897509902347 + 1.21445837087478*I)]

[videlec]

comment:3

Replying to @rwst:

We are following Maxima's lead here:

<SNIP>

Note that our documentation states: Whenever possible, answers will be symbolic, but with systems of equations, at times approximations will be given by Maxima, due to the underlying algorithm.

Though it make no sense in the present situation! And the meaning of "symbolic" is anyway very vague.

The workaround would be to use SymPy:

<SNIP>

Much better :-)

[rwst]

comment:4

Finally, see polynomial_element.pyx the roots() member where I inserted fast code that avoids calling Maxima for degree 2. I would be interested in taking over solving every polynomial that's possible in Pynac if I had support for it. It would make QQbar faster, because QQbar uses that code in roots() too.

[videlec]

comment:5

Replying to @rwst:

Finally, see polynomial_element.pyx the roots() member where I inserted fast code that avoids calling Maxima for degree 2. I would be interested in taking over solving every polynomial that's possible in Pynac if I had support for it. It would make QQbar faster, because QQbar uses that code in roots() too.

This code in polynomial_element.pyx is numerically unstable. You should not try to reinvent the wheel here

sage: coeffs = [85.37r, -59.22r, 10.27r]
sage: RDFx = RDF['x']
sage: pRDF = RDFx(coeffs)
sage: r1, r2 = pRDF.roots(multiplicities=False)
sage: coeffs[2] * r1^2 + coeffs[1] * r1 + coeffs[0]
0.0

the same with SR

sage: coeffs = [85.37r, -59.22r, 10.27r]
sage: SRx = SR['x']
sage: pSR = SRx(coeffs)
sage: r1, r2 = pSR.roots(multiplicities=False)
sage: coeffs[2] * r1^2 + coeffs[1] * r1 + coeffs[0]
1.4210854715202004e-14

Moreover having import SR any time roots is called is slowing down everything and is a very bad practice. You should have used the special method _roots_univariate_polynomial that has to be implemented directly into the ring. See #24942.

[videlec]

comment:6

Replying to @videlec:

Replying to @rwst:
Moreover having import SR any time roots is called is slowing down everything and is a very bad practice. You should have used the special method _roots_univariate_polynomial that has to be implemented directly into the ring. See #24942.

To be clear: anything specific to SR has to go in symbolic/ and not invade the whole Sage code. This code was badly intrusive.

[rwst]

comment:7

Replying to @videlec:

To be clear: anything specific to SR has to go in symbolic/ and not invade the whole Sage code. This code was badly intrusive.

To be clear as well: that specialization was not introduced by me. I merely added a shortcut that no longer used Maxima.

[videlec]

comment:8

Replying to @rwst:

Replying to @videlec:

To be clear: anything specific to SR has to go in symbolic/ and not invade the whole Sage code. This code was badly intrusive.

To be clear as well: that specialization was not introduced by me. I merely added a shortcut that no longer used Maxima.

Actually, this is not only true for SR but for some other rings as well! This roots code is a mess...

[rwst]

comment:9

Coming back to the original problem the solution would be to insert a shortcut for inexact polynomial input in solve that uses other methods. The default algorithm of solve at the moment is Maxima so the documentation should explain that with such input other methods are used. Specifically, I would not try to implement this by turning a switch on/off in Maxima because Sage's own root finding methods are probably better suited for all the possible inexact Sage number types.