Repository: sympy/sympy

Issue:
normal_ordered_form works incorrectly with multiple fermions
Normal ordering of fermionic operators does not follow anticommutation rules of different fermions.
For example,
```py
from sympy.physics.quantum.fermion import FermionOp
from sympy.physics.quantum.operatorordering import normal_ordered_form

F = FermionOp("F")
G = FermionOp("G")

normal_ordered_form(G * F, True)
```
should give the result `-FG` or `GF`. However, the code returns `FG`.

The following example 
```py
normal_ordered_form(G * F + F * G, True)
```
should be $0$, but it returns `2FG`.


Discussion:

Patch:
diff --git a/.mailmap b/.mailmap
index cc0769a17b..f1e4d265e5 100644
--- a/.mailmap
+++ b/.mailmap
@@ -667,6 +667,7 @@ Ilya Pchelintsev <ilya.pchelintsev@outlook.com> ILLIA PCHALINTSAU <ilya.pchelint
 Imran Ahmed Manzoor <imran.manzoor31@gmail.com>
 Ishan Joshi <ishanaj98@gmail.com>
 Ishan Pandhare <ishan9096137017@gmail.com>
+Isidora Araya <iarayaday@gmail.com>
 Islam Mansour <is3mansour@gmail.com>
 Islem BOUZENIA <fi_bouzenia@esi.dz> islem-esi <fi_bouzenia@esi.dz>
 Isuru Fernando <isuruf@gmail.com> <isuru.11@cse.mrt.ac.lk>
diff --git a/sympy/physics/quantum/operatorordering.py b/sympy/physics/quantum/operatorordering.py
index 48f4591a1c..d6ba3dd83b 100644
--- a/sympy/physics/quantum/operatorordering.py
+++ b/sympy/physics/quantum/operatorordering.py
@@ -6,7 +6,7 @@
 from sympy.core.mul import Mul
 from sympy.core.numbers import Integer
 from sympy.core.power import Pow
-from sympy.physics.quantum import Operator, Commutator, AntiCommutator
+from sympy.physics.quantum import Commutator, AntiCommutator
 from sympy.physics.quantum.boson import BosonOp
 from sympy.physics.quantum.fermion import FermionOp
 
@@ -35,7 +35,6 @@ def _expand_powers(factors):
 
     return new_factors
 
-
 def _normal_ordered_form_factor(product, independent=False, recursive_limit=10,
                                 _recursive_depth=0):
     """
@@ -51,86 +50,45 @@ def _normal_ordered_form_factor(product, independent=False, recursive_limit=10,
     new_factors = []
     n = 0
     while n < len(factors) - 1:
-
-        if isinstance(factors[n], BosonOp):
-            # boson
-            if not isinstance(factors[n + 1], BosonOp):
-                new_factors.append(factors[n])
-
-            elif factors[n].is_annihilation == factors[n + 1].is_annihilation:
-                if (independent and
-                        str(factors[n].name) > str(factors[n + 1].name)):
-                    new_factors.append(factors[n + 1])
-                    new_factors.append(factors[n])
-                    n += 1
-                else:
-                    new_factors.append(factors[n])
-
-            elif not factors[n].is_annihilation:
-                new_factors.append(factors[n])
-
-            else:
-                if factors[n + 1].is_annihilation:
-                    new_factors.append(factors[n])
-                else:
-                    if factors[n].args[0] != factors[n + 1].args[0]:
-                        if independent:
-                            c = 0
-                        else:
-                            c = Commutator(factors[n], factors[n + 1])
-                        new_factors.append(factors[n + 1] * factors[n] + c)
+        current, next = factors[n], factors[n + 1]
+        if any(not isinstance(f, (FermionOp, BosonOp)) for f in (current, next)):
+            new_factors.append(current)
+            n += 1
+            continue
+
+        key_1 = (current.is_annihilation, str(current.name))
+        key_2 = (next.is_annihilation, str(next.name))
+
+        if key_1 <= key_2:
+            new_factors.append(current)
+            n += 1
+            continue
+
+        n += 2
+        if current.is_annihilation and not next.is_annihilation:
+            if isinstance(current, BosonOp) and isinstance(next, BosonOp):
+                if current.args[0] != next.args[0]:
+                    if independent:
+                        c = 0
                     else:
-                        c = Commutator(factors[n], factors[n + 1])
-                        new_factors.append(
-                            factors[n + 1] * factors[n] + c.doit())
-                    n += 1
-
-        elif isinstance(factors[n], FermionOp):
-            # fermion
-            if not isinstance(factors[n + 1], FermionOp):
-                new_factors.append(factors[n])
-
-            elif factors[n].is_annihilation == factors[n + 1].is_annihilation:
-                if (independent and
-                        str(factors[n].name) > str(factors[n + 1].name)):
-                    new_factors.append(factors[n + 1])
-                    new_factors.append(factors[n])
-                    n += 1
-                else:
-                    new_factors.append(factors[n])
-
-            elif not factors[n].is_annihilation:
-                new_factors.append(factors[n])
-
-            else:
-                if factors[n + 1].is_annihilation:
-                    new_factors.append(factors[n])
+                        c = Commutator(current, next)
+                    new_factors.append(next * current + c)
                 else:
-                    if factors[n].args[0] != factors[n + 1].args[0]:
-                        if independent:
-                            c = 0
-                        else:
-                            c = AntiCommutator(factors[n], factors[n + 1])
-                        new_factors.append(-factors[n + 1] * factors[n] + c)
+                    new_factors.append(next * current + 1)
+            elif isinstance(current, FermionOp) and isinstance(next, FermionOp):
+                if current.args[0] != next.args[0]:
+                    if independent:
+                        c = 0
                     else:
-                        c = AntiCommutator(factors[n], factors[n + 1])
-                        new_factors.append(
-                            -factors[n + 1] * factors[n] + c.doit())
-                    n += 1
-
-        elif isinstance(factors[n], Operator):
-
-            if isinstance(factors[n + 1], (BosonOp, FermionOp)):
-                new_factors.append(factors[n + 1])
-                new_factors.append(factors[n])
-                n += 1
-            else:
-                new_factors.append(factors[n])
-
+                        c = AntiCommutator(current, next)
+                    new_factors.append(-next * current + c)
+                else:
+                    new_factors.append(-next * current + 1)
+        elif (current.is_annihilation == next.is_annihilation and
+            isinstance(current, FermionOp) and isinstance(next, FermionOp)):
+            new_factors.append(-next * current)
         else:
-            new_factors.append(factors[n])
-
-        n += 1
+            new_factors.append(next * current)
 
     if n == len(factors) - 1:
         new_factors.append(factors[-1])
@@ -177,7 +135,10 @@ def normal_ordered_form(expr, independent=False, recursive_limit=10,
 
     expr : expression
         The expression write on normal ordered form.
-
+    independent : bool (default False)
+        Whether to consider operator with different names as operating in
+        different Hilbert spaces. If False, the (anti-)commutation is left
+        explicit.
     recursive_limit : int (default 10)
         The number of allowed recursive applications of the function.
 


Test patch:
diff --git a/sympy/physics/quantum/tests/test_operatorordering.py b/sympy/physics/quantum/tests/test_operatorordering.py
index b4703d846c..f5255d555d 100644
--- a/sympy/physics/quantum/tests/test_operatorordering.py
+++ b/sympy/physics/quantum/tests/test_operatorordering.py
@@ -21,8 +21,10 @@ def test_normal_order():
 
 def test_normal_ordered_form():
     a = BosonOp('a')
+    b = BosonOp('b')
 
     c = FermionOp('c')
+    d = FermionOp('d')
 
     assert normal_ordered_form(Dagger(a) * a) == Dagger(a) * a
     assert normal_ordered_form(a * Dagger(a)) == 1 + Dagger(a) * a
@@ -36,3 +38,13 @@ def test_normal_ordered_form():
     assert normal_ordered_form(c ** 2 * Dagger(c)) == Dagger(c) * c ** 2
     assert normal_ordered_form(c ** 3 * Dagger(c)) == \
         c ** 2 - Dagger(c) * c ** 3
+
+    assert normal_ordered_form(a * Dagger(b), True) == Dagger(b) * a
+    assert normal_ordered_form(Dagger(a) * b, True) == Dagger(a) * b
+    assert normal_ordered_form(b * a, True) == a * b
+    assert normal_ordered_form(Dagger(b) * Dagger(a), True) == Dagger(a) * Dagger(b)
+
+    assert normal_ordered_form(c * Dagger(d), True) == -Dagger(d) * c
+    assert normal_ordered_form(Dagger(c) * d, True) == Dagger(c) * d
+    assert normal_ordered_form(d * c, True) == -c * d
+    assert normal_ordered_form(Dagger(d) * Dagger(c), True) == -Dagger(c) * Dagger(d)