diff --git a/src/openfermion/ops/representations/polynomial_tensor.py b/src/openfermion/ops/representations/polynomial_tensor.py
index 9b3f20b47..60e74b6e3 100644
--- a/src/openfermion/ops/representations/polynomial_tensor.py
+++ b/src/openfermion/ops/representations/polynomial_tensor.py
@@ -31,8 +31,7 @@ def general_basis_change(general_tensor, rotation_matrix, key):
     r"""Change the basis of a general interaction tensor.
 
     M'^{p_1p_2...p_n} = R^{p_1}_{a_1} R^{p_2}_{a_2} ...
-                        R^{p_n}_{a_n} M^{a_1a_2...a_n} R^{p_n}_{a_n}^T ...
-                        R^{p_2}_{a_2}^T R_{p_1}_{a_1}^T
+                        R^{p_n}_{a_n} M^{a_1a_2...a_n}
 
     where R is the rotation matrix, M is the general tensor, M' is the
     transformed general tensor, and a_k and p_k are indices. The formula uses
@@ -67,7 +66,7 @@ def general_basis_change(general_tensor, rotation_matrix, key):
 
     # Do the basis change through a single call of numpy.einsum. For example,
     # for the (1, 1, 0, 0) tensor, the call is:
-    #     numpy.einsum('abcd,aA,bB,cC,dD',
+    #     numpy.einsum('abcd,Aa,Bb,Cc,Dd',
     #                  general_tensor,
     #                  rotation_matrix.conj(),
     #                  rotation_matrix.conj(),
@@ -78,7 +77,7 @@ def general_basis_change(general_tensor, rotation_matrix, key):
     subscripts_first = ''.join(chr(ord('a') + i) for i in range(order))
 
     # The 'Aa,Bb,Cc,Dd' part of the subscripts
-    subscripts_rest = ','.join(chr(ord('a') + i) + chr(ord('A') + i) for i in range(order))
+    subscripts_rest = ','.join(chr(ord('A') + i) + chr(ord('a') + i) for i in range(order))
 
     subscripts = subscripts_first + ',' + subscripts_rest
 
diff --git a/src/openfermion/ops/representations/polynomial_tensor_test.py b/src/openfermion/ops/representations/polynomial_tensor_test.py
index 8b166111b..3a8a59ed6 100644
--- a/src/openfermion/ops/representations/polynomial_tensor_test.py
+++ b/src/openfermion/ops/representations/polynomial_tensor_test.py
@@ -474,6 +474,40 @@ def test_rotate_basis_reverse(self):
         polynomial_tensor.rotate_basis(rotation_matrix_reverse)
         self.assertEqual(polynomial_tensor, want_polynomial_tensor)
 
+    def test_rotate_basis_90deg(self):
+        rotation_matrix_90deg = numpy.zeros((self.n_qubits, self.n_qubits))
+        rotation_matrix_90deg[0, 1] = -1
+        rotation_matrix_90deg[1, 0] = 1
+
+        one_body = numpy.zeros((self.n_qubits, self.n_qubits))
+        two_body = numpy.zeros((self.n_qubits, self.n_qubits, self.n_qubits, self.n_qubits))
+        one_body_90deg = numpy.zeros((self.n_qubits, self.n_qubits))
+        two_body_90deg = numpy.zeros((self.n_qubits, self.n_qubits, self.n_qubits, self.n_qubits))
+        i = 0
+        j = 0
+        i_90deg = pow(self.n_qubits, 2) - 1
+        j_90deg = pow(self.n_qubits, 4) - 1
+        for p in range(self.n_qubits):
+            for q in range(self.n_qubits):
+                one_body[p, q] = i
+                i = i + 1
+                one_body_90deg[p, q] = (-1) ** ([p, q].count(0)) * i_90deg
+                i_90deg = i_90deg - 1
+                for r in range(self.n_qubits):
+                    for s in range(self.n_qubits):
+                        two_body[p, q, r, s] = j
+                        j = j + 1
+                        two_body_90deg[p, q, r, s] = (-1) ** ([p, q, r, s].count(0)) * j_90deg
+                        j_90deg = j_90deg - 1
+        polynomial_tensor = PolynomialTensor(
+            {(): self.constant, (1, 0): one_body, (1, 1, 0, 0): two_body}
+        )
+        want_polynomial_tensor = PolynomialTensor(
+            {(): self.constant, (1, 0): one_body_90deg, (1, 1, 0, 0): two_body_90deg}
+        )
+        polynomial_tensor.rotate_basis(rotation_matrix_90deg)
+        self.assertEqual(polynomial_tensor, want_polynomial_tensor)
+
     def test_rotate_basis_quadratic_hamiltonian_real(self):
         self.do_rotate_basis_quadratic_hamiltonian(True)
 
@@ -492,7 +526,7 @@ def do_rotate_basis_quadratic_hamiltonian(self, real):
 
         # Rotate a basis where the Hamiltonian is diagonal
         _, diagonalizing_unitary, _ = quad_ham.diagonalizing_bogoliubov_transform()
-        quad_ham.rotate_basis(diagonalizing_unitary.T)
+        quad_ham.rotate_basis(diagonalizing_unitary)
 
         # Check that the rotated Hamiltonian is diagonal with the correct
         # orbital energies
@@ -515,6 +549,8 @@ def test_rotate_basis_max_order(self):
             self.assertEqual(tensor, want_tensor)
         # I originally wanted to test 25 and 26, but it turns out that
         # numpy.einsum complains "too many subscripts in einsum" before 26.
+        # Currently, the sum of the number of output labels and combined labels
+        # can't not exceed NPY_MAXDIMS(=32).
 
         for order in [27, 28]:
             with self.assertRaises(ValueError):