Repository: sympy/sympy

Issue:
Bug in indefinite integrals
I have found a curious bug in `integrate`. It may or may not be related to #28596. I say "may not" because, back to 1.8, I can't find any version of sympy that does it right. Just what it exactly does differs from version to version.

Below I calculate the so-called noise integrals of second-order lowpass and bandpass filter functions. The correct value for both is $\frac{\omega_p q_p}{4}$.

The integration with symbolic $\omega_p$ and $q_p$ gives an incorrect result.

For numerical values, it depends on the actual value whether the result is correct or not. Interestingly, increasing the parameter $q_p$ lets it go from correct to wrong, but it does not happen at the same value of $q_p$ for both funktions even though the two functions have the same denominator. 

Here is my test code:

```python
from sympy import *

s = symbols('s')
f = symbols('f', real=True)
wp, qp = symbols('omega_p, q_p', positive=True)
tf = [
    wp**2 / (s**2 + wp / qp * s + wp**2),
    wp*s / (s**2 + wp / qp * s + wp**2)]
print(pretty(tf))
b2 = [
    (abs(t.subs(s, I*2*pi*f))**2).simplify()
    for t in tf
]
print(pretty(b2))
for q in [1, 2, 10, 100]:
    expected = S(1)*q/4
    ni = [
        integrate(b.subs(wp, 1).subs(qp, q), (f, 0, oo))
        for b in b2
    ]
    print(pretty([expected, ni]))
```

This is the last part of the output:

```
>python int_bug.py
[1/4, [1/4, 1/4]]

[1/2, [1/2, 1/2]]

⎡     ⎡     5⋅(-2⋅√160797⋅π + 798⋅π) ⎤⎤
⎢5/2, ⎢5/2, ────────────────────────⎥⎥
⎣     ⎣              1596⋅π         ⎦⎦

⎡    ⎡     25⋅√1600079997  25⋅(-2⋅√1600079997⋅π + 79998⋅π)  ⎤⎤
⎢25, ⎢25 + ──────────────, ───────────────────────────────⎥⎥
⎣    ⎣         39999                   79998⋅π            ⎦⎦

```

Discussion hints:

Solution patch:
diff --git a/sympy/integrals/rationaltools.py b/sympy/integrals/rationaltools.py
index a37a587e7e62..8a33a4b954a7 100644
--- a/sympy/integrals/rationaltools.py
+++ b/sympy/integrals/rationaltools.py
@@ -340,16 +340,20 @@ def _roots_real_complex(poly):
     reals = {}
     complexes = {}
 
-    for r, m in list(rs.items()):
+    remaining = list(rs)
+
+    while remaining:
+        r = remaining.pop()
         r_c = r.conjugate()
         if r != r_c and r_c in rs:
-            assert rs.pop(r_c) == m
+            assert rs[r_c] == rs[r]
+            remaining.remove(r_c)
             r_re, r_im = r.as_real_imag()
             _, r_im = ordered([r_im, -r_im])
-            complexes[(r_re, r_im)] = m
+            complexes[(r_re, r_im)] = rs[r]
         else:
             # If we didn't find a conjugate just treat as real
-            reals[r] = m
+            reals[r] = rs[r]
 
     return reals, complexes
 


Test patch:
diff --git a/sympy/integrals/tests/test_rationaltools.py b/sympy/integrals/tests/test_rationaltools.py
index 9c8d4a3cc5c1..2380dca9aae8 100644
--- a/sympy/integrals/tests/test_rationaltools.py
+++ b/sympy/integrals/tests/test_rationaltools.py
@@ -1,4 +1,4 @@
-from sympy.core.numbers import (I, Rational)
+from sympy.core.numbers import (I, Rational, pi)
 from sympy.core.singleton import S
 from sympy.core.symbol import (Dummy, symbols)
 from sympy.functions.elementary.exponential import log
@@ -201,3 +201,9 @@ def test_issue_28186():
     )
     assert integrate(f, x) == F
     assert (F.diff(x) - f).ratsimp() == 0
+
+    # https://github.com/sympy/sympy/issues/28657
+    p = 400*pi**2*x**2/(1600*pi**4*x**4 - 796*pi**2*x**2 + 100)
+    P = ratint(p, x)
+    res = P.evalf(subs={x:1e100}) - P.evalf(subs={x:0})
+    assert abs(res - 2.5) < 1e-10