Opened 3 years ago

# solve should not convert floating point to rationals when solving

Reported by: Owned by: vdelecroix minor sage-8.2 basic arithmetic tmonteil, rws N/A

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]
```

### comment:1 Changed 3 years ago by vdelecroix

• Description modified (diff)

### comment:2 follow-up: ↓ 3 Changed 3 years ago by rws

• Priority changed from major to minor

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)]
```

### comment:3 in reply to: ↑ 2 Changed 3 years ago by vdelecroix

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 :-)

### comment:4 follow-up: ↓ 5 Changed 3 years ago by rws

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.

### comment:5 in reply to: ↑ 4 ; follow-up: ↓ 6 Changed 3 years ago by vdelecroix

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 * r1^2 + coeffs * r1 + coeffs
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 * r1^2 + coeffs * r1 + coeffs
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.

### comment:6 in reply to: ↑ 5 ; follow-up: ↓ 7 Changed 3 years ago by vdelecroix

Replying to rws: 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.

### comment:7 in reply to: ↑ 6 ; follow-up: ↓ 8 Changed 3 years ago by rws

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.

### comment:8 in reply to: ↑ 7 Changed 3 years ago by vdelecroix

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...

### comment:9 Changed 3 years ago by rws

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.

### comment:10 Changed 3 years ago by rws

See `Expression.is_exact()`.

Note: See TracTickets for help on using tickets.