Integration of certain cubic rational functions is incorrect

[ethankward]

>>> from sympy import *
>>> import mpmath
>>> x = Symbol('x')
>>> a, b = 2, 3
>>> n1 = N(integrate(1/(-28*x**3 - 46*x**2 - 25*x - 10), [x, a, b]))
>>> n2 = mpmath.quad(lambda x: 1/(-28*x**3 - 46*x**2 - 25*x - 10), [a, b])
>>> print(abs(n1 - n2))
0.00836361801336797

More examples: 1/(40*x**3 + 34*x**2 - 30*x - 28), 1/(-(x**3) - 15*x**2 + 17*x - 18), 1/(-(x**3) - 15*x**2 + 23*x + 40).

I've verified the numerical result with other tools, so I'm pretty sure the error isn't in mpmath.

[skirpichev]

Apparently, issue is not in the definite integration:

In [3]: r = integrate(1/(-28*x**3 - 46*x**2 - 25*x - 10), x)

In [4]: n, d = r.diff(x).as_numer_denom()

In [5]: n.is_number
Out[5]: True

In [6]: Poly(d, x).degree()
Out[6]: 1

[jksuom]

The antiderivative should be a sum of three logarithms since the denominator has three simple roots. It seems that log_to_real fails to recognize the two complex roots because 'chop' does not work as expected on this line:

(Pdb) p D.evalf(chop=True)
-4.63871372897167  # should be 0
(Pdb) p D.evalf()
-0.e-124  # even this is better

[asmeurer]

I would be careful about using chop=True in library code. I think the option only works when you know the value you are working with and know that it should go to 0. As far as I understand, it uses a method based on absolute value, which can easily lead to false positives if you actually do have a small but nonzero value (e.g., (1/2**S(100)).evalf(chop=True)).

[smichr]

-0.e-124 # even this is better

That looks like a 0 with no precision -- and like using chop we must not make decisions based upon precisionless evaluations.

[skirpichev]

Probably, best you can do here to test for zero - call evalf with strict=True and "pass" if it throw a PrecisionExhausted exception.

Similar heuristics I use in diofant/diofant#430. We can exactly solve zero-decision problem for algebraic numbers, but that seems to slow for me, at least with current polys module.