Repository: sympy/sympy

Issue:
Wrong result in simplify(1.0**(x+1)/1.0**x)
If I write:
```python
In [1]: simplify(1.0**(x+1)/1.0**x)
Out[1]: 
   2⋅x + 3
1.0    
```

I expect the result to be `1.0`.
It works as expected with any other numbers like `2.0`, `3.0`, etc...
If I write `simplify(1**(x+1)/1**x)` it works correctly.

Discussion:

Patch:
diff --git a/.mailmap b/.mailmap
index a14247066b..7116e10e89 100644
--- a/.mailmap
+++ b/.mailmap
@@ -770,6 +770,7 @@ Jason Ross <jasonross1024@gmail.com>
 Jason Siefken <siefkenj@gmail.com>
 Jason Tokayer <jason.tokayer@gmail.com>
 Jason Tokayer <jason.tokayer@gmail.com> <jason.tokayer@capitalone.com>
+Jatin Bhardwaj <bhardwajjatin093@gmail.com> <148186488+Jatinbhardwaj-093@users.noreply.github.com>
 Jatin Yadav <jatinyadav25@gmail.com>
 Javed Nissar <javednissar@gmail.com>
 Jay Patankar <patankarjays@gmail.com> Jay-Patankar <107458263+Jay-Patankar@users.noreply.github.com>
@@ -923,7 +924,6 @@ Louis Abraham <louis.abraham@yahoo.fr> louisabraham <louis.abraham@yahoo.fr>
 Luca Weihs <astronomicalcuriosity@gmail.com>
 Lucas Gallindo <lgallindo@gmail.com>
 Lucas Jones <lucas@lucasjones.co.uk>
-Lucas Kletzander <lucas.kletzander@gmail.com>
 Lucas Wiman <lucas.wiman@gmail.com>
 Lucy Mountain <lucymountain1@icloud.com> LucyMountain <lucymountain1@icloud.com>
 Luis Garcia <ppn.online@me.com>
diff --git a/doc/src/modules/solvers/solvers.rst b/doc/src/modules/solvers/solvers.rst
index 012bb7e119..7e7661586d 100644
--- a/doc/src/modules/solvers/solvers.rst
+++ b/doc/src/modules/solvers/solvers.rst
@@ -124,8 +124,3 @@ Linear Programming (Optimization)
 .. autofunction:: lpmin
 
 .. autofunction:: linprog
-
-Types
------
-
-.. py:class:: collections.abc.Sequence
diff --git a/sympy/concrete/products.py b/sympy/concrete/products.py
index e6a51f105a..c23c7212d7 100644
--- a/sympy/concrete/products.py
+++ b/sympy/concrete/products.py
@@ -285,6 +285,11 @@ def doit(self, **hints):
         else:
             return powsimp(f)
 
+    def _eval_adjoint(self):
+        if self.is_commutative:
+            return self.func(self.function.adjoint(), *self.limits)
+        return None
+
     def _eval_conjugate(self):
         return self.func(self.function.conjugate(), *self.limits)
 
diff --git a/sympy/core/basic.py b/sympy/core/basic.py
index 69334d9b01..e56118041a 100644
--- a/sympy/core/basic.py
+++ b/sympy/core/basic.py
@@ -1,7 +1,6 @@
 """Base class for all the objects in SymPy"""
 from __future__ import annotations
 
-from collections import Counter
 from collections.abc import Mapping, Iterable
 from itertools import zip_longest
 from functools import cmp_to_key
@@ -1804,7 +1803,16 @@ def find(self, query, group=False):
 
         if not group:
             return set(results)
-        return dict(Counter(results))
+        else:
+            groups = {}
+
+            for result in results:
+                if result in groups:
+                    groups[result] += 1
+                else:
+                    groups[result] = 1
+
+            return groups
 
     def count(self, query):
         """Count the number of matching subexpressions."""
diff --git a/sympy/core/expr.py b/sympy/core/expr.py
index 7e78484c2f..4f4d63ed74 100644
--- a/sympy/core/expr.py
+++ b/sympy/core/expr.py
@@ -1059,7 +1059,7 @@ def dir(self, x, cdir):
 
     def _eval_transpose(self):
         from sympy.functions.elementary.complexes import conjugate
-        if self.is_commutative:
+        if (self.is_complex or self.is_infinite):
             return self
         elif self.is_hermitian:
             return conjugate(self)
diff --git a/sympy/matrices/expressions/blockmatrix.py b/sympy/matrices/expressions/blockmatrix.py
index 0125d6233b..4b5b840010 100644
--- a/sympy/matrices/expressions/blockmatrix.py
+++ b/sympy/matrices/expressions/blockmatrix.py
@@ -1,6 +1,7 @@
 from sympy.assumptions.ask import (Q, ask)
 from sympy.core import Basic, Add, Mul, S
 from sympy.core.sympify import _sympify
+from sympy.functions import adjoint
 from sympy.functions.elementary.complexes import re, im
 from sympy.strategies import typed, exhaust, condition, do_one, unpack
 from sympy.strategies.traverse import bottom_up
@@ -186,9 +187,13 @@ def _eval_transpose(self):
         return BlockMatrix(M)
 
     def _eval_adjoint(self):
-        return BlockMatrix(
-            Matrix(self.blockshape[0], self.blockshape[1], self.blocks).adjoint()
-        )
+        # Adjoint all the individual matrices
+        matrices = [adjoint(matrix) for matrix in self.blocks]
+        # Make a copy
+        M = Matrix(self.blockshape[0], self.blockshape[1], matrices)
+        # Transpose the block structure
+        M = M.transpose()
+        return BlockMatrix(M)
 
     def _eval_trace(self):
         if self.rowblocksizes == self.colblocksizes:
diff --git a/sympy/matrices/matrixbase.py b/sympy/matrices/matrixbase.py
index f994f362bf..442f8c37fe 100644
--- a/sympy/matrices/matrixbase.py
+++ b/sympy/matrices/matrixbase.py
@@ -2029,7 +2029,7 @@ def iter_items(self):
         return self._eval_iter_items()
 
     def _eval_adjoint(self):
-        return self.transpose().applyfunc(lambda x: x.adjoint())
+        return self.transpose().conjugate()
 
     def _eval_applyfunc(self, f):
         cols = self.cols
@@ -2191,7 +2191,7 @@ def H(self):
         conjugate: By-element conjugation
         sympy.matrices.matrixbase.MatrixBase.D: Dirac conjugation
         """
-        return self.adjoint()
+        return self.T.C
 
     def permute(self, perm, orientation='rows', direction='forward'):
         r"""Permute the rows or columns of a matrix by the given list of
diff --git a/sympy/matrices/solvers.py b/sympy/matrices/solvers.py
index 1fba990df8..01a85b29be 100644
--- a/sympy/matrices/solvers.py
+++ b/sympy/matrices/solvers.py
@@ -357,7 +357,7 @@ def _LUsolve(M, rhs, iszerofunc=_iszero):
     for i in range(m):
         for j in range(min(i, n)):
             scale = A[i, j]
-            b.zip_row_op(i, j, lambda x, y: dps(x - scale * y))
+            b.zip_row_op(i, j, lambda x, y: dps(x - y * scale))
 
     # consistency check for overdetermined systems
     if m > n:
@@ -372,10 +372,10 @@ def _LUsolve(M, rhs, iszerofunc=_iszero):
     for i in range(n - 1, -1, -1):
         for j in range(i + 1, n):
             scale = A[i, j]
-            b.zip_row_op(i, j, lambda x, y: dps(x - scale * y))
+            b.zip_row_op(i, j, lambda x, y: dps(x - y * scale))
 
         scale = A[i, i]
-        b.row_op(i, lambda x, _: dps(scale**-1 * x))
+        b.row_op(i, lambda x, _: dps(x / scale))
 
     return rhs.__class__(b)
 
diff --git a/sympy/polys/polytools.py b/sympy/polys/polytools.py
index 807215733b..bd4f703f5b 100644
--- a/sympy/polys/polytools.py
+++ b/sympy/polys/polytools.py
@@ -164,9 +164,6 @@ class Poly(Basic):
     is_Poly = True
     _op_priority = 10.001
 
-    rep: DMP
-    gens: tuple[Expr, ...]
-
     def __new__(cls, rep, *gens, **args) -> Poly:
         """Create a new polynomial instance out of something useful. """
         opt = options.build_options(gens, args)
@@ -3360,7 +3357,7 @@ def sqf_list_include(f, all=False):
 
         return [(f.per(g), k) for g, k in factors]
 
-    def factor_list(f) -> tuple[Expr, list[tuple[Poly, int]]]:
+    def factor_list(f):
         """
         Returns a list of irreducible factors of ``f``.
 
@@ -6264,12 +6261,6 @@ def sqf_part(f, *gens, **args):
         return result
 
 
-def _poly_sort_key(poly):
-    """Sort a list of polys."""
-    rep = poly.rep.to_list()
-    return (len(rep), len(poly.gens), str(poly.domain), rep)
-
-
 def _sorted_factors(factors, method):
     """Sort a list of ``(expr, exp)`` pairs. """
     if method == 'sqf':
diff --git a/sympy/printing/pycode.py b/sympy/printing/pycode.py
index f07785d612..9d768b8c0a 100644
--- a/sympy/printing/pycode.py
+++ b/sympy/printing/pycode.py
@@ -557,7 +557,7 @@ def _print_sign(self, e):
 
     def _print_Not(self, expr):
         PREC = precedence(expr)
-        return self._operators['not'] + ' ' + self.parenthesize(expr.args[0], PREC)
+        return self._operators['not'] + self.parenthesize(expr.args[0], PREC)
 
     def _print_IndexedBase(self, expr):
         return expr.name
diff --git a/sympy/simplify/powsimp.py b/sympy/simplify/powsimp.py
index 51baa56c97..f72dfeb072 100644
--- a/sympy/simplify/powsimp.py
+++ b/sympy/simplify/powsimp.py
@@ -5,7 +5,7 @@
 from sympy.core.function import expand_log, count_ops, _coeff_isneg
 from sympy.core import sympify, Basic, Dummy, S, Add, Mul, Pow, expand_mul, factor_terms
 from sympy.core.sorting import ordered, default_sort_key
-from sympy.core.numbers import Integer, Rational
+from sympy.core.numbers import Integer, Rational, equal_valued
 from sympy.core.mul import _keep_coeff
 from sympy.core.rules import Transform
 from sympy.functions import exp_polar, exp, log, root, polarify, unpolarify
@@ -188,6 +188,10 @@ def recurse(arg, **kwargs):
             bpos = b.is_positive or b.is_polar
             if bpos:
                 binv = 1/b
+                #Special case for float 1
+                if b.is_Float and equal_valued(b, 1):
+                    c_powers[b] = S.One
+                    continue
                 if b != binv and binv in c_powers:
                     if b.as_numer_denom()[0] is S.One:
                         c_powers.pop(b)
diff --git a/sympy/solvers/ode/nonhomogeneous.py b/sympy/solvers/ode/nonhomogeneous.py
index ae39d55664..715de97000 100644
--- a/sympy/solvers/ode/nonhomogeneous.py
+++ b/sympy/solvers/ode/nonhomogeneous.py
@@ -8,7 +8,7 @@
 instances in other classes for using them it is better to keep it here as separate helpers.
 
 """
-from collections import Counter
+from collections import defaultdict
 from sympy.core import Add, S
 from sympy.core.function import diff, expand, _mexpand, expand_mul
 from sympy.core.relational import Eq
@@ -83,7 +83,9 @@ def _get_euler_characteristic_eq_sols(eq, func, match_obj):
     constants.reverse()
 
     # Create a dict root: multiplicity or charroots
-    charroots = Counter(chareqroots)
+    charroots = defaultdict(int)
+    for root in chareqroots:
+        charroots[root] += 1
     gsol = S.Zero
     ln = log
     for root, multiplicity in charroots.items():
@@ -202,7 +204,9 @@ def _get_const_characteristic_eq_sols(r, func, order):
     chareq_is_complex = not all(i.is_real for i in chareq.all_coeffs())
 
     # Create a dict root: multiplicity or charroots
-    charroots = Counter(chareqroots)
+    charroots = defaultdict(int)
+    for root in chareqroots:
+        charroots[root] += 1
     # We need to keep track of terms so we can run collect() at the end.
     # This is necessary for constantsimp to work properly.
     collectterms = []
diff --git a/sympy/solvers/polysys.py b/sympy/solvers/polysys.py
index 2e8340b62d..83026664c8 100644
--- a/sympy/solvers/polysys.py
+++ b/sympy/solvers/polysys.py
@@ -1,21 +1,12 @@
 """Solvers of systems of polynomial equations. """
-
-from __future__ import annotations
-
-from typing import Any
-from collections.abc import Sequence, Iterable
-
 import itertools
+from collections import defaultdict
 
 from sympy import Dummy
 from sympy.core import S
-from sympy.core.expr import Expr
-from sympy.core.exprtools import factor_terms
 from sympy.core.sorting import default_sort_key
-from sympy.logic.boolalg import Boolean
 from sympy.polys import Poly, groebner, roots
-from sympy.polys.domains import ZZ
-from sympy.polys.polytools import parallel_poly_from_expr, sqf_part
+from sympy.polys.polytools import parallel_poly_from_expr
 from sympy.polys.polyerrors import (
     ComputationFailed,
     PolificationFailed,
@@ -25,8 +16,8 @@
 )
 from sympy.simplify import rcollect
 from sympy.utilities import postfixes
-from sympy.utilities.iterables import cartes
 from sympy.utilities.misc import filldedent
+from sympy.utilities.iterables import cartes
 from sympy.logic.boolalg import Or, And
 from sympy.core.relational import Eq
 
@@ -450,7 +441,7 @@ def solve_triangulated(polys, *gens, **args):
     return sorted((s for s, _ in solutions), key=default_sort_key)
 
 
-def factor_system(eqs: Sequence[Expr | complex], gens: Sequence[Expr] = (), **kwargs: Any) -> list[list[Expr]]:
+def factor_system(eqs, *gens, **kwargs):
     """
     Factorizes a system of polynomial equations into
     irreducible subsystems.
@@ -462,7 +453,7 @@ def factor_system(eqs: Sequence[Expr | complex], gens: Sequence[Expr] = (), **kw
         List of expressions to be factored.
         Each expression is assumed to be equal to zero.
 
-    gens : list, optional
+    *gens : Symbol or sequence of Symbols, optional
         Generator(s) of the polynomial ring.
         If not provided, all free symbols will be used.
 
@@ -472,80 +463,79 @@ def factor_system(eqs: Sequence[Expr | complex], gens: Sequence[Expr] = (), **kw
     Returns
     =======
 
-    list[list[Expr]]
-        A list of lists of expressions, where each sublist represents
-        an irreducible subsystem. When solved, each subsystem gives
-        one component of the solution. Only generic solutions are
-        returned (cases not requiring parameters to be zero).
+    tuple
+        A pair (systems, condition) where:
+
+        * systems is a list of lists of expressions,
+          where each sublist when solved gives one
+          component of the solution
+
+        * condition is a Boolean expression that must
+          be satisfied for the system to be solvable
 
     Examples
     ========
 
-    >>> from sympy.solvers.polysys import factor_system, factor_system_cond
+    >>> from sympy.solvers.polysys import factor_system
     >>> from sympy.abc import x, y, a, b, c
 
     A simple system with multiple solutions:
 
-    >>> factor_system([x**2 - 1, y - 1])
-    [[x + 1, y - 1], [x - 1, y - 1]]
+    >>> systems, cond = factor_system([x**2 - 1, y - 1])
+    >>> systems
+    [[x - 1, y - 1], [x + 1, y - 1]]
+    >>> cond
+    True
 
     A system with no solution:
 
-    >>> factor_system([x, 1])
+    >>> systems, cond = factor_system([x, 1])
+    >>> systems
     []
+    >>> cond
+    True
 
-    A system where any value of the symbol(s) is a solution:
-
-    >>> factor_system([x - x, (x + 1)**2 - (x**2 + 2*x + 1)])
-    [[]]
-
-    A system with no generic solution:
-
-    >>> factor_system([a*x*(x-1), b*y, c], [x, y])
-    []
-
-    If c is added to the unknowns then the system has a generic solution:
-
-    >>> factor_system([a*x*(x-1), b*y, c], [x, y, c])
-    [[x - 1, y, c], [x, y, c]]
+    >>> systems, cond = factor_system([a*x*(x-1), b*y, c], [x, y])
+    >>> systems
+    [[x, y], [x - 1, y]]
+    >>> cond
+    Eq(c, 0)
 
-    Alternatively :func:`factor_system_cond` can be used to get degenerate
-    cases as well:
+    The above cond signifies that the symbol c must be equal to
+    zero for the system to be solvable.
 
-    >>> factor_system_cond([a*x*(x-1), b*y, c], [x, y])
-    [[x - 1, y, c], [x, y, c], [x - 1, b, c], [x, b, c], [y, a, c], [a, b, c]]
+    A system with infinite solutions:
 
-    Each of the above cases is only satisfiable in the degenerate case `c = 0`.
+    >>> systems, cond = factor_system([x - x, (x + 1)**2 - (x**2 + 2*x + 1)])
+    >>> systems
+    [[]]
+    >>> cond
+    True
 
     The solution set of the original system represented
     by eqs is the union of the solution sets of the
     factorized systems.
 
-    An empty list [] means no generic solution exists.
-    A list containing an empty list [[]] means any value of
-    the symbol(s) is a solution.
+    A return of a list containing an empty list [[]]
+    in the tuple means any value of the symbol(s)
+    is a solution whereas a return of an empty list []
+    in the tuple means no solutions for the system exists.
 
     See Also
     ========
 
-    factor_system_cond : Returns both generic and degenerate solutions
-    factor_system_bool : Returns a Boolean combination representing all solutions
+    factor_system_cond : Returns both factors and conditions
+    factor_system_bool : Returns a Boolean combination representing the solution
     sympy.polys.polytools.factor : Factors a polynomial into irreducible factors
                                    over the rational numbers
     """
 
-    systems = _factor_system_poly_from_expr(eqs, gens, **kwargs)
-    systems_generic = [sys for sys in systems if not _is_degenerate(sys)]
-    systems_expr = [[p.as_expr() for p in system] for system in systems_generic]
-    return systems_expr
+    systems, conds = factor_system_cond(eqs, *gens, **kwargs)
+    systems = [[f for f, c in system] for system in systems]
+    return systems, And(*[Eq(c, 0) for c in conds]) if conds else True
 
 
-def _is_degenerate(system: list[Poly]) -> bool:
-    """Helper function to check if a system is degenerate"""
-    return any(p.is_ground for p in system)
-
-
-def factor_system_bool(eqs: Sequence[Expr | complex], gens: Sequence[Expr] = (), **kwargs: Any) -> Boolean:
+def factor_system_bool(eqs, *gens, **kwargs):
     """
     Factorizes a system of polynomial equations into irreducible DNF.
 
@@ -560,7 +550,7 @@ def factor_system_bool(eqs: Sequence[Expr | complex], gens: Sequence[Expr] = (),
        List of expressions to be factored.
        Each expression is assumed to be equal to zero.
 
-    gens : list, optional
+    *gens : Symbol or sequence of Symbols, optional
        Generator(s) of the polynomial ring.
        If not provided, all free symbols will be used.
 
@@ -588,11 +578,12 @@ def factor_system_bool(eqs: Sequence[Expr | complex], gens: Sequence[Expr] = (),
     (Eq(x - 1, 0) & Eq(y - 1, 0)) | (Eq(x + 1, 0) & Eq(y - 1, 0))
 
     >>> eqs = [a * (x - 1), b]
-    >>> factor_system_bool([a*(x - 1), b])
-    (Eq(a, 0) & Eq(b, 0)) | (Eq(b, 0) & Eq(x - 1, 0))
+    >>> factor_system_bool(eqs, [x])
+    Eq(b, 0) & (Eq(a, 0) | Eq(x - 1, 0))
 
-    >>> factor_system_bool([a*x**2 - a, b*(x + 1), c], [x])
-    (Eq(c, 0) & Eq(x + 1, 0)) | (Eq(a, 0) & Eq(b, 0) & Eq(c, 0)) | (Eq(b, 0) & Eq(c, 0) & Eq(x - 1, 0))
+    >>> eqs = [a * x ** 2 - a, b * (x + 1), c]
+    >>> factor_system_bool(eqs, [x])
+    Eq(c, 0) & (Eq(x + 1, 0) | (Eq(a, 0) & Eq(b, 0)))
 
     >>> factor_system_bool([x**2 + 2*x + 1 - (x + 1)**2])
     True
@@ -610,14 +601,26 @@ def factor_system_bool(eqs: Sequence[Expr | complex], gens: Sequence[Expr] = (),
 
     """
 
-    systems = factor_system_cond(eqs, gens, **kwargs)
-    return Or(*[And(*[Eq(eq, 0) for eq in sys]) for sys in systems])
+    systems, conds = factor_system_cond(eqs, *gens, **kwargs)
+    sys = Or(*[And(*[_eq2bool(eq) for eq in sys]) for sys in systems])
+    if conds:
+        sys &= And(*[Eq(c, 0) for c in conds])
+    return sys
+
+
+def _eq2bool(eq):
+    """Convert a (factor, conditions) pair to a Boolean equation."""
+    f, cs = eq
+    b = Eq(f, 0)
+    if cs:
+        b |= And(*[Eq(c, 0) for c in cs])
+    return b
 
 
-def factor_system_cond(eqs: Sequence[Expr | complex], gens: Sequence[Expr] = (), **kwargs: Any) -> list[list[Expr]]:
+def factor_system_cond(eqs, *gens, **kwargs):
     """
     Factorizes a polynomial system into irreducible components and returns
-    both generic and degenerate solutions.
+    factors with their associated conditions.
 
     Parameters
     ==========
@@ -626,7 +629,7 @@ def factor_system_cond(eqs: Sequence[Expr | complex], gens: Sequence[Expr] = (),
         List of expressions to be factored.
         Each expression is assumed to be equal to zero.
 
-    gens : list, optional
+    *gens : Symbol or sequence of Symbols, optional
         Generator(s) of the polynomial ring.
         If not provided, all free symbols will be used.
 
@@ -636,10 +639,15 @@ def factor_system_cond(eqs: Sequence[Expr | complex], gens: Sequence[Expr] = (),
     Returns
     =======
 
-    list[list[Expr]]
-        A list of lists of expressions, where each sublist represents
-        an irreducible subsystem. Includes both generic solutions and
-        degenerate cases requiring equality conditions on parameters.
+    tuple
+        A pair (systems, conditions) where:
+
+        * systems is a list of lists of (factor, conditions) pairs,
+          where factor is an expression and conditions is a list
+          of expressions that must be zero for the factor to be relevant
+
+        * conditions is a list of expressions that must be zero
+          for any solution to exist
 
     Examples
     ========
@@ -647,39 +655,30 @@ def factor_system_cond(eqs: Sequence[Expr | complex], gens: Sequence[Expr] = (),
     >>> from sympy.solvers.polysys import factor_system_cond
     >>> from sympy.abc import x, y, a, b, c
 
-    >>> factor_system_cond([x**2 - 4, a*y, b], [x, y])
-    [[x + 2, y, b], [x - 2, y, b], [x + 2, a, b], [x - 2, a, b]]
+    >>> systems, conds = factor_system_cond([x**2 - 4, a*y, b], [x, y])
+    >>> systems
+    [[(x - 2, []), (y, [a])], [(x + 2, []), (y, [a])]]
+    >>> conds
+    [b]
 
-    >>> factor_system_cond([a*x*(x-1), b*y, c], [x, y])
-    [[x - 1, y, c], [x, y, c], [x - 1, b, c], [x, b, c], [y, a, c], [a, b, c]]
+    >>> systems, conds = factor_system_cond([a*x*(x-1), b*y, c], [x, y])
+    >>> systems
+    [[(x, [a]), (y, [b])], [(x - 1, [a]), (y, [b])]]
+    >>> conds
+    [c]
 
-    An empty list [] means no solution exists.
-    A list containing an empty list [[]] means any value of
-    the symbol(s) is a solution.
+    In the return of the (systems, conds) pair, systems is [[]]
+    when the system implies tautology. Eg. x - x = 0,
+    and [] when the system of equations is unatisfiable. Eg. 1 = 0
 
     See Also
     ========
 
-    factor_system : Returns only generic solutions
-    factor_system_bool : Returns a Boolean combination representing all solutions
+    factor_system : Returns factors and solvability condition separately
+    factor_system_bool : Returns a Boolean combination representing the solution
     sympy.polys.polytools.factor : Factors a polynomial into irreducible factors
                                    over the rational numbers
     """
-    systems_poly = _factor_system_poly_from_expr(eqs, gens, **kwargs)
-    systems = [[p.as_expr() for p in system] for system in systems_poly]
-    return systems
-
-
-def _factor_system_poly_from_expr(
-        eqs: Sequence[Expr | complex], gens: Sequence[Expr], **kwargs: Any
-) -> list[list[Poly]]:
-    """
-    Convert expressions to polynomials and factor the system.
-
-    Takes a sequence of expressions, converts them to
-    polynomials, and factors the resulting system. Handles both regular
-    polynomial systems and purely numerical cases.
-    """
     try:
         polys, opts = parallel_poly_from_expr(eqs, *gens, **kwargs)
         only_numbers = False
@@ -690,30 +689,43 @@ def _factor_system_poly_from_expr(
         only_numbers = True
 
     if only_numbers:
-        return [[]] if all(p == 0 for p in polys) else []
+        if all(p == 0 for p in polys):
+            systems = [[]]
+            conditions = []
+        else:
+            systems = []
+            conditions = []
+    else:
+        systems, conditions = factor_system_poly(polys)
+        systems = [[(f.as_expr(), [c.as_expr() for c in cs]) for f, cs in system] for system in systems]
+        conditions = [c.as_expr() for c in conditions]
 
-    return factor_system_poly(polys)
+    return systems, conditions
 
 
-def factor_system_poly(polys: list[Poly]) -> list[list[Poly]]:
+def factor_system_poly(polys):
     """
-    Factors a system of polynomial equations into irreducible subsystems
-
+    Factors a system of polynomials into irreducible factors with conditions.
     Core implementation that works directly with Poly instances.
 
     Parameters
     ==========
 
-    polys : list[Poly]
+    polys : list
         A list of Poly instances to be factored.
 
     Returns
     =======
 
-    list[list[Poly]]
-        A list of lists of polynomials, where each sublist represents
-        an irreducible component of the solution. Includes both
-        generic and degenerate cases.
+    tuple
+        A pair (systems, constant_conds) where:
+
+        * systems is a list of lists of (factor, conditions) pairs, with each
+          factor being a Poly instance and conditions being a tuple of Poly
+          instances
+
+        * constant_conds is a list of Poly instances representing conditions that
+          must be satisfied for any solution to exist
 
     Examples
     ========
@@ -729,120 +741,144 @@ def factor_system_poly(polys: list[Poly]) -> list[list[Poly]]:
     of the two conditions on the parameters ``a`` and ``b`` is nonzero
     and the constant parameter ``c`` should be zero.
 
-    >>> sys1, sys2 = factor_system_poly([p1, p2, p3])
-    >>> sys1
-    [Poly(x - 2, x, domain='ZZ[a,b,c]'),
-     Poly(c, x, domain='ZZ[a,b,c]')]
-    >>> sys2
-    [Poly(a - 1, x, domain='ZZ[a,b,c]'),
-     Poly(b - 3, x, domain='ZZ[a,b,c]'),
-     Poly(c, x, domain='ZZ[a,b,c]')]
+    >>> systems, constant_conds = factor_system_poly([p1, p2, p3])
+    >>> systems
+    [[(Poly(x - 2, x, domain='ZZ[a,b,c]'), (Poly(a - 1, x, domain='ZZ[a,b,c]'), Poly(b - 3, x, domain='ZZ[a,b,c]')))]]
+
+    >>> constant_conds
+    [Poly(c, x, domain='ZZ[a,b,c]')]
 
-     An empty list [] when returned means no solution exists.
-     Whereas a list containing an empty list [[]] means any value is a solution.
+    This is the core routine used by higher-level functions factor_system_cond,
+    factor_system, and factor_system_bool. It Returns empty systems list
+    ([]) when no solution exists and [[]] (list containing empty list) when any
+    value is a solution
 
     See Also
     ========
 
-    factor_system : Returns only generic solutions
-    factor_system_bool : Returns a Boolean combination representing the solutions
-    factor_system_cond : Returns both generic and degenerate solutions
+    factor_system : Returns factors and solvability condition separately
+    factor_system_bool : Returns a Boolean combination representing the solution
+    factor_system_cond : Returns both factors and conditions
     sympy.polys.polytools.factor : Factors a polynomial into irreducible factors
                                    over the rational numbers
     """
     if not all(isinstance(poly, Poly) for poly in polys):
         raise TypeError("polys should be a list of Poly instances")
     if not polys:
-        return [[]]
+        return [[]], []
+    else:
+        base_domain = polys[0].domain
+        base_gens = polys[0].gens
+        if not all(poly.domain == base_domain and poly.gens == base_gens for poly in polys[1:]):
+            raise DomainError("All polynomials must have the same domain and generators")
 
-    base_domain = polys[0].domain
-    base_gens = polys[0].gens
-    if not all(poly.domain == base_domain and poly.gens == base_gens for poly in polys[1:]):
-        raise DomainError("All polynomials must have the same domain and generators")
+    constant_eqs = []
+    eqs_factors = []
+    conds_factor = defaultdict(list)
 
-    factor_sets = []
     for poly in polys:
         constant, factors_mult = poly.factor_list()
 
-        if constant.is_zero is True:
-            continue
+        factors = [f for f, m in factors_mult]
+
+        if factors:
+            eqs_factors.append(factors)
+            if constant.is_zero is not False:
+                constp = constant.as_poly(base_gens, domain=base_domain)
+                for f in factors:
+                    if constp not in conds_factor[f]:
+                        conds_factor[f].append(constp)
+        elif constant.is_zero is True:
+            pass
         elif constant.is_zero is False:
-            if not factors_mult:
-                return []
-            factor_sets.append([f for f, _ in factors_mult])
+            return ([], [])
         else:
-            constant = sqf_part(factor_terms(constant).as_coeff_Mul()[1])
-            constp = Poly(constant, base_gens, domain=base_domain)
-            factors = [f for f, _ in factors_mult]
-            factors.append(constp)
-            factor_sets.append(factors)
+            constp = constant.as_poly(base_gens, domain=base_domain)
+            constant_eqs.append(constp)
 
-    if not factor_sets:
-        return [[]]
+    fac2conds = {f: tuple(conds) for f, conds in conds_factor.items()}
 
-    result = _factor_sets(factor_sets)
-    return _sort_systems(result)
+    cnf = [[(f, fac2conds.get(f, ())) for f in eq] for eq in eqs_factors]
 
+    dnf = _cnf2dnf(cnf)
 
-def _factor_sets_slow(eqs: list[list]) -> set[frozenset]:
-    """
-    Helper to find the minimal set of factorised subsystems that is
-    equivalent to the original system.
+    return dnf, constant_eqs
 
-    The result is in DNF.
-    """
-    if not eqs:
-        return {frozenset()}
-    systems_set = {frozenset(sys) for sys in cartes(*eqs)}
-    return {s1 for s1 in systems_set if not any(s1 > s2 for s2 in systems_set)}
+
+def _has_common_variables(poly1, poly2):
+    """Helper function to check if two polynomials have nonzero degree in any common variables"""
+    return any(d1 and d2 for d1, d2 in zip(poly1.degree_list(), poly2.degree_list()))
 
 
-def _factor_sets(eqs: list[list]) -> set[frozenset]:
+def _cnf2dnf(eqs):
     """
-    Helper that builds factor combinations.
+    Given a list of lists of (factor, conditions) pairs from the factorization of a
+    polynomial system, returns the minimal DNF sufficient to
+    satisfy the CNF. Only includes terms necessary for
+    satisfying the CNF, omitting redundant factors, hence different from a simple
+    mechanical rewrite of CNF to DNF.
+
+    The input is a list of lists of (Poly, tuple) pairs, where each Poly represents a
+    factor and the tuple contains conditions on that factor..
     """
-    if not eqs:
-        return {frozenset()}
-
-    current_set = min(eqs, key=len)
-    other_sets = [s for s in eqs if s is not current_set]
-
-    stack = [(factor, [s for s in other_sets if factor not in s], {factor})
-             for factor in current_set]
-
-    result = set()
-
-    while stack:
-        factor, remaining_sets, current_solution = stack.pop()
-
-        if not remaining_sets:
-            result.add(frozenset(current_solution))
-            continue
-
-        next_set = min(remaining_sets, key=len)
-        next_remaining = [s for s in remaining_sets if s is not next_set]
-
-        for next_factor in next_set:
-            valid_remaining = [s for s in next_remaining if next_factor not in s]
-            new_solution = current_solution | {next_factor}
-            stack.append((next_factor, valid_remaining, new_solution))
-
-    return {s1 for s1 in result if not any(s1 > s2 for s2 in result)}
-
-
-def _sort_systems(systems: Iterable[Iterable[Poly]]) -> list[list[Poly]]:
-    """Sorts a list of lists of polynomials"""
-    systems_list = [sorted(s, key=_poly_sort_key, reverse=True) for s in systems]
-    return sorted(systems_list, key=_sys_sort_key, reverse=True)
+    # remove equations that are independent of all others
+    # and generate the Cartesian product for them
+    eqs = [list(i) for i in eqs]
+    indep = set()
+
+    for i in range(len(eqs)):
+        is_independent = True
+        for j in range(len(eqs)):
+            if i != j:
+                for expr1 in eqs[i]:
+                    for expr2 in eqs[j]:
+                        if _has_common_variables(expr1[0], expr2[0]):
+                            is_independent = False
+                            break
+                    if not is_independent:
+                        break
+            if not is_independent:
+                break
+        if is_independent:
+            indep.add(i)
+
+    if not indep:
+        return list(_dnf(eqs))
+
+    indep_eqs = []
+    dep_eqs = []
+    for i, e in enumerate(eqs):
+        if i in indep:
+            indep_eqs.append(e)
+        else:
+            dep_eqs.append(e)
 
+    result = []
+    for t in cartes(*indep_eqs):
+        t_list = list(t)
+        for d in _dnf(dep_eqs):
+            result.append(t_list + d)
+    return result
 
-def _poly_sort_key(poly):
-    """Sort key for polynomials"""
-    if poly.domain.is_FF:
-        poly = poly.set_domain(ZZ)
-    return poly.degree_list(), poly.rep.to_list()
 
+def _dnf(eqs):
+    # helper for _cnf2dnf that recursively enumerates the
+    # minimal dnf args that satisfy the cnf expression;
+    if not eqs:
+        return [[]]
+    elif len(eqs) == 1:
+        return [[f] for f in eqs[0]]
+    else:
+        f = eqs[0][0]
+        eqs_f_zero = [eq for eq in eqs if f not in eq]
+        dnf_f = _dnf(eqs_f_zero)
+
+        f_free_eqs = [[x for x in eq if x != f] for eq in eqs]
+        if not all(f_free_eqs):
+            dnf = [[f] + s for s in dnf_f]
+        else:
+            eqs_f_nonzero = list(filter(None, f_free_eqs))
+            dnf_no_f = _dnf(eqs_f_nonzero)
+            dnf = dnf_no_f + [[f] + s for s in dnf_f if s not in dnf_no_f]
 
-def _sys_sort_key(sys):
-    """Sort key for lists of polynomials"""
-    return list(zip(*map(_poly_sort_key, sys)))
+        return dnf


Test patch:
diff --git a/sympy/functions/elementary/tests/test_complexes.py b/sympy/functions/elementary/tests/test_complexes.py
index 699c0fef96..be2a45cd9b 100644
--- a/sympy/functions/elementary/tests/test_complexes.py
+++ b/sympy/functions/elementary/tests/test_complexes.py
@@ -1,3 +1,4 @@
+from sympy.core.expr import Expr
 from sympy.core.function import (Derivative, Function, Lambda, expand, PoleError)
 from sympy.core.numbers import (E, I, Rational, comp, nan, oo, pi, zoo)
 from sympy.core.relational import Eq
@@ -694,11 +695,11 @@ def test_adjoint():
 
     x, y = symbols('x y')
     assert adjoint(adjoint(x)) == x
-    assert adjoint(x + y) == conjugate(x) + conjugate(y)
-    assert adjoint(x - y) == conjugate(x) - conjugate(y)
-    assert adjoint(x * y) == conjugate(x) * conjugate(y)
-    assert adjoint(x / y) == conjugate(x) / conjugate(y)
-    assert adjoint(-x) == -conjugate(x)
+    assert adjoint(x + y) == adjoint(x) + adjoint(y)
+    assert adjoint(x - y) == adjoint(x) - adjoint(y)
+    assert adjoint(x * y) == adjoint(x) * adjoint(y)
+    assert adjoint(x / y) == adjoint(x) / adjoint(y)
+    assert adjoint(-x) == -adjoint(x)
 
     x, y = symbols('x y', commutative=False)
     assert adjoint(adjoint(x)) == x
@@ -736,7 +737,7 @@ def test_conjugate():
 
 
 def test_conjugate_transpose():
-    x = Symbol('x', commutative=False)
+    x = Symbol('x')
     assert conjugate(transpose(x)) == adjoint(x)
     assert transpose(conjugate(x)) == adjoint(x)
     assert adjoint(transpose(x)) == conjugate(x)
@@ -744,7 +745,16 @@ def test_conjugate_transpose():
     assert adjoint(conjugate(x)) == transpose(x)
     assert conjugate(adjoint(x)) == transpose(x)
 
-    x = Symbol('x')
+    class Symmetric(Expr):
+        def _eval_adjoint(self):
+            return None
+
+        def _eval_conjugate(self):
+            return None
+
+        def _eval_transpose(self):
+            return self
+    x = Symmetric()
     assert conjugate(x) == adjoint(x)
     assert transpose(x) == x
 
@@ -756,11 +766,11 @@ def test_transpose():
 
     x, y = symbols('x y')
     assert transpose(transpose(x)) == x
-    assert transpose(x + y) == x + y
-    assert transpose(x - y) == x - y
-    assert transpose(x * y) == x * y
-    assert transpose(x / y) == x / y
-    assert transpose(-x) == -x
+    assert transpose(x + y) == transpose(x) + transpose(y)
+    assert transpose(x - y) == transpose(x) - transpose(y)
+    assert transpose(x * y) == transpose(x) * transpose(y)
+    assert transpose(x / y) == transpose(x) / transpose(y)
+    assert transpose(-x) == -transpose(x)
 
     x, y = symbols('x y', commutative=False)
     assert transpose(transpose(x)) == x
diff --git a/sympy/functions/elementary/tests/test_exponential.py b/sympy/functions/elementary/tests/test_exponential.py
index ee8c311d01..9ababee840 100644
--- a/sympy/functions/elementary/tests/test_exponential.py
+++ b/sympy/functions/elementary/tests/test_exponential.py
@@ -161,7 +161,6 @@ def test_exp_subs():
 
 
 def test_exp_adjoint():
-    x = Symbol('x', commutative=False)
     assert adjoint(exp(x)) == exp(adjoint(x))
 
 
diff --git a/sympy/matrices/tests/test_matrices.py b/sympy/matrices/tests/test_matrices.py
index 4fb69a9afe..a5a1bf86ee 100644
--- a/sympy/matrices/tests/test_matrices.py
+++ b/sympy/matrices/tests/test_matrices.py
@@ -2639,18 +2639,6 @@ def test_adjoint():
     for cls in classes:
         assert ans == cls(dat).adjoint()
 
-
-def test_adjoint_with_operator():
-    # Regression test for issue 25130: adjoint() should propagate to operators
-    import sympy.physics.quantum
-    a = sympy.physics.quantum.operator.Operator('a')
-    a_dag = sympy.physics.quantum.Dagger(a)
-    dat = [[0, I * a], [0, a_dag]]
-    ans = Matrix([[0, 0], [-I * a_dag, a]])
-    for cls in classes:
-        assert ans == cls(dat).adjoint()
-
-
 def test_simplify_immutable():
     assert simplify(ImmutableMatrix([[sin(x)**2 + cos(x)**2]])) == \
                     ImmutableMatrix([[1]])
diff --git a/sympy/matrices/tests/test_solvers.py b/sympy/matrices/tests/test_solvers.py
index c1347062c0..5f5152e462 100644
--- a/sympy/matrices/tests/test_solvers.py
+++ b/sympy/matrices/tests/test_solvers.py
@@ -121,20 +121,6 @@ def test_LUsolve():
     raises(NonInvertibleMatrixError, lambda: A.LUsolve(b))
 
 
-def test_LUsolve_noncommutative():
-    a0, a1, a2, a3 = symbols("a:4", commutative=False)
-    b0, b1 = symbols("b:2", commutative=False)
-    A = Matrix([[a0, a1], [a2, a3]])
-    check = A * A.LUsolve(Matrix([b0, b1]))
-    assert check[0, 0].expand() == b0
-    # Because sympy simplification is very limited with noncommutative expressions,
-    # perform an explicit check with the second element
-    assert check[1, 0] == (
-        a2*a0**(-1)*(-a1*(-a2*a0**(-1)*a1 + a3)**(-1)*(-a2*a0**(-1)*b0 + b1) + b0)
-        + a3*(-a2*a0**(-1)*a1 + a3)**(-1)*(-a2*a0**(-1)*b0 + b1)
-    )
-
-
 def test_QRsolve():
     A = Matrix([[2, 3, 5],
                 [3, 6, 2],
diff --git a/sympy/physics/quantum/tests/test_dagger.py b/sympy/physics/quantum/tests/test_dagger.py
index 1357c9320a..8aff96506d 100644
--- a/sympy/physics/quantum/tests/test_dagger.py
+++ b/sympy/physics/quantum/tests/test_dagger.py
@@ -20,7 +20,7 @@ def test_scalars():
     assert Dagger(i) == i
 
     p = symbols('p')
-    assert isinstance(Dagger(p), conjugate)
+    assert isinstance(Dagger(p), adjoint)
 
     i = Integer(3)
     assert Dagger(i) == i
@@ -89,7 +89,7 @@ def test_unknown():
     Objects without adjoint or conjugate/transpose methods
     are sympified and wrapped in dagger.
     """
-    x = symbols("x", commutative=False)
+    x = symbols("x")
     result = Dagger(x)
     assert result.args == (x,) and isinstance(result, adjoint)
 
diff --git a/sympy/physics/quantum/tests/test_printing.py b/sympy/physics/quantum/tests/test_printing.py
index ce4004cee2..9e7c485af5 100644
--- a/sympy/physics/quantum/tests/test_printing.py
+++ b/sympy/physics/quantum/tests/test_printing.py
@@ -220,7 +220,7 @@ def test_constants():
 
 
 def test_dagger():
-    x = symbols('x', commutative=False)
+    x = symbols('x')
     expr = Dagger(x)
     assert str(expr) == 'Dagger(x)'
     ascii_str = \
@@ -236,7 +236,7 @@ def test_dagger():
     assert pretty(expr) == ascii_str
     assert upretty(expr) == ucode_str
     assert latex(expr) == r'x^{\dagger}'
-    sT(expr, "Dagger(Symbol('x', commutative=False))")
+    sT(expr, "Dagger(Symbol('x'))")
 
 
 @XFAIL
diff --git a/sympy/printing/tests/test_mathml.py b/sympy/printing/tests/test_mathml.py
index 3d0016f471..88caad269c 100644
--- a/sympy/printing/tests/test_mathml.py
+++ b/sympy/printing/tests/test_mathml.py
@@ -1351,7 +1351,6 @@ def test_print_Contains():
 
 
 def test_print_Dagger():
-    x = symbols('x', commutative=False)
     assert mpp.doprint(Dagger(x)) == '<msup><mi>x</mi>&#x2020;</msup>'
 
 
diff --git a/sympy/printing/tests/test_pycode.py b/sympy/printing/tests/test_pycode.py
index 932b703c28..84ac7c1c87 100644
--- a/sympy/printing/tests/test_pycode.py
+++ b/sympy/printing/tests/test_pycode.py
@@ -1,4 +1,3 @@
-from sympy import Not
 from sympy.codegen import Assignment
 from sympy.codegen.ast import none
 from sympy.codegen.cfunctions import expm1, log1p
@@ -41,7 +40,6 @@ def test_PythonCodePrinter():
     assert prntr.doprint(And(x, y)) == 'x and y'
     assert prntr.doprint(Or(x, y)) == 'x or y'
     assert prntr.doprint(1/(x+y)) == '1/(x + y)'
-    assert prntr.doprint(Not(x)) == 'not x'
     assert not prntr.module_imports
 
     assert prntr.doprint(pi) == 'math.pi'
diff --git a/sympy/simplify/tests/test_powsimp.py b/sympy/simplify/tests/test_powsimp.py
index fdae6bfc1b..61bdc93d05 100644
--- a/sympy/simplify/tests/test_powsimp.py
+++ b/sympy/simplify/tests/test_powsimp.py
@@ -122,6 +122,8 @@ def test_powsimp():
     assert powsimp(x**2*y**3*(x*y**2)**Rational(3, 2)
         ) == x*y*(x*y**2)**Rational(5, 2)
 
+    #issue 27380
+    assert powsimp(1.0**(x+1)/1.0**x) == 1.0
 
 def test_powsimp_negated_base():
     assert powsimp((-x + y)/sqrt(x - y)) == -sqrt(x - y)
diff --git a/sympy/simplify/tests/test_simplify.py b/sympy/simplify/tests/test_simplify.py
index f4392b6693..926b0380cc 100644
--- a/sympy/simplify/tests/test_simplify.py
+++ b/sympy/simplify/tests/test_simplify.py
@@ -232,6 +232,8 @@ def test_issue_5652():
     n = symbols('n', commutative=False)
     assert simplify(n + n**(-n)) == n + n**(-n)
 
+def test_issue_27380():
+    assert simplify(1.0**(x+1)/1.0**x) == 1.0
 
 def test_simplify_fail1():
     x = Symbol('x')
diff --git a/sympy/solvers/tests/test_polysys.py b/sympy/solvers/tests/test_polysys.py
index b90803873b..04c1884231 100644
--- a/sympy/solvers/tests/test_polysys.py
+++ b/sympy/solvers/tests/test_polysys.py
@@ -17,7 +17,7 @@
                                    solve_biquadratic, SolveFailed,
                                    solve_generic, factor_system_bool,
                                    factor_system_cond, factor_system_poly,
-                                   factor_system, _factor_sets, _factor_sets_slow)
+                                   factor_system)
 from sympy.polys.polytools import parallel_poly_from_expr
 from sympy.testing.pytest import raises
 from sympy.core.relational import Eq
@@ -187,271 +187,231 @@ def test_solve_issue_3686():
 
 def test_factor_system():
 
-    assert factor_system([x**2 + 2*x + 1]) ==  [[x + 1]]
-    assert factor_system([x**2 + 2*x + 1, y**2 + 2*y + 1]) ==  [[x + 1, y + 1]]
-    assert factor_system([x**2 + 1]) ==  [[x**2 + 1]]
-    assert factor_system([]) == [[]]
-
-    assert factor_system([x**2 + y**2 + 2*x*y, x**2 - 2], extension=sqrt(2)) == [
-        [x + y, x + sqrt(2)],
-        [x + y, x - sqrt(2)],
-    ]
-
-    assert factor_system([x**2 + 1, y**2 + 1], gaussian=True) == [
-        [x + I, y + I],
-        [x + I, y - I],
-        [x - I, y + I],
-        [x - I, y - I],
-    ]
-
-    assert factor_system([x**2 + 1, y**2 + 1], domain=QQ_I) == [
-        [x + I, y + I],
-        [x + I, y - I],
-        [x - I, y + I],
-        [x - I, y - I],
-    ]
-
-    assert factor_system([0]) == [[]]
-    assert factor_system([1]) == []
-    assert factor_system([0 , x]) == [[x]]
-    assert factor_system([1, 0, x]) == []
-
-    assert factor_system([x**4 - 1, y**6 - 1]) == [
-        [x**2 + 1, y**2 + y + 1],
-        [x**2 + 1, y**2 - y + 1],
-        [x**2 + 1, y + 1],
-        [x**2 + 1, y - 1],
-        [x + 1, y**2 + y + 1],
-        [x + 1, y**2 - y + 1],
-        [x - 1, y**2 + y + 1],
-        [x - 1, y**2 - y + 1],
-        [x + 1, y + 1],
-        [x + 1, y - 1],
-        [x - 1, y + 1],
-        [x - 1, y - 1],
-    ]
-
-    assert factor_system([(x - 1)*(y - 2), (y - 2)*(z - 3)]) == [
-        [x - 1, z - 3],
-        [y - 2]
-    ]
-
-    assert factor_system([sin(x)**2 + cos(x)**2 - 1, x]) == [
-        [x, sin(x)**2 + cos(x)**2 - 1],
-    ]
-
-    assert factor_system([sin(x)**2 + cos(x)**2 - 1]) == [
-        [sin(x)**2 + cos(x)**2 - 1]
-    ]
-
-    assert factor_system([sin(x)**2 + cos(x)**2]) == [
-        [sin(x)**2 + cos(x)**2]
-    ]
-
-    assert factor_system([a*x, y, a]) == [[y, a]]
-
-    assert factor_system([a*x, y, a], [x, y]) == []
-
-    assert factor_system([a ** 2 * x, y], [x, y]) == [[x, y]]
-
-    assert factor_system([a*x*(x - 1), b*y, c], [x, y]) == []
-
-    assert factor_system([a*x*(x - 1), b*y, c], [x, y, c]) == [
-        [x - 1, y, c],
-        [x, y, c],
-    ]
-
-    assert factor_system([a*x*(x - 1), b*y, c]) == [
-        [x - 1, y, c],
-        [x, y, c],
-        [x - 1, b, c],
-        [x, b, c],
-        [y, a, c],
-        [a, b, c],
-    ]
-
-    assert factor_system([x**2 - 2], [y]) == []
-
-    assert factor_system([x**2 - 2], [x]) == [[x**2 - 2]]
-
-    assert factor_system([cos(x)**2 - sin(x)**2, cos(x)**2 + sin(x)**2 - 1]) == [
-        [sin(x)**2 + cos(x)**2 - 1, sin(x) + cos(x)],
-        [sin(x)**2 + cos(x)**2 - 1, -sin(x) + cos(x)],
-    ]
-
-    assert factor_system([(cos(x) + sin(x))**2 - 1, cos(x)**2 - sin(x)**2 - cos(2*x)]) == [
-        [sin(x)**2 - cos(x)**2 + cos(2*x), sin(x) + cos(x) + 1],
-        [sin(x)**2 - cos(x)**2 + cos(2*x), sin(x) + cos(x) - 1],
-    ]
-
-    assert factor_system([(cos(x) + sin(x))*exp(y) - 1, (cos(x) - sin(x))*exp(y) - 1]) == [
-        [exp(y)*sin(x) + exp(y)*cos(x) - 1, -exp(y)*sin(x) + exp(y)*cos(x) - 1]
-    ]
-
-
-def test_factor_system_poly():
-
-    px = lambda e: Poly(e, x)
-    pxab = lambda e: Poly(e, x, domain=ZZ[a, b])
-    pxI = lambda e: Poly(e, x, domain=QQ_I)
-    pxyz = lambda e: Poly(e, (x, y, z))
-
-    assert factor_system_poly([px(x**2 - 1), px(x**2 - 4)]) == [
-        [px(x + 2), px(x + 1)],
-        [px(x + 2), px(x - 1)],
-        [px(x + 1), px(x - 2)],
-        [px(x - 1), px(x - 2)],
-    ]
-
-    assert factor_system_poly([px(x**2 - 1)]) == [[px(x + 1)], [px(x - 1)]]
-
-    assert factor_system_poly([pxyz(x**2*y - y), pxyz(x**2*z - z)]) == [
-        [pxyz(x + 1)],
-        [pxyz(x - 1)],
-        [pxyz(y), pxyz(z)],
-    ]
+    assert factor_system([x**2 + 2*x + 1]) ==  ([[x + 1]], True)
+    assert factor_system([x**2 + 2*x + 1, y**2 + 2*y + 1]) ==  ([[x + 1, y + 1]], True)
+    assert factor_system([x**2 + 1]) ==  ([[x**2 + 1]], True)
+    assert factor_system([]) == ([[]], True)
+
+    result = factor_system([x**2 + y**2 + 2*x*y, x**2 - 2], extension=sqrt(2))
+    expected = ([[x + y, x - sqrt(2)], [x + y, x + sqrt(2)]], True)
+    assert result == expected
+
+    result = factor_system([x ** 2 + 1, y ** 2 + 1], gaussian=True)
+    expcted = ([[x - I, y - I], [x - I, y + I], [x + I, y - I], [x + I, y + I]], True)
+    assert result == expcted
+
+    result = factor_system([x ** 2 + 1, y ** 2 + 1], domain=QQ_I)
+    expected = ([[x - I, y - I], [x - I, y + I], [x + I, y - I], [x + I, y + I]], True)
+    assert result == expected
+
+    assert factor_system([0]) == ([[]], True)
+    assert factor_system([1]) == ([], True)
+    assert factor_system([0 , x]) == ([[x]], True)
+    assert factor_system([1, 0, x]) == ([], True)
+
+    result = factor_system([x**4 - 1, y**6 - 1])
+    expected = (
+        [
+            [x - 1, y - 1],
+            [x - 1, y**2 - y + 1],
+            [x - 1, y + 1],
+            [x - 1, y**2 + y + 1],
+            [x + 1, y - 1],
+            [x + 1, y**2 - y + 1],
+            [x + 1, y + 1],
+            [x + 1, y**2 + y + 1],
+            [x**2 + 1, y - 1],
+            [x**2 + 1, y**2 - y + 1],
+            [x**2 + 1, y + 1],
+            [x**2 + 1, y**2 + y + 1]
+        ],
+        True
+    )
+    assert result == expected
+
+    result = factor_system([(x - 1)*(y - 2), (y  - 2)*(z - 3)])
+    expected = ([[x - 1, z - 3], [y - 2]], True)
+    assert result == expected
+
+    result = factor_system([sin(x) ** 2 + cos(x) ** 2 - 1, x])
+    expected = ([[sin(x)**2 + cos(x)**2 - 1, x]], True)
+    assert result == expected
+
+    result = factor_system([sin(x)**2 + cos(x)**2 - 1])
+    expected = ([[sin(x)**2 + cos(x)**2 - 1]], True)
+    assert result == expected
+
+    result = factor_system([sin(x) ** 2 + cos(x) ** 2])
+    expected = ([[sin(x)**2 + cos(x)**2]], True)
+    assert result == expected
+
+    result = factor_system([a*x*(x-1), b*y, c], [x, y])
+    expected = ([[x, y], [x - 1, y]], Eq(c, 0))
+    assert result == expected
+
+    result = factor_system([x**2 - 2], [y])
+    expected = ([[]], Eq(x**2 - 2, 0))
+    assert result == expected
+
+    # Generators not symbols
+    expr1 = cos(x) ** 2 - sin(x) ** 2
+    expr2 = cos(x) ** 2 + sin(x) ** 2 - 1
+    result = factor_system([expr1, expr2])
+    expected = (
+        [
+            [sin(x) + cos(x), sin(x)**2 + cos(x)**2 - 1],
+            [-sin(x) + cos(x), sin(x)**2 + cos(x)**2 - 1]],
+        True
+    )
+    assert result == expected
+
+    expr3 = (cos(x) + sin(x)) ** 2 - 1
+    expr4 = cos(x) ** 2 - sin(x) ** 2 - cos(2 * x)
+    result = factor_system([expr3, expr4])
+    expected = (
+        [
+            [sin(x) + cos(x) + 1, sin(x)**2 - cos(x)**2 + cos(2*x)],
+            [sin(x) + cos(x) - 1, sin(x)**2 - cos(x)**2 + cos(2*x)]],
+        True
+    )
+    assert result == expected
 
-    assert factor_system_poly([px(x**2*(x - 1)**2), px(x*(x - 1))]) == [
-        [px(x)],
-        [px(x - 1)],
-    ]
+    expr5 = (cos(x) + sin(x)) * exp(y) - 1
+    expr6 = (cos(x) - sin(x)) * exp(y) - 1
+    result = factor_system([expr5, expr6])
+    expected = ([[exp(y)*sin(x) + exp(y)*cos(x) - 1, -exp(y)*sin(x) + exp(y)*cos(x) - 1]], True)
+    assert result == expected
 
-    assert factor_system_poly([pxyz(x**2 + y*x), pxyz(x**2 + z*x)]) == [
-        [pxyz(x + y), pxyz(x + z)],
-        [pxyz(x)],
-    ]
 
-    assert factor_system_poly([pxab((a - 1)*(x - 2)), pxab((b - 3)*(x - 2))]) == [
-        [pxab(x - 2)],
-        [pxab(a - 1), pxab(b - 3)],
-    ]
+def test_factor_system_poly():
 
-    assert factor_system_poly([pxI(x**2 + 1)]) == [[pxI(x + I)], [pxI(x - I)]]
+    p1 = Poly(x ** 2 - 1, x)
+    p2 = Poly(x ** 2 - 4, x)
+    result = factor_system_poly([p1, p2])
+    expected = ([[(Poly(x + 1, x, domain='ZZ'), ()), (Poly(x - 2, x, domain='ZZ'), ())],
+                 [(Poly(x + 1, x, domain='ZZ'), ()), (Poly(x + 2, x, domain='ZZ'), ())],
+                 [(Poly(x - 1, x, domain='ZZ'), ()), (Poly(x - 2, x, domain='ZZ'), ())],
+                 [(Poly(x - 1, x, domain='ZZ'), ()), (Poly(x + 2, x, domain='ZZ'), ())]], [])
+    assert result == expected
+
+    p = Poly(x ** 2 - 1, x)
+    result = factor_system_poly([p])
+    expected = ([[(Poly(x - 1, x, domain='ZZ'), ())], [(Poly(x + 1, x, domain='ZZ'), ())]], [])
+    assert result == expected
+
+    p1 = Poly(x ** 2 * y - y, x, y, z)
+    p2 = Poly(x ** 2 * z - z, x, y, z)
+    result = factor_system_poly([p1, p2])
+    expected = ([[(Poly(x + 1, x, y, z, domain='ZZ'), ())],
+                 [(Poly(x - 1, x, y, z, domain='ZZ'), ())],
+                 [(Poly(y, x, y, z, domain='ZZ'), ()), (Poly(z, x, y, z, domain='ZZ'), ())]], [])
+    assert result == expected
+
+    p1 = Poly(x ** 2 * (x - 1) ** 2, x)
+    p2 = Poly(x * (x - 1), x)
+    result = factor_system_poly([p1, p2])
+    expected = ([[(Poly(x, x, domain='ZZ'), ())], [(Poly(x - 1, x, domain='ZZ'), ())]], [])
+    assert result == expected
+
+    p1 = Poly(x ** 2 + y * x, x, y, z)
+    p2 = Poly(x ** 2 + z * x, x, y, z)
+    result = factor_system_poly([p1, p2])
+    expected = ([[(Poly(x + y, x, y, z, domain='ZZ'), ()), (Poly(x + z, x, y, z, domain='ZZ'), ())],
+                 [(Poly(x, x, y, z, domain='ZZ'), ())]], [])
+    assert result == expected
+
+    p1 = Poly((a - 1) * (x - 2), x, domain=ZZ[a, b])
+    p2 = Poly((b - 3) * (x - 2), x, domain=ZZ[a, b])
+    result = factor_system_poly([p1, p2])
+    expected = (
+        [[
+            (
+                Poly(x - 2, x, domain='ZZ[a,b]'),
+                (Poly(a - 1, x, domain='ZZ[a,b]'), Poly(b - 3, x, domain='ZZ[a,b]'))
+            )
+        ]],
+        []
+    )
+    assert result == expected
 
-    assert factor_system_poly([]) == [[]]
+    result = factor_system_poly([Poly(x ** 2 + 1, x, domain=QQ_I)])
+    expected = ([[(Poly(x - I, x, domain='QQ_I'), ())], [(Poly(x + I, x, domain='QQ_I'), ())]], [])
+    assert result == expected
 
-    assert factor_system_poly([px(1)]) == []
-    assert factor_system_poly([px(0), px(x)]) == [[px(x)]]
+    assert factor_system_poly([]) == ([[]], [])
+    assert factor_system_poly([Poly(1, x), Poly(x, x)]) == ([], [])
+    assert factor_system_poly([Poly(0, x), Poly(x, x)]) == ([[(Poly(x, x, domain='ZZ'), ())]], [])
 
 
 def test_factor_system_cond():
 
-    assert factor_system_cond([x ** 2 - 1, x ** 2 - 4]) == [
-        [x + 2, x + 1],
-        [x + 2, x - 1],
-        [x + 1, x - 2],
-        [x - 1, x - 2],
-    ]
-
-    assert factor_system_cond([1]) == []
-    assert factor_system_cond([0]) == [[]]
-    assert factor_system_cond([1, x]) == []
-    assert factor_system_cond([0, x]) == [[x]]
-    assert factor_system_cond([]) == [[]]
-
-    assert factor_system_cond([x**2 + y*x]) == [[x + y], [x]]
-
-    assert factor_system_cond([(a - 1)*(x - 2), (b - 3)*(x - 2)], [x]) == [
-        [x - 2],
-        [a - 1, b - 3],
-    ]
-
-    assert factor_system_cond([a * (x - 1), b], [x]) == [[x - 1, b], [a, b]]
+    result = factor_system_cond([x ** 2 - 1, x ** 2 - 4])
+    expected = (
+        [
+            [(x + 1, []), (x - 2, [])],
+            [(x + 1, []), (x + 2, [])],
+            [(x - 1, []), (x - 2, [])],
+            [(x - 1, []), (x + 2, [])]
+        ],
+        []
+    )
+    assert result == expected
 
-    assert factor_system_cond([a*x*(x-1), b*y, c], [x, y]) == [
-        [x - 1, y, c],
-        [x, y, c],
-        [x - 1, b, c],
-        [x, b, c],
-        [y, a, c],
-        [a, b, c],
-    ]
+    assert factor_system_cond([1]) == ([], [])
+    assert factor_system_cond([0]) == ([[]], [])
+    assert factor_system_cond([1, x]) == ([], [])
+    assert factor_system_cond([0, x]) == ([[(x, [])]], [])
+    assert factor_system_cond([]) == ([[]], [])
 
-    assert factor_system_cond([x*(x-1), y], [x, y]) == [[x - 1, y], [x, y]]
+    result = factor_system_cond([x ** 2 + y * x])
+    expected = ([[(x, [])], [(x + y, [])]], [])
+    assert result == expected
 
-    assert factor_system_cond([a*x, y, a], [x, y]) == [[y, a]]
+    result = factor_system_cond([(a - 1) * (x - 2), (b - 3) * (x - 2)], [x])
+    expected = ([[(x - 2, [a - 1, b - 3])]], [])
+    assert result == expected
 
-    assert factor_system_cond([a*x, b*x], [x, y]) == [[x], [a, b]]
+    result = factor_system_cond([a * (x - 1), b], [x])
+    expected = ([[(x - 1, [a])]], [b])
+    assert result == expected
 
-    assert factor_system_cond([a*b*x, y], [x, y]) == [[x, y], [y, a*b]]
+    result = factor_system_cond([a*x*(x-1), b*y, c], [x, y])
+    expected = ([[(x, [a]), (y, [b])], [(x - 1, [a]), (y, [b])]], [c])
+    assert result == expected
 
-    assert factor_system_cond([a*b*x, y]) == [[x, y], [y, a], [y, b]]
+    result = factor_system_cond([x*(x-1), y], [x, y])
+    expected = ([[(x, []), (y, [])], [(x - 1, []), (y, [])]], [])
+    assert result == expected
 
-    assert factor_system_cond([a**2*x, y], [x, y]) == [[x, y], [y, a]]
 
 def test_factor_system_bool():
 
-    eqs = [a*(x - 1)*(y - 1), b*(x - 2)*(y - 1)*(y - 2)]
-    assert factor_system_bool(eqs, [x, y]) == (
-        Eq(y - 1, 0)
-        | (Eq(a, 0) & Eq(b, 0))
-        | (Eq(a, 0) & Eq(x - 2, 0))
-        | (Eq(a, 0) & Eq(y - 2, 0))
-        | (Eq(b, 0) & Eq(x - 1, 0))
-        | (Eq(x - 2, 0) & Eq(x - 1, 0))
-        | (Eq(x - 1, 0) & Eq(y - 2, 0))
-    )
+    eqs = [a * (x - 1) * (y - 1), b * (x - 2) * (y - 1) * (y - 2)]
+    expected = (Eq(y - 1, 0) |
+                (Eq(a, 0) & Eq(b, 0)) |
+                ((Eq(a, 0) | Eq(x - 1, 0)) & (Eq(b, 0) | Eq(x - 2, 0))) |
+                ((Eq(a, 0) | Eq(x - 1, 0)) & (Eq(b, 0) | Eq(y - 2, 0))))
+    assert factor_system_bool(eqs, [x, y]) == expected
 
-    assert factor_system_bool([x - 1], [x]) == Eq(x - 1, 0)
+    eqs = [x - 1]
+    assert factor_system_bool(eqs, [x]) == Eq(x - 1, 0)
 
-    assert factor_system_bool([(x - 1)*(x - 2)], [x]) == Eq(x - 2, 0) | Eq(x - 1, 0)
+    eqs = [(x - 1) * (x - 2)]
+    assert factor_system_bool(eqs, [x]) == Eq(x - 2, 0) | Eq(x - 1, 0)
 
     assert factor_system_bool([], [x]) == True
     assert factor_system_bool([0], [x]) == True
     assert factor_system_bool([1], [x]) == False
     assert factor_system_bool([a], [x]) == Eq(a, 0)
 
-    assert factor_system_bool([a * x, y, a], [x, y]) == Eq(a, 0) & Eq(y, 0)
-
-    assert (factor_system_bool([a*x, b*y*x, a], [x, y]) == (
-        Eq(a, 0) & Eq(b, 0))
-        | (Eq(a, 0) & Eq(x, 0))
-        | (Eq(a, 0) & Eq(y, 0)))
-
-    assert (factor_system_bool([a*x, b*x], [x, y]) == Eq(x, 0) |
-            (Eq(a, 0) & Eq(b, 0)))
-
-    assert (factor_system_bool([a*b*x, y], [x, y]) == (
-        Eq(x, 0) & Eq(y, 0)) |
-        (Eq(y, 0) & Eq(a*b, 0)))
-
-    assert (factor_system_bool([a**2*x, y], [x, y]) == (
-        Eq(a, 0) & Eq(y, 0)) |
-        (Eq(x, 0) & Eq(y, 0)))
-
-    assert factor_system_bool([a*x*y, b*y*z], [x, y, z]) == (
+    eqs = [a * x * y, b * y * z]
+    result = factor_system_bool(eqs, [x, y, z])
+    expected = (
         Eq(y, 0)
         | (Eq(a, 0) & Eq(b, 0))
-        | (Eq(a, 0) & Eq(z, 0))
-        | (Eq(b, 0) & Eq(x, 0))
-        | (Eq(x, 0) & Eq(z, 0))
+        | ((Eq(a, 0) | Eq(x, 0)) & (Eq(b, 0) | Eq(z, 0)))
     )
+    assert result == expected
 
-    assert factor_system_bool([a*(x - 1), b], [x]) == (
-        (Eq(a, 0) & Eq(b, 0))
-        | (Eq(x - 1, 0) & Eq(b, 0))
-    )
-
-
-def test_factor_sets():
-    #
-    from random import randint
-
-    def generate_random_system(n_eqs=3, n_factors=2, max_val=10):
-        return [
-            [randint(0, max_val) for _ in range(randint(1, n_factors))]
-            for _ in range(n_eqs)
-        ]
-
-    test_cases = [
-        [[1, 2], [1, 3]],
-        [[1, 2], [3, 4]],
-        [[1], [1, 2], [2]],
-    ]
-
-    for case in test_cases:
-        assert _factor_sets(case) == _factor_sets_slow(case)
-
-    for _ in range(100):
-        system = generate_random_system()
-        assert _factor_sets(system) == _factor_sets_slow(system)
+    eqs = [a * (x - 1), b]
+    assert factor_system_bool(eqs, [x]) == Eq(b, 0) & (Eq(a, 0) | Eq(x - 1, 0))