Opened 15 months ago

Last modified 15 months ago

#24939 new defect

solve should not convert floating point to rationals when solving

Reported by: vdelecroix Owned by:
Priority: minor Milestone: sage-8.2
Component: basic arithmetic Keywords:
Cc: tmonteil, rws Merged in:
Authors: Reviewers:
Report Upstream: N/A Work issues:
Branch: Commit:
Dependencies: Stopgaps:

Description (last modified by vdelecroix)

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]

Change History (10)

comment:1 Changed 15 months ago by vdelecroix

  • Description modified (diff)

comment:2 follow-up: Changed 15 months 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 15 months ago by vdelecroix

Replying to rws:

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: Changed 15 months 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: Changed 15 months ago by vdelecroix

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

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.

comment:6 in reply to: ↑ 5 ; follow-up: Changed 15 months ago by vdelecroix

Replying to 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: Changed 15 months ago by rws

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

comment:8 in reply to: ↑ 7 Changed 15 months ago by vdelecroix

Replying to rws:

Replying to 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 15 months 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 15 months ago by rws

See Expression.is_exact().

Note: See TracTickets for help on using tickets.