diff --git a/src/openfermion/ops/representations/polynomial_tensor.py b/src/openfermion/ops/representations/polynomial_tensor.py
index 0d06bebb..20942590 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(),
@@ -79,7 +78,7 @@ def general_basis_change(general_tensor, rotation_matrix, key):
 
     # 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))
+        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 c608309f..3852e72f 100644
--- a/src/openfermion/ops/representations/polynomial_tensor_test.py
+++ b/src/openfermion/ops/representations/polynomial_tensor_test.py
@@ -526,6 +526,47 @@ 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)
 
@@ -545,7 +586,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
@@ -567,6 +608,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):