Add regression test

Closes sympy/sympy#26497

diofant/tests/polys/test_polytools.py

@@ -3389,3 +3389,23 @@ def test_sympyissue_25723():
         assert gcd(x - i, (x - i)*(2*x + 1)) == x - i
 
     assert gcd(x - i, (x - i)*(2*x - 1)) == x - i
+
+
+def test_sympyissue_26497():
+    a_n, m = symbols('a_n m', real=True, positive=True)
+    expr = (-a_n**4*m**2/(a_n**2 + 1) - 2*a_n**4*m/(a_n**2 + 1) -
+            a_n**4 - a_n**4/(a_n**2 + 1) + 2*I*a_n**3*m**2/(a_n**2 + 1) +
+            4*I*a_n**3*m/(a_n**2 + 1) + 2*I*a_n**3 +
+            2*I*a_n**3/(a_n**2 + 1) + 2*I*a_n*m**2/(a_n**2 + 1) +
+            4*I*a_n*m/(a_n**2 + 1) + 2*I*a_n + 2*I*a_n/(a_n**2 + 1) +
+            m**2/(a_n**2 + 1) + 2*m/(a_n**2 + 1) + 1 + 1/(a_n**2 + 1))
+    expr2 = (-a_n**2 + 2*I*a_n + 1)*(a_n**2 + m**2 + 2*m + 2)
+
+    assert expr.factor() == expr2
+    assert expr.subs({a_n: 1, m: 1}) == expr2.subs({a_n: 1, m: 1}) == 12*I
+
+    expr = ((b - I)**2*(b + I)*(a + 1))
+    expr2 = expr.expand()
+
+    assert expr2.factor() == (a + 1)*(b - I)*(b**2 + 1)
+    assert expr2.factor(extension=True) == expr

docs/release/notes-0.15.rst

@@ -95,3 +95,4 @@ These Sympy issues also were addressed:
 * :sympyissue:`26343`: TypeError: Invalid NaN Comparison using dsolve for ODE with ics={v(0) : 0}
 * :sympyissue:`26313`: Error result for limit of a piecewise
 * :sympyissue:`26477`: Error in integral result using hyper
+* :sympyissue:`26497`: factor produces wrong output