diff --git a/src/openfermion/ops/__init__.py b/src/openfermion/ops/__init__.py
index d604a0845..202a65ba0 100644
--- a/src/openfermion/ops/__init__.py
+++ b/src/openfermion/ops/__init__.py
@@ -13,6 +13,7 @@
 from ._binary_polynomial import BinaryPolynomial
 from ._diagonal_coulomb_hamiltonian import DiagonalCoulombHamiltonian
 from ._indexing import down_index, up_index
+from ._majorana_operator import MajoranaOperator
 from ._polynomial_tensor import PolynomialTensor, general_basis_change
 from ._quadratic_hamiltonian import QuadraticHamiltonian
 from ._symbolic_operator import SymbolicOperator
diff --git a/src/openfermion/ops/_majorana_operator.py b/src/openfermion/ops/_majorana_operator.py
new file mode 100644
index 000000000..a4488fc8e
--- /dev/null
+++ b/src/openfermion/ops/_majorana_operator.py
@@ -0,0 +1,328 @@
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+
+"""The MajoranaOperator data structure."""
+
+from __future__ import division
+from future.utils import viewkeys
+
+import numpy
+
+
+class MajoranaOperator:
+    """A linear combination of products of Majorana operators.
+
+    A system of N fermionic modes can be described using 2N Majorana operators
+    :math:`\gamma_1, \ldots, \gamma_{2N}`
+    as an alternative to using N fermionic annihilation operators. The algebra
+    of Majorana operators amounts to the relation
+
+    .. math::
+        \{\gamma_i, \gamma_j\}
+         = \gamma_i \gamma_j + \gamma_j \gamma_i
+         = 2 \delta_{ij}
+
+    Note that this implies :math:`\gamma_i^2 = 1`.
+
+    The MajoranaOperator class stores a linear combination of products
+    of Majorana operators. Each product is represented as a tuple of
+    integers representing the indices of the operators. As an example,
+    `MajoranaOperator((2, 3, 5), -1.5)` initializes an operator with
+    a single term which represents the operator
+    :math:`-1.5 \gamma_2 \gamma_3 \gamma_5`. MajoranaOperators can be
+    added, subtracted, multiplied, and divided by scalars. They can be
+    compared for approximate numerical equality using `==`.
+
+    Attributes:
+        terms: A dictionary from term, represented by a tuple of integers,
+        to the coefficient of the term in the linear combination.
+    """
+
+    def __init__(self, term=None, coefficient=1.0):
+        """Initialize a MajoranaOperator with a single term.
+
+        Args:
+            term (Tuple[int]): The indices of a Majorana operator term
+                to start off with
+            coefficient (complex): The coefficient of the term
+
+        Returns:
+            MajoranaOperator
+        """
+        self.terms = {}
+        if term is not None:
+            term, parity = _sort_majorana_term(term)
+            self.terms[term] = coefficient * (-1)**parity
+
+    @staticmethod
+    def from_dict(terms):
+        """Initialize a MajoranaOperator from a terms dictionary.
+
+        WARNING: The given dictionary is not validated whatsoever. It's up
+        to you to ensure that it is properly formed.
+
+        Args:
+            terms: A dictionary from Majorana term to coefficient
+        """
+        op = MajoranaOperator()
+        op.terms = terms
+        return op
+
+    def commutes_with(self, other):
+        """Test commutation with another MajoranaOperator"""
+        if not isinstance(other, type(self)):
+            raise TypeError(
+                'Can only test commutation with another MajoranaOperator.')
+
+        if len(self.terms) == 1 and len(other.terms) == 1:
+            return _majorana_terms_commute(list(self.terms.keys())[0],
+                                           list(other.terms.keys())[0])
+        return self*other == other*self
+
+    def __iadd__(self, other):
+        if not isinstance(other, type(self)):
+            return NotImplemented
+
+        for term, coefficient in other.terms.items():
+            if term in self.terms:
+                self.terms[term] += coefficient
+            else:
+                self.terms[term] = coefficient
+
+        return self
+
+    def __add__(self, other):
+        if not isinstance(other, type(self)):
+            return NotImplemented
+
+        terms = {}
+        terms.update(self.terms)
+
+        for term, coefficient in other.terms.items():
+            if term in terms:
+                terms[term] += coefficient
+            else:
+                terms[term] = coefficient
+
+        return MajoranaOperator.from_dict(terms)
+
+    def __isub__(self, other):
+        if not isinstance(other, type(self)):
+            return NotImplemented
+
+        for term, coefficient in other.terms.items():
+            if term in self.terms:
+                self.terms[term] -= coefficient
+            else:
+                self.terms[term] = coefficient
+
+        return self
+
+    def __sub__(self, other):
+        if not isinstance(other, type(self)):
+            return NotImplemented
+
+        terms = {}
+        terms.update(self.terms)
+        for term, coefficient in other.terms.items():
+            if term in terms:
+                terms[term] -= coefficient
+            else:
+                terms[term] = -coefficient
+        return MajoranaOperator.from_dict(terms)
+
+    def __mul__(self, other):
+        if not isinstance(other, (type(self), int, float, complex)):
+            return NotImplemented
+
+        if isinstance(other, (int, float, complex)):
+            terms = {term: coefficient*other
+                     for term, coefficient in self.terms.items()}
+            return MajoranaOperator.from_dict(terms)
+
+        terms = {}
+        for left_term, left_coefficient in self.terms.items():
+            for right_term, right_coefficient in other.terms.items():
+                new_term, parity = _merge_majorana_terms(left_term, right_term)
+                coefficient = left_coefficient*right_coefficient*(-1)**parity
+                if new_term in terms:
+                    terms[new_term] += coefficient
+                else:
+                    terms[new_term] = coefficient
+        return MajoranaOperator.from_dict(terms)
+
+    def __imul__(self, other):
+        if not isinstance(other, (type(self), int, float, complex)):
+            return NotImplemented
+
+        if isinstance(other, (int, float, complex)):
+            for term in self.terms:
+                self.terms[term] *= other
+            return self
+
+        return self * other
+
+    def __rmul__(self, other):
+        if not isinstance(other, (int, float, complex)):
+            return NotImplemented
+        return self * other
+
+    def __truediv__(self, other):
+        if not isinstance(other, (int, float, complex)):
+            return NotImplemented
+
+        terms = {term: coefficient/other
+                 for term, coefficient in self.terms.items()}
+        return MajoranaOperator.from_dict(terms)
+
+    def __itruediv__(self, other):
+        if not isinstance(other, (int, float, complex)):
+            return NotImplemented
+
+        for term in self.terms:
+            self.terms[term] /= other
+        return self
+
+    def __div__(self, divisor):
+        """ For compatibility with Python 2. """
+        return self.__truediv__(divisor)
+
+    def __pow__(self, other):
+        if not isinstance(other, int):
+            return NotImplemented
+
+        if other < 0:
+            raise TypeError('Cannot raise MajoranaOperator to negative power.')
+
+        result = MajoranaOperator(())
+        for _ in range(other):
+            result *= self
+        return result
+
+    def __neg__(self):
+        return -1 * self
+
+    def __eq__(self, other):
+        """Approximate numerical equality."""
+        if not isinstance(other, type(self)):
+            return NotImplemented
+
+        for term in viewkeys(self.terms) | viewkeys(other.terms):
+            if term in self.terms and term in other.terms:
+                if not numpy.isclose(self.terms[term], other.terms[term]):
+                    return False
+            elif term in self.terms:
+                if not numpy.isclose(self.terms[term], 0.0):
+                    return False
+            else:
+                if not numpy.isclose(other.terms[term], 0.0):
+                    return False
+        return True
+
+    def __ne__(self, other):
+        return not self == other
+
+    def __str__(self):
+        if not self.terms:
+            return '0'
+        lines = []
+        for term, coeff in sorted(self.terms.items()):
+            if numpy.isclose(coeff, 0.0):
+                continue
+            lines.append('{} {} +'.format(coeff, term))
+        if not lines:
+            return '0'
+        return '\n'.join(lines)[:-2]
+
+    def __repr__(self):
+        return 'MajoranaOperator.from_dict(terms={!r})'.format(self.terms)
+
+
+def _sort_majorana_term(term):
+    """Sort a Majorana term.
+
+    Args:
+        term (Tuple[int]): The indices of a Majorana operator term
+
+    Returns:
+        Tuple[Tuple[int], int]. The first object returned is a sorted list
+        representing the indices acted upon. The second object is the parity
+        of the term. A parity of 1 indicates that the term should include
+        a minus sign.
+    """
+    if len(term) < 2:
+        return term, 0
+    center = len(term) // 2
+    left_term, left_parity = _sort_majorana_term(term[:center])
+    right_term, right_parity = _sort_majorana_term(term[center:])
+    merged_term, merge_parity = _merge_majorana_terms(left_term, right_term)
+    return merged_term, (left_parity + right_parity + merge_parity) % 2
+
+
+def _merge_majorana_terms(left_term, right_term):
+    """Merge two Majorana terms.
+
+    Args:
+        left_term (Tuple[int]): The left-hand term
+        right_term (Tuple[int]): The right-hand term
+
+    Returns:
+        Tuple[Tuple[int], int]. The first object returned is a sorted list
+        representing the indices acted upon. The second object is the parity
+        of the term. A parity of 1 indicates that the term should include
+        a minus sign.
+    """
+    merged_term = []
+    parity = 0
+    i, j = 0, 0
+    while i < len(left_term) and j < len(right_term):
+        if left_term[i] < right_term[j]:
+            merged_term.append(left_term[i])
+            i += 1
+        elif left_term[i] > right_term[j]:
+            merged_term.append(right_term[j])
+            j += 1
+            parity += len(left_term) - i
+        else:
+            parity += len(left_term) - i - 1
+            i += 1
+            j += 1
+    if i == len(left_term):
+        merged_term.extend(right_term[j:])
+    else:
+        merged_term.extend(left_term[i:])
+    return tuple(merged_term), parity % 2
+
+
+def _majorana_terms_commute(term_a, term_b):
+    """Whether two Majorana terms commute.
+
+    Args:
+        term_a (Tuple[int]): The indices of a Majorana operator term
+        term_b (Tuple[int]): The indices of a Majorana operator term
+
+    Returns:
+        bool. Whether The terms commute.
+    """
+    intersection = 0
+    i, j = 0, 0
+    while i < len(term_a) and j < len(term_b):
+        if term_a[i] < term_b[j]:
+            i += 1
+        elif term_a[i] > term_b[j]:
+            j += 1
+        else:
+            intersection += 1
+            i += 1
+            j += 1
+    parity = (len(term_a)*len(term_b) - intersection) % 2
+    return not parity
diff --git a/src/openfermion/ops/_majorana_operator_test.py b/src/openfermion/ops/_majorana_operator_test.py
new file mode 100644
index 000000000..37d35f83e
--- /dev/null
+++ b/src/openfermion/ops/_majorana_operator_test.py
@@ -0,0 +1,200 @@
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+
+import pytest
+
+from openfermion import MajoranaOperator
+
+
+def test_majorana_operator_init():
+    op = MajoranaOperator((0, 7, 4, 1, 10))
+    assert op.terms == {(0, 1, 4, 7, 10): -1.0}
+
+    op = MajoranaOperator((0, 1, 0, 1, 0))
+    assert op.terms == {(0,): -1.0}
+
+    op = MajoranaOperator((3, 2, 5, 1, 3, 4))
+    assert op.terms == {(1, 2, 4, 5): 1.0}
+
+    op = MajoranaOperator((5, 10, 4, 3, 6, 9, 6))
+    assert op.terms == {(3, 4, 5, 9, 10): -1.0}
+
+
+def test_majorana_operator_commutes_with():
+    a = MajoranaOperator((0, 1, 5))
+    b = MajoranaOperator((1, 2, 7))
+    c = MajoranaOperator((2, 3, 4))
+    d = MajoranaOperator((0, 3, 6))
+
+    assert a.commutes_with(b)
+    assert not a.commutes_with(c)
+    assert a.commutes_with(d)
+    assert b.commutes_with(c)
+    assert not b.commutes_with(d)
+    assert c.commutes_with(d)
+    assert (a+c).commutes_with(b+d)
+
+    e = MajoranaOperator((0, 1, 1, 1, 4, 5))
+    f = MajoranaOperator((0, 1, 1, 4))
+
+    assert e.commutes_with(f)
+
+    with pytest.raises(TypeError):
+        _ = e.commutes_with(0)
+
+
+
+def test_majorana_operator_add_subtract():
+    a = MajoranaOperator((0, 2, 3), -1.25)
+    b = MajoranaOperator((0, 2, 3), -0.5)
+    c = MajoranaOperator((1, 5, 7), 4.75)
+    d = MajoranaOperator((3, 5, 7), 2.25)
+
+    a += c
+    assert a.terms == {(0, 2, 3): -1.25,
+                       (1, 5, 7): 4.75}
+
+    a -= d
+    assert a.terms == {(0, 2, 3): -1.25,
+                       (1, 5, 7): 4.75,
+                       (3, 5, 7): 2.25}
+
+    a += MajoranaOperator((0, 2, 3), 0.5)
+    assert a.terms == {(0, 2, 3): -.75,
+                       (1, 5, 7): 4.75,
+                       (3, 5, 7): 2.25}
+
+    a -= MajoranaOperator((0, 2, 3), 0.25)
+    assert a.terms == {(0, 2, 3): -1.0,
+                       (1, 5, 7): 4.75,
+                       (3, 5, 7): 2.25}
+
+    assert (a + b).terms == {(0, 2, 3): -1.5,
+                             (1, 5, 7): 4.75,
+                             (3, 5, 7): 2.25}
+
+    assert (a - b).terms == {(0, 2, 3): -0.5,
+                             (1, 5, 7): 4.75,
+                             (3, 5, 7): 2.25}
+
+    with pytest.raises(TypeError):
+        _ = a + 0
+
+    with pytest.raises(TypeError):
+        a += 0
+
+    with pytest.raises(TypeError):
+        _ = a - 0
+
+    with pytest.raises(TypeError):
+        a -= 0
+
+
+def test_majorana_operator_multiply():
+    a = MajoranaOperator((0, 1, 5), 1.5) + MajoranaOperator((1, 2, 7), -0.5)
+    b = MajoranaOperator((2, 3, 4), 1.75) - MajoranaOperator((0, 3, 6), 0.25)
+    assert (a * a).terms == {(): -2.5,
+                             (0, 2, 5, 7): -1.5}
+    assert (a * b).terms == {(1, 3, 4, 7): 0.875,
+                             (0, 1, 2, 3, 6, 7): -0.125,
+                             (1, 3, 5, 6): 0.375,
+                             (0, 1, 2, 3, 4, 5): -2.625}
+    assert (2 * a).terms == (a * 2).terms == {(0, 1, 5): 3.0,
+                                              (1, 2, 7): -1.0}
+
+    a *= 2
+    a *= MajoranaOperator(())
+    assert a.terms == {(0, 1, 5): 3.0,
+                       (1, 2, 7): -1.0}
+
+    with pytest.raises(TypeError):
+        _ = a * 'a'
+
+    with pytest.raises(TypeError):
+        _ = 'a' * a
+
+    with pytest.raises(TypeError):
+        a *= 'a'
+
+
+def test_majorana_operator_pow():
+    a = MajoranaOperator((0, 1, 5), 1.5) + MajoranaOperator((1, 2, 7), -0.5)
+    assert (a**2).terms == {(): -2.5,
+                            (0, 2, 5, 7): -1.5}
+
+    with pytest.raises(TypeError):
+        _ = a**-1
+
+    with pytest.raises(TypeError):
+        _ = a**'a'
+
+
+def test_majorana_operator_divide():
+    a = MajoranaOperator((0, 1, 5), 1.5) + MajoranaOperator((1, 2, 7), -0.5)
+    assert (a / 2).terms == {(0, 1, 5): 0.75,
+                             (1, 2, 7): -0.25}
+
+    a /= 2
+    assert a.terms == {(0, 1, 5): 0.75,
+                       (1, 2, 7): -0.25}
+
+    with pytest.raises(TypeError):
+        _ = a / 'a'
+
+    with pytest.raises(TypeError):
+        a /= 'a'
+
+
+def test_majorana_operator_neg():
+    a = MajoranaOperator((0, 1, 5), 1.5) + MajoranaOperator((1, 2, 7), -0.5)
+    assert (-a).terms == {(0, 1, 5): -1.5,
+                          (1, 2, 7): 0.5}
+
+
+def test_majorana_operator_eq():
+    a = MajoranaOperator((0, 1, 5), 1.5) + MajoranaOperator((1, 2, 7), -0.5)
+    b = (MajoranaOperator((0, 1, 5), 1.5) +
+         MajoranaOperator((1, 2, 7), -0.5) +
+         MajoranaOperator((3, 4, 5), 0.0))
+    c = (MajoranaOperator((0, 1, 5), 1.5) +
+         MajoranaOperator((1, 2, 7), -0.5) +
+         MajoranaOperator((3, 4, 5), 0.1))
+    d = MajoranaOperator((0, 1, 5), 1.75) + MajoranaOperator((1, 2, 7), -0.75)
+    e = MajoranaOperator((0, 1, 5), 1.5) - MajoranaOperator((0, 3, 6), 0.25)
+
+    assert a == b
+    assert a != c
+    assert a != d
+    assert a != e
+
+    assert a != 0
+
+
+def test_majorana_operator_str():
+    zero = MajoranaOperator()
+    one = MajoranaOperator(())
+    still_zero = MajoranaOperator((0, 1, 5), 0.0)
+    a = MajoranaOperator((0, 1, 5), 1.5)
+    b = MajoranaOperator((1, 2, 7), -0.5)
+
+    assert str(zero) == '0'
+    assert str(one) == '1.0 ()'
+    assert str(still_zero) == '0'
+    assert str(a) == '1.5 (0, 1, 5)'
+    assert str(b) == '-0.5 (1, 2, 7)'
+    assert str(a+b) == """1.5 (0, 1, 5) +
+-0.5 (1, 2, 7)"""
+
+
+def test_majorana_operator_repr():
+    a = MajoranaOperator((0, 1, 5), 1.5)
+    assert repr(a) == 'MajoranaOperator.from_dict(terms={(0, 1, 5): 1.5})'