Cancel with floating point numbers

[Bartosz]

Hi all,

I have been trying to simplify expressions with floating point numbers
and I noticed that the results are quite unpredictable. In many cases
the calculation took indefinite time (until I killed it) or produced
wrong results (infinity or other incorrect expressions). Here is an
example:

simplify of:

2
2.1 - 2.1⋅y + 2⋅x
───────────
3
x

returns:

2
─
x

While:

2
2.1⋅x⋅(1 - y) + 2⋅y⋅x
──────────────────────
3
x

returns infinity.

My workaround is now to use rational numbers instead:

21⋅x⋅(1 - y) 2
──────────── + 2⋅y⋅x
10
─────────────────────
3
x

but it would be better if this could be done automatically. Is there a
simple fix to this problem?

Thanks for your help and sympy!

Cheers,

Bartosz

[Andy Ray Terrel]

Quite strange indeed.

On my Mac Python 2.7, sympy.__version__ = 0.6.7 I get:

In [17]: expr = (2.1 - 2.1*y**2 + 2*x)/x**3

In [18]: sympy.simplify(expr)
Out[18]: -inf

In [19]: expr = (2.1*x*(1-y)**2+2*y*x)/ x**3

In [20]: sympy.simplify(expr)
Out[20]: (2.1 - 2.2*y + 2.1*y**2)/x**2

In [21]: expr = (sympy.Rational(21,10)*x*(1-y)**2 + 2*y*x) / X**3

In [22]: sympy.simplify(expr)
Out[22]: (105*x - 110*x*y + 105*x*y**2)/(50*X**3)

In the dev version I just get completely wrong results.

-- Andy

> --
> You received this message because you are subscribed to the Google Groups "sympy" group.
> To post to this group, send email to sy...@googlegroups.com.
> To unsubscribe from this group, send email to sympy+un...@googlegroups.com.
> For more options, visit this group at http://groups.google.com/group/sympy?hl=en.
>
>

[Chris Smith]

Andy Ray Terrel wrote:
> Quite strange indeed.
>
> On my Mac Python 2.7, sympy.__version__ = 0.6.7 I get:
>
> In [17]: expr = (2.1 - 2.1*y**2 + 2*x)/x**3
>
> In [18]: sympy.simplify(expr)
> Out[18]: -inf
>

This raises an error for me:

raise ValueError("can't raise polynomial to a negative power")
ValueError: can't raise polynomial to a negative power

But if I nsimplify(..., rational=True) everthing works ok:

h[12] >>> e1= (2.1 - 2.1*y**2 + 2*x)/x**3
h[12] >>> e2= (2.1*x*(1-y)**2+2*y*x)/ x**3
h[12] >>> e3= (sympy.Rational(21,10)*x*(1-y)**2 + 2*y*x) / X**3
h[12] >>> for e in [e1,e2,e3]:
... print simplify(nsimplify(e,rational=True))
...
(20*x - 21*y**2 + 21)/(10*x**3)
(21*y**2 - 22*y + 21)/(10*x**2)
x*(21*y**2 - 22*y + 21)/(10*X**3)

CAUTION: you used an X instead of an x in the last equation.

[Tom Bachmann]

Presumably mateusz (or aaron) know better than me, but the problem is
simply that polys does not really work well with Floats.

There is also an open issue (2384) related to this.

[Aaron S. Meurer]

That's basically correct. The problem is that certain polynomial algorithms (in this case, polynomial division) rely on the ability to determine zero equivalence in the coefficients. But floating point numbers, as we all know, have the problem where they will often give something like 1e-14 instead of 0.

The solution is to replace the RR domain with RR(eps) domain, where eps is the cutoff for zero (if 0 <= a < eps, then consider a == 0).

By the way, if anyone is interested, the problem comes from while loop at line 1268 of densearith.py. The loop breaks when the degree of r is less than the degree of g, but in this case, on the final iteration where it should be less, the leading term of r has a coefficient somewhere on the order of 1e-13 instead of 0, so it goes through the loop one too many times. Since the loop decrements N each time, this ends up becoming negative. This is the "negative power" it complains about later.

Mateusz, would a suitable temporary fix for this be to put chop=True in key places in the RR domain? I think a lot of people use floating point coefficients, and they won't be happy to see errors like this one.

Aaron Meurer

[Bartosz Telenczuk]

Hi all,

Thanks for the quick reply. I thought that a sensible solution might
be to introduce a new floating point type (for example, Float) which
does the round off by itself (with given precision). Basically, one
could write the equation like this:

expr = (Float(2.1,digits=1)*(1-y**2)+2*x)/x**3

Is it now possible in sympy? If not, would it be hard to implement it?

BTW Your are doing great job with sympy! Thanks a lot!

Cheers,

Bartosz

[Mateusz Paprocki]

Hi,

One problem was with PRS GCD algorithm and was fixed in:

 https://github.com/sympy/sympy/commit/2f6031a364cc33612285e704e2615e3cd5682aa3

 

The other problem (mainly "can't raise polynomial to a negative power") is related to poor floating point number support in sympy.polys. The exception is being raised in dmp_prem() which fails to recognize a zero polynomial and keeps iterating. I prepared a patch in which I extended RR domain to support eps (tolerance) parameter and robustified low-level polynomial arithmetic functions to take advantage of this, but this doesn't solve all problems with floating point numbers (e.g. Groebner bases algorithm over RR doesn't work as expected). It seems that the right way to go is (as suggested by Bartosz) to replace mpf dtype of RR with a class that can manage tolerance on its own. I'm now experimenting with mpf, because it may already have all features that we need to implement this.

CAUTION: you used an X instead of an x in the last equation.

--

You received this message because you are subscribed to the Google Groups "sympy" group.
To post to this group, send email to sy...@googlegroups.com.
To unsubscribe from this group, send email to sympy+un...@googlegroups.com.
For more options, visit this group at http://groups.google.com/group/sympy?hl=en.

Mateusz