ronantakizawa/codereview-bench huggingface dataset:
Repository: sympy/sympy
File: sympy/polys/matrices/sdm.py (Python)

Code before:
            raise DMDomainError
        m, n = A.shape
        n2, o = B.shape
        if n != n2:
            raise DMShapeError
        C = sdm_matmul(A, B, A.domain, m, o)
        return A.new(C, (m, o), A.domain)

    def mul(A, b):
        """
        Multiplies each element of A with a scalar b

        Examples
        ========

        >>> from sympy import ZZ
        >>> from sympy.polys.matrices.sdm import SDM
        >>> A = SDM({0:{1: ZZ(2)}, 1:{0:ZZ(1)}}, (2, 2), ZZ)
        >>> A.mul(ZZ(3))
        {0: {1: 6}, 1: {0: 3}}

        """
        return sdm_scalar_mul(A, b, lambda x, y: x * y)

    def rmul(A, b):
        return sdm_scalar_mul(A, b, lambda x, y: y * x)

    def mul_elementwise(A, B):
        if A.domain != B.domain:
            raise DMDomainError
        if A.shape != B.shape:
            raise DMShapeError

        K = A.domain
        zero = K.zero
        if K.is_EXRAW:
            fmul_zero_a = lambda e: e * zero
            fmul_zero_b = lambda e: zero * e
            Csdm = binop_dict(A, B, mul, fmul_zero_a, fmul_zero_b)
        else:
            fzero = lambda e: zero
            Csdm = binop_dict(A, B, mul, fzero, fzero)

        return A.new(Csdm, A.shape, A.domain)

    def add(A, B):
        """

        Adds two :py:class:`~.SDM` matrices

        Examples

Review comment (suggestion):
The call to `A.new` should be here rather than in `sdm_scalar_mul` so that the `sdm_*` functions return consistent types.

Code after:
            raise DMDomainError
        m, n = A.shape
        n2, o = B.shape
        if n != n2:
            raise DMShapeError
        C = sdm_matmul(A, B, A.domain, m, o)
        return A.new(C, (m, o), A.domain)

    def mul(A, b):
        """
        Multiplies each element of A with a scalar b

        Examples
        ========

        >>> from sympy import ZZ
        >>> from sympy.polys.matrices.sdm import SDM
        >>> A = SDM({0:{1: ZZ(2)}, 1:{0:ZZ(1)}}, (2, 2), ZZ)
        >>> A.mul(ZZ(3))
        {0: {1: 6}, 1: {0: 3}}

        """
        Csdm = sdm_scalar_mul(A, b, lambda x, y: x * y, A.domain)
        return A.new(Csdm, A.shape, A.domain)

    def rmul(A, b):
        Csdm = sdm_scalar_mul(A, b, lambda x, y: y * x, A.domain)
        return A.new(Csdm, A.shape, A.domain)

    def mul_elementwise(A, B):
        if A.domain != B.domain:
            raise DMDomainError
        if A.shape != B.shape:
            raise DMShapeError

        K = A.domain
        zero = K.zero
        if K.is_EXRAW:
            fmul_zero_a = lambda e: e * zero
            fmul_zero_b = lambda e: zero * e
            Csdm = binop_dict(A, B, mul, fmul_zero_a, fmul_zero_b)
        else:
            fzero = lambda e: zero
            Csdm = binop_dict(A, B, mul, fzero, fzero)

        return A.new(Csdm, A.shape, A.domain)

    def add(A, B):
        """

        Adds two :py:class:`~.SDM` matrices
Repository: sympy/sympy
File: sympy/matrices/tests/test_eigen.py (Python)

Code before:
                [1, 0, 0],
                [1, 1, 1]])
    assert M.eigenvals() == {2*S.One: 1, -S.One: 1, S.Zero: 1}

    m = Matrix([
        [3,  0,  0, 0, -3],
        [0, -3, -3, 0,  3],
        [0,  3,  0, 3,  0],
        [0,  0,  3, 0,  3],
        [3,  0,  0, 3,  0]])

    # XXX Used dry-run test because arbitrary symbol that appears in
    # CRootOf may not be unique.
    assert m.eigenvals()

    A = Matrix([
        [1, 1, 0, 1, -1, 0, 1],
        [0, 2, 1, 0, 1, -1, 0],
        [1, 0, 1, 1, 0, 1, -1],
        [1, -1, 0, 0, 1, 0, 1],
        [0, 1, 1, -1, 1, 1, 0],
        [1, 0, -1, 1, 0, 1, 1],
        [0, 1, 0, 1, -1, 0, 2]
    ])
    eigs = A.eigenvals()
    assert sum(eigs.values()) == A.charpoly().degree()


def test_eigenvects():
    M = Matrix([[0, 1, 1],
                [1, 0, 0],
                [1, 1, 1]])
    vecs = M.eigenvects()
    for val, mult, vec_list in vecs:
        assert len(vec_list) == 1
        assert M*vec_list[0] == val*vec_list[0]


def test_left_eigenvects():
    M = Matrix([[0, 1, 1],
                [1, 0, 0],
                [1, 1, 1]])
    vecs = M.left_eigenvects()
    for val, mult, vec_list in vecs:
        assert len(vec_list) == 1
        assert vec_list[0]*M == val*vec_list[0]


@slow
def test_bidiagonalize():
    M = Matrix([[1, 0, 0],

Review comment (suggestion):
It is not obvious to me why this equality should hold.

I think it would be better to compare with `A.trace()` rather than `A.charpoly().degree()`. For any square matrix the sum of the eigenvalues (counted with multiplicity) should be equal to the trace. That would seem more obviously correct to me.

Code after:
                [1, 0, 0],
                [1, 1, 1]])
    assert M.eigenvals() == {2*S.One: 1, -S.One: 1, S.Zero: 1}

    m = Matrix([
        [3,  0,  0, 0, -3],
        [0, -3, -3, 0,  3],
        [0,  3,  0, 3,  0],
        [0,  0,  3, 0,  3],
        [3,  0,  0, 3,  0]])

    # XXX Used dry-run test because arbitrary symbol that appears in
    # CRootOf may not be unique.
    assert m.eigenvals()

    A = Matrix([
        [1, 1, 0, 1, -1, 0, 1],
        [0, 2, 1, 0, 1, -1, 0],
        [1, 0, 1, 1, 0, 1, -1],
        [1, -1, 0, 0, 1, 0, 1],
        [0, 1, 1, -1, 1, 1, 0],
        [1, 0, -1, 1, 0, 1, 1],
        [0, 1, 0, 1, -1, 0, 2]
    ])
    eigs = A.eigenvals()
    diff = (A.trace() - sum(eig * mult for eig, mult in eigs.items())).evalf()
    assert abs(diff) < 1e-10


def test_eigenvects():
    M = Matrix([[0, 1, 1],
                [1, 0, 0],
                [1, 1, 1]])
    vecs = M.eigenvects()
    for val, mult, vec_list in vecs:
        assert len(vec_list) == 1
        assert M*vec_list[0] == val*vec_list[0]


def test_left_eigenvects():
    M = Matrix([[0, 1, 1],
                [1, 0, 0],
                [1, 1, 1]])
    vecs = M.left_eigenvects()
    for val, mult, vec_list in vecs:
        assert len(vec_list) == 1
        assert vec_list[0]*M == val*vec_list[0]


@slow
def test_bidiagonalize():
Repository: sympy/sympy
File: sympy/matrices/tests/test_eigen.py (Python)

Code before:
                [1, 0, 0],
                [1, 1, 1]])
    assert M.eigenvals() == {2*S.One: 1, -S.One: 1, S.Zero: 1}

    m = Matrix([
        [3,  0,  0, 0, -3],
        [0, -3, -3, 0,  3],
        [0,  3,  0, 3,  0],
        [0,  0,  3, 0,  3],
        [3,  0,  0, 3,  0]])

    # XXX Used dry-run test because arbitrary symbol that appears in
    # CRootOf may not be unique.
    assert m.eigenvals()

    A = Matrix([
        [1, 1, 0, 1, -1, 0, 1],
        [0, 2, 1, 0, 1, -1, 0],
        [1, 0, 1, 1, 0, 1, -1],
        [1, -1, 0, 0, 1, 0, 1],
        [0, 1, 1, -1, 1, 1, 0],
        [1, 0, -1, 1, 0, 1, 1],
        [0, 1, 0, 1, -1, 0, 2]
    ])
    eigs = A.eigenvals()
    assert A.trace() == nsimplify(sum(eig * mult for eig, mult in eigs.items()))


def test_eigenvects():
    M = Matrix([[0, 1, 1],
                [1, 0, 0],
                [1, 1, 1]])
    vecs = M.eigenvects()
    for val, mult, vec_list in vecs:
        assert len(vec_list) == 1
        assert M*vec_list[0] == val*vec_list[0]


def test_left_eigenvects():
    M = Matrix([[0, 1, 1],
                [1, 0, 0],
                [1, 1, 1]])
    vecs = M.left_eigenvects()
    for val, mult, vec_list in vecs:
        assert len(vec_list) == 1
        assert vec_list[0]*M == val*vec_list[0]


@slow
def test_bidiagonalize():
    M = Matrix([[1, 0, 0],

Review comment (refactor):
Rather than using `nsimplify` it would be better to check something like `abs(ans - expected) < 1e-20`.

Code after:
                [1, 0, 0],
                [1, 1, 1]])
    assert M.eigenvals() == {2*S.One: 1, -S.One: 1, S.Zero: 1}

    m = Matrix([
        [3,  0,  0, 0, -3],
        [0, -3, -3, 0,  3],
        [0,  3,  0, 3,  0],
        [0,  0,  3, 0,  3],
        [3,  0,  0, 3,  0]])

    # XXX Used dry-run test because arbitrary symbol that appears in
    # CRootOf may not be unique.
    assert m.eigenvals()

    A = Matrix([
        [1, 1, 0, 1, -1, 0, 1],
        [0, 2, 1, 0, 1, -1, 0],
        [1, 0, 1, 1, 0, 1, -1],
        [1, -1, 0, 0, 1, 0, 1],
        [0, 1, 1, -1, 1, 1, 0],
        [1, 0, -1, 1, 0, 1, 1],
        [0, 1, 0, 1, -1, 0, 2]
    ])
    eigs = A.eigenvals()
    diff = (A.trace() - sum(eig * mult for eig, mult in eigs.items())).evalf()
    assert abs(diff) < 1e-10


def test_eigenvects():
    M = Matrix([[0, 1, 1],
                [1, 0, 0],
                [1, 1, 1]])
    vecs = M.eigenvects()
    for val, mult, vec_list in vecs:
        assert len(vec_list) == 1
        assert M*vec_list[0] == val*vec_list[0]


def test_left_eigenvects():
    M = Matrix([[0, 1, 1],
                [1, 0, 0],
                [1, 1, 1]])
    vecs = M.left_eigenvects()
    for val, mult, vec_list in vecs:
        assert len(vec_list) == 1
        assert vec_list[0]*M == val*vec_list[0]


@slow
def test_bidiagonalize():
Repository: sympy/sympy
File: sympy/vector/tests/test_coordsysrect.py (Python)

Code before:
from sympy.testing.pytest import raises
from sympy.vector.coordsysrect import CoordSys3D
from sympy.vector.operators import gradient
from sympy.vector.scalar import BaseScalar
from sympy.core.function import expand
from sympy.core.numbers import pi
from sympy.core.symbol import symbols
from sympy.functions.elementary.hyperbolic import (cosh, sinh)
from sympy.functions.elementary.miscellaneous import sqrt
from sympy.functions.elementary.trigonometric import (acos, atan2, cos, sin)
from sympy.matrices.dense import zeros
from sympy.matrices.immutable import ImmutableDenseMatrix as Matrix
from sympy.simplify.simplify import simplify
from sympy.vector.functions import express
from sympy.vector.point import Point
from sympy.vector.vector import Vector
from sympy.vector.orienters import (AxisOrienter, BodyOrienter,
                                    SpaceOrienter, QuaternionOrienter)


x, y, z = symbols('x y z')
a, b, c, q = symbols('a b c q')
q1, q2, q3, q4 = symbols('q1 q2 q3 q4')


def test_func_args():
    A = CoordSys3D('A')
    assert A.x.func(*A.x.args) == A.x
    expr = 3*A.x + 4*A.y
    assert expr.func(*expr.args) == expr

Review comment (suggestion):
```suggestion
from sympy.core.function import expand, Function
```

Code after:
from sympy.testing.pytest import raises
from sympy.vector.coordsysrect import CoordSys3D
from sympy.vector.operators import gradient
from sympy.vector.scalar import BaseScalar
from sympy.core.function import expand, Function
from sympy.core.numbers import pi
from sympy.core.symbol import symbols
from sympy.functions.elementary.hyperbolic import (cosh, sinh)
from sympy.functions.elementary.miscellaneous import sqrt
from sympy.functions.elementary.trigonometric import (acos, atan2, cos, sin)
from sympy.matrices.dense import zeros
from sympy.matrices.immutable import ImmutableDenseMatrix as Matrix
from sympy.simplify.simplify import simplify
from sympy.vector.functions import express
from sympy.vector.point import Point
from sympy.vector.vector import Vector
from sympy.vector.orienters import (AxisOrienter, BodyOrienter,
                                    SpaceOrienter, QuaternionOrienter)


x, y, z = symbols('x y z')
a, b, c, q = symbols('a b c q')
q1, q2, q3, q4 = symbols('q1 q2 q3 q4')


def test_func_args():
    A = CoordSys3D('A')
    assert A.x.func(*A.x.args) == A.x
    expr = 3*A.x + 4*A.y
    assert expr.func(*expr.args) == expr
Repository: sympy/sympy
File: sympy/vector/tests/test_coordsysrect.py (Python)

Code before:
        (x, y, z), (x*sin(y/2)*cos(z), x*sin(y)*sin(z), x*cos(y)))))


def test_rotation_trans_equations():
    a = CoordSys3D('a')
    from sympy.core.symbol import symbols
    q0 = symbols('q0')
    assert a._rotation_trans_equations(a._parent_rotation_matrix, a.base_scalars()) == (a.x, a.y, a.z)
    assert a._rotation_trans_equations(a._inverse_rotation_matrix(), a.base_scalars()) == (a.x, a.y, a.z)
    b = a.orient_new_axis('b', 0, -a.k)
    assert b._rotation_trans_equations(b._parent_rotation_matrix, b.base_scalars()) == (b.x, b.y, b.z)
    assert b._rotation_trans_equations(b._inverse_rotation_matrix(), b.base_scalars()) == (b.x, b.y, b.z)
    c = a.orient_new_axis('c', q0, -a.k)
    assert c._rotation_trans_equations(c._parent_rotation_matrix, c.base_scalars()) == \
           (-sin(q0) * c.y + cos(q0) * c.x, sin(q0) * c.x + cos(q0) * c.y, c.z)
    assert c._rotation_trans_equations(c._inverse_rotation_matrix(), c.base_scalars()) == \
           (sin(q0) * c.y + cos(q0) * c.x, -sin(q0) * c.x + cos(q0) * c.y, c.z)


def test_issue_28559():
    R = CoordSys3D('S', transformation=((x, y, z), (x, y, z)),
        variable_names=('a', 'b', 'c'))

    # Test gradient simplification
    eq = gradient(f(R.a, R.b, R.c))
    result = eq.simplify()
    # Should not be 0
    assert result != 0

    # Test simple derivative simplification
    simple_eq = f(R.a).diff(R.a)
    simple_result = simple_eq.simplify()
    assert simple_result == simple_eq


Review comment (suggestion):
```suggestion
    f = Function('f')

    # Test gradient simplification
```

Code after:
        (x, y, z), (x*sin(y/2)*cos(z), x*sin(y)*sin(z), x*cos(y)))))


def test_rotation_trans_equations():
    a = CoordSys3D('a')
    from sympy.core.symbol import symbols
    q0 = symbols('q0')
    assert a._rotation_trans_equations(a._parent_rotation_matrix, a.base_scalars()) == (a.x, a.y, a.z)
    assert a._rotation_trans_equations(a._inverse_rotation_matrix(), a.base_scalars()) == (a.x, a.y, a.z)
    b = a.orient_new_axis('b', 0, -a.k)
    assert b._rotation_trans_equations(b._parent_rotation_matrix, b.base_scalars()) == (b.x, b.y, b.z)
    assert b._rotation_trans_equations(b._inverse_rotation_matrix(), b.base_scalars()) == (b.x, b.y, b.z)
    c = a.orient_new_axis('c', q0, -a.k)
    assert c._rotation_trans_equations(c._parent_rotation_matrix, c.base_scalars()) == \
           (-sin(q0) * c.y + cos(q0) * c.x, sin(q0) * c.x + cos(q0) * c.y, c.z)
    assert c._rotation_trans_equations(c._inverse_rotation_matrix(), c.base_scalars()) == \
           (sin(q0) * c.y + cos(q0) * c.x, -sin(q0) * c.x + cos(q0) * c.y, c.z)


def test_issue_28559():
    R = CoordSys3D('S', transformation=((x, y, z), (x, y, z)),
        variable_names=('a', 'b', 'c'))

    f = Function('f')

    # Test gradient simplification
    eq = gradient(f(R.a, R.b, R.c))
    result = eq.simplify()
    # Should not be 0
    assert result != 0

    # Test simple derivative simplification
    simple_eq = f(R.a).diff(R.a)
    simple_result = simple_eq.simplify()
    assert simple_result == simple_eq
Repository: sympy/sympy
File: sympy/matrices/expressions/determinant.py (Python)

Code before:

    >>> from sympy import MatrixSymbol, Permanent, ones
    >>> A = MatrixSymbol('A', 3, 3)
    >>> Permanent(A)
    Permanent(A)
    >>> Permanent(ones(3, 3)).doit()
    6
    """

    def __new__(cls, mat):
        mat = sympify(mat)
        if not mat.is_Matrix:
            raise TypeError("Input to Permanent, %s, not a matrix" % str(mat))

        return Basic.__new__(cls, mat)

    @property
    def arg(self):
        return self.args[0]

    def doit(self, expand=False, **hints):
        if isinstance(self.arg, MatrixBase):
            if self.arg.shape[0] == 0 or self.arg.shape[1] == 0:
                return S.One
            else:
                return self.arg.per()
        else:
            return self

def per(matexpr):
    """ Matrix Permanent

    Examples
    ========

    >>> from sympy import MatrixSymbol, Matrix, per, ones
    >>> A = MatrixSymbol('A', 3, 3)
    >>> per(A)
    Permanent(A)
    >>> per(ones(5, 5))
    120
    >>> M = Matrix([1, 2, 5])
    >>> per(M)
    8
    """

    return Permanent(matexpr).doit()

from sympy.assumptions.ask import ask, Q
from sympy.assumptions.refine import handlers_dict


Review comment (bug):
The same thing applies as before that this is not the correct place for the check. Checking here means that the Matrix.per method still returns an incorrect result.

The logic that checks this needs to be moved to the place where there is the actual code that computes the permanant.

Code after:
    Represents the permanent of a matrix expression.

    Examples
    ========

    >>> from sympy import MatrixSymbol, Permanent, ones
    >>> A = MatrixSymbol('A', 3, 3)
    >>> Permanent(A)
    Permanent(A)
    >>> Permanent(ones(3, 3)).doit()
    6
    """

    def __new__(cls, mat):
        mat = sympify(mat)
        if not mat.is_Matrix:
            raise TypeError("Input to Permanent, %s, not a matrix" % str(mat))

        return Basic.__new__(cls, mat)

    @property
    def arg(self):
        return self.args[0]

    def doit(self, expand=False, **hints):
        if isinstance(self.arg, MatrixBase):
            return self.arg.per()
        else:
            return self

def per(matexpr):
    """ Matrix Permanent

    Examples
    ========

    >>> from sympy import MatrixSymbol, Matrix, per, ones
    >>> A = MatrixSymbol('A', 3, 3)
    >>> per(A)
    Permanent(A)
    >>> per(ones(5, 5))
    120
    >>> M = Matrix([1, 2, 5])
    >>> per(M)
    8
    """

    return Permanent(matexpr).doit()

from sympy.assumptions.ask import ask, Q
from sympy.assumptions.refine import handlers_dict
Repository: sympy/sympy
File: sympy/ntheory/residue_ntheory.py (Python)

Code before:
    >>> from sympy.ntheory.residue_ntheory import _sqrt_mod_tonelli_shanks
    >>> _sqrt_mod_tonelli_shanks(2, 17) in [6, 11]
    True

    References
    ==========

    .. [1] Carl Pomerance, Richard Crandall, Prime Numbers: A Computational Perspective,
           2nd Edition (2005), page 101, ISBN:978-0387252827

    """
    s = bit_scan1(p - 1)
    t = p >> s
    # find a non-quadratic residue
    if p % 12 == 5:
        # Legendre symbol (3/p) == -1 if p % 12 in [5, 7]
        d = 3
    elif p % 5 in [2, 3]:
        # Legendre symbol (5/p) == -1 if p % 5 in [2, 3]
        d = 5
    else:
        while 1:
            d = randint(6, p - 1)
            if jacobi(d, p) == -1:
                break
    #assert legendre_symbol(d, p) == -1
    A = pow(a, t, p)
    D = pow(d, t, p)
    m = 0
    for i in range(s):
        adm = A*pow(D, m, p) % p
        adm = pow(adm, 2**(s - 1 - i), p)
        if adm % p == p - 1:
            m += 2**i
    #assert A*pow(D, m, p) % p == 1
    x = pow(a, (t + 1)//2, p)*pow(D, m//2, p) % p
    return x


@overload
def sqrt_mod(
    a: SupportsIndex, p: SupportsIndex, all_roots: Literal[False] = False
) -> int | None: ...


@overload
def sqrt_mod(
    a: SupportsIndex, p: SupportsIndex, all_roots: Literal[True]
) -> list[int]: ...



Review comment (bug):
Using `assert False` is problematic because asserts can be disabled with Python's `-O` (optimize) flag, which would cause undefined behavior if the loop completes without finding a suitable value. 

Use `raise RuntimeError("Failed to find a quadratic non-residue")` or similar instead to ensure the error is always raised.
```suggestion
            raise RuntimeError("Failed to find a quadratic non-residue")
```

Code after:
    >>> from sympy.ntheory.residue_ntheory import _sqrt_mod_tonelli_shanks
    >>> _sqrt_mod_tonelli_shanks(2, 17) in [6, 11]
    True

    References
    ==========

    .. [1] Carl Pomerance, Richard Crandall, Prime Numbers: A Computational Perspective,
           2nd Edition (2005), page 101, ISBN:978-0387252827

    """
    s = bit_scan1(p - 1)
    t = p >> s
    # find a non-quadratic residue
    # Noting that we are assuming p=1 (mod 8), we have (d/p) = (p/d)
    # for any odd prime d. Here, (d/p) denotes the Legendre symbol.
    if p % 3 == 2:
        d = 3
    elif p % 5 in [2, 3]:
        d = 5
    else:
        for d in primerange(7, p):
            if jacobi(p, d) == -1:
                break
        else:
            assert False
    #assert legendre_symbol(d, p) == -1
    A = pow(a, t, p)
    D = pow(d, t, p)
    m = 0
    for i in range(s):
        adm = A*pow(D, m, p) % p
        adm = pow(adm, 2**(s - 1 - i), p)
        if adm % p == p - 1:
            m += 2**i
    #assert A*pow(D, m, p) % p == 1
    x = pow(a, (t + 1)//2, p)*pow(D, m//2, p) % p
    return x


@overload
def sqrt_mod(
    a: SupportsIndex, p: SupportsIndex, all_roots: Literal[False] = False
) -> int | None: ...


@overload
def sqrt_mod(
    a: SupportsIndex, p: SupportsIndex, all_roots: Literal[True]
) -> list[int]: ...
Repository: sympy/sympy
File: sympy/logic/boolalg.py (Python)

Code before:
    """
    Converts a term in the expansion of a function from binary to its
    variable form (for ANF).

    Parameters
    ==========

    term : list of 1's and 0's (complementation pattern)
    variables : list of variables

    """
    temp = [variables[n] for n, t in enumerate(term) if t == 1]

    if not temp:
        return true

    return And(*temp)


def _get_odd_parity_terms(n):
    """
    Returns a list of lists, with all possible combinations of n zeros and ones
    with an odd number of ones.
    """
    return [[1 if (mask >> i) & 1 else 0 for i in range(n)]
            for mask in range(1 << n) if bin(mask).count("1") % 2 == 1]


def _get_even_parity_terms(n):
    """
    Returns a list of lists, with all possible combinations of n zeros and ones
    with an even number of ones.
    """
    return [[1 if (mask >> i) & 1 else 0 for i in range(n)]
            for mask in range(1 << n) if bin(mask).count("1") % 2 == 0]


def _simplified_pairs(terms):
    """
    Reduces a set of minterms, if possible, to a simplified set of minterms
    with one less variable in the terms using QM method.
    """
    if not terms:
        return []

    simplified_terms = []
    todo = list(range(len(terms)))

    # Count number of ones as _check_pair can only potentially match if there
    # is at most a difference of a single one
    termdict = defaultdict(list)

Review comment (suggestion):
There should be a function for this that checks the Python version. Then it can be removed when Python 3.9 is not supported any more.

Code after:
    """
    temp = [variables[n] if val == 0 else Not(variables[n])
            for n, val in enumerate(maxterm) if val != 3]
    return Or(*temp)


def _convert_to_varsANF(term, variables):
    """
    Converts a term in the expansion of a function from binary to its
    variable form (for ANF).

    Parameters
    ==========

    term : list of 1's and 0's (complementation pattern)
    variables : list of variables

    """
    temp = [variables[n] for n, t in enumerate(term) if t == 1]

    if not temp:
        return true

    return And(*temp)


def _get_odd_parity_terms(n):
    """
    Returns a list of lists, with all possible combinations of n zeros and ones
    with an odd number of ones.
    """
    return [[1 if (mask >> i) & 1 else 0 for i in range(n)]
            for mask in range(1 << n) if _bit_count(mask) % 2 == 1]


def _get_even_parity_terms(n):
    """
    Returns a list of lists, with all possible combinations of n zeros and ones
    with an even number of ones.
    """
    return [[1 if (mask >> i) & 1 else 0 for i in range(n)]
            for mask in range(1 << n) if _bit_count(mask) % 2 == 0]


def _simplified_pairs(terms):
    """
    Reduces a set of minterms, if possible, to a simplified set of minterms
    with one less variable in the terms using QM method.
    """
    if not terms:
        return []