diff --git a/Makefile b/Makefile
index 15b3afa..e503456 100644
--- a/Makefile
+++ b/Makefile
@@ -42,12 +42,15 @@ endif
 CFLAGS += $(shell $(PYTHON) tools/python_include_dirs.py)
 
 # Set external projects and their include directories
-DEPS := $(addprefix deps/,eigen spectra parallel-hashmap pybind11)
-CFLAGS += $(addprefix -Ideps/,eigen spectra/include parallel-hashmap pybind11/include)
+DEPS := $(addprefix deps/,eigen spectra parallel-hashmap pybind11 rapidhash)
+CFLAGS += $(addprefix -Ideps/,eigen spectra/include parallel-hashmap pybind11/include rapidhash)
 
 # This C++ compile flag is needed in order for Macs to find system libraries
 ifeq ($(shell uname -s),Darwin)
 CFLAGS += -undefined dynamic_lookup
+PYCI_EXTENSION := dylib
+else
+PYCI_EXTENSION := so
 endif
 
 # Set PyCI version number
@@ -61,6 +64,9 @@ DEFS := -D_PYCI_VERSION='$(PYCI_VERSION)'
 DEFS += -D_GIT_BRANCH='$(shell git rev-parse --abbrev-ref HEAD)'
 DEFS += -D_BUILD_TIME='$(shell date -u +%F\ %T)'
 DEFS += -D_COMPILER_VERSION='$(shell $(CXX) --version | head -n 1)'
+ifeq ($(shell uname -s),Darwin)
+DEFS += -D_USE_RAPIDHASH='1'
+endif
 
 # Set objects
 OBJECTS := $(patsubst %.cpp,%.o,$(wildcard pyci/src/*.cpp))
@@ -70,7 +76,7 @@ OBJECTS := $(patsubst %.cpp,%.o,$(wildcard pyci/src/*.cpp))
 # -------------
 
 .PHONY: all
-all: pyci/_pyci.so.$(PYCI_VERSION) pyci/_pyci.so.$(VERSION_MAJOR) pyci/_pyci.so
+all: pyci/_pyci.$(PYCI_EXTENSION).$(PYCI_VERSION) pyci/_pyci.$(PYCI_EXTENSION).$(VERSION_MAJOR) pyci/_pyci.$(PYCI_EXTENSION)
 
 .PHONY: test
 test:
@@ -78,29 +84,30 @@ test:
 
 .PHONY: clean
 clean:
-	rm -rf pyci/src/*.o pyci/_pyci.so*
+	rm -rf pyci/src/*.o pyci/_pyci.$(PYCI_EXTENSION)*
 
 .PHONY: cleandeps
 cleandeps:
 	rm -rf deps
 
+.PHONY: compile_flags.txt
+compile_flags.txt:
+	echo "$(CFLAGS)" | tr ' ' '\n' > $(@)
+
 
 # Make targets
 # ------------
 
-compile_flags.txt:
-	echo "$(CFLAGS)" | tr ' ' '\n' > $(@)
-
 pyci/src/%.o: pyci/src/%.cpp pyci/include/pyci.h $(DEPS)
 	$(CXX) $(CFLAGS) $(DEFS) -c $(<) -o $(@)
 
-pyci/_pyci.so.$(PYCI_VERSION): $(OBJECTS)
+pyci/_pyci.$(PYCI_EXTENSION).$(PYCI_VERSION): $(OBJECTS)
 	$(CXX) $(CFLAGS) $(DEFS) -shared $(^) -o $(@)
 
-pyci/_pyci.so.$(VERSION_MAJOR): pyci/_pyci.so.$(PYCI_VERSION)
+pyci/_pyci.$(PYCI_EXTENSION).$(VERSION_MAJOR): pyci/_pyci.$(PYCI_EXTENSION).$(PYCI_VERSION)
 	ln -sf $(notdir $(<)) $(@)
 
-pyci/_pyci.so: pyci/_pyci.so.$(PYCI_VERSION)
+pyci/_pyci.$(PYCI_EXTENSION): pyci/_pyci.$(PYCI_EXTENSION).$(PYCI_VERSION)
 	ln -sf $(notdir $(<)) $(@)
 
 deps/eigen:
@@ -114,3 +121,6 @@ deps/parallel-hashmap:
 
 deps/pybind11:
 	[ -d $@ ] || git clone https://github.com/pybind/pybind11.git $@
+
+deps/rapidhash:
+	[ -d $@ ] || git clone https://github.com/nicoshev/rapidhash.git $@
diff --git a/pyci/cepa0.py b/pyci/cepa0.py
new file mode 100644
index 0000000..c1a4049
--- /dev/null
+++ b/pyci/cepa0.py
@@ -0,0 +1,232 @@
+# This file is part of PyCI.
+#
+# PyCI is free software: you can redistribute it and/or modify it under
+# the terms of the GNU General Public License as published by the Free
+# Software Foundation, either version 3 of the License, or (at your
+# option) any later version.
+#
+# PyCI is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+# for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with PyCI. If not, see <http://www.gnu.org/licenses/>.
+
+r"""PyCI CEPA0 module."""
+
+import numpy as np
+
+from scipy.sparse.linalg import lgmres, lsqr
+
+from . import pyci
+
+__all__ = [
+    "CEPA0",
+]
+
+
+class CEPA0:
+    r"""CEPA0 solver class."""
+
+    @property
+    def shape(self) -> tuple[int, int]:
+        r"""
+        Shape of the linear operator.
+
+        Returns
+        -------
+        nrow : int
+            Number of rows.
+        ncol : int
+            Number of columns.
+
+        """
+        return self.op.shape
+
+    @property
+    def dtype(self) -> np.dtype:
+        r"""
+        Data type of the linear operator.
+
+        Returns
+        -------
+        dtype : np.dtype
+            Data type of the linear operator.
+
+        """
+        return self.op.dtype
+
+    def __init__(
+        self,
+        *args: tuple[pyci.wavefunction, pyci.secondquant_op] | pyci.sparse_op,
+        ref: int = 0,
+        damp: float = np.inf,
+    ) -> None:
+        r"""
+        Initialize the CEPA0 solver.
+
+        Parameters
+        ----------
+        args : (pyci.sparse_op,) or (pyci.secondquant_op, pyci.wavefunction)
+            System to solve.
+        ref : int, default=0
+            Index of reference determinant.
+        damp : float, default=np.inf
+            Damping parameter.
+
+        """
+        if len(args) == 1:
+            self.op = args[0]
+            if self.op.symmetric:
+                raise TypeError("Cannot run CEPA0 on symmetric matrix operator")
+            elif self.op.shape[0] < self.op.shape[1]:
+                raise ValueError("Cannot solve underdetermined system (nrow < ncol)")
+        elif len(args) == 2:
+            self.op = pyci.sparse_op(*args, symmetric=False)
+        else:
+            raise ValueError("must pass `ham, wfn` or `op`")
+        self.ref = ref
+        self.h_ref = self.op.get_element(ref, ref)
+        self.damp = -1.0 / damp ** 2
+
+    def __call__(self, x: np.ndarray) -> np.ndarray:
+        r"""
+        Apply the linear operator to a vector ``x``.
+
+        Parameters
+        ----------
+        x : np.ndarray
+            Vector to which linear operator will be applied.
+
+        Returns
+        -------
+        y : np.ndarray
+            Result of applying linear operator to ``x``.
+
+        """
+        return self.matvec(x)
+
+    def matvec(self, x: np.ndarray) -> np.ndarray:
+        r"""
+        Apply the linear operator to a vector ``x``.
+
+        Parameters
+        ----------
+        x : np.ndarray
+            Vector to which linear operator will be applied.
+
+        Returns
+        -------
+        y : np.ndarray
+            Result of applying linear operator to ``x``.
+
+        """
+        op = self.op
+        ref = self.ref
+        h_ref = self.h_ref
+        damp = self.damp
+        y = np.zeros(op.shape[0], dtype=pyci.c_double)
+        for i in range(op.shape[0]):
+            h_ref_damp = np.exp(damp * x[i] ** 2) * h_ref
+            if i == ref:
+                for j in range(self.op.shape[1]):
+                    if j == ref:
+                        y[i] += x[j]
+                    else:
+                        y[i] -= op.get_element(i, j) * x[j]
+            else:
+                for j in range(self.op.shape[1]):
+                    if i == j:
+                        y[i] += (h_ref_damp - op.get_element(i, j)) * x[j]
+                    elif j != ref:
+                        y[i] -= op.get_element(i, j) * x[j]
+        return y
+
+    def rmatvec(self, x: np.ndarray) -> np.ndarray:
+        r"""
+        Apply the transpose of the linear operator to a vector ``x``.
+
+        Parameters
+        ----------
+        x : np.ndarray
+            Vector to which linear operator will be applied.
+
+        Returns
+        -------
+        y : np.ndarray
+            Result of applying linear operator to ``x``.
+
+        """
+        op = self.op
+        h_ref = self.h_ref
+        damp = self.damp
+        y = np.zeros(op.shape[1], dtype=pyci.c_double)
+        for i in range(op.shape[0]):
+            h_ref_damp = np.exp(damp * x[i] ** 2) * h_ref
+            for j in range(self.op.shape[1]):
+                if j == i:
+                    y[j] += (h_ref_damp - op.get_element(i, j)) * x[i]
+                else:
+                    y[j] -= op.get_element(i, j) * x[i]
+        y[self.ref] = x[self.ref]
+        return y
+
+    def rhs(self, x: np.ndarray) -> np.ndarray:
+        r"""
+        Construct the right-hand side vector for the linear system of equations.
+
+        Parameters
+        ----------
+        x : np.ndarray
+            Initial guess.
+
+        Returns
+        -------
+        rhs : np.ndarray
+            Right-hand side vector.
+
+        """
+        op = self.op
+        ref = self.ref
+        y = np.zeros(op.shape[0], dtype=pyci.c_double)
+        for i in range(op.shape[0]):
+            y[i] = op.get_element(i, ref)
+        y -= self.matvec(x)
+        return y
+
+    def solve(
+        self, c0: np.ndarray = None, e0: float = None, maxiter: int = 1000, tol: float = 1.0e-12,
+    ) -> tuple[np.ndarray, np.ndarray]:
+        r"""
+        Solve the CEPA0 problem.
+
+        Parameters
+        ----------
+        c0 : np.ndarray, optional
+            Initial guess for lowest eigenvector.
+        e0 : float, optional
+            Initial guess for energy.
+        maxiter : int, default=1000
+            Maximum number of iterations to perform.
+        tol : float, default=1.0e-12
+            Convergence tolerance.
+
+        Returns
+        -------
+        es : np.ndarray
+            Energy.
+        cs : np.ndarray
+            Coefficient vector.
+
+        """
+        c0 = np.zeros(self.op.shape[1], dtype=pyci.c_double) if c0 is None else c0 / c0[self.ref]
+        c0[self.ref] = self.op.get_element(self.ref, self.ref) if e0 is None else e0 - self.op.ecore
+        if self.op.shape[0] == self.op.shape[1]:
+            cs = lgmres(self, self.rhs(c0), maxiter=maxiter, tol=tol, atol=tol)[0]
+        else:
+            cs = lsqr(self, self.rhs(c0), maxiter=maxiter, atol=tol, btol=tol)[0]
+        cs += c0
+        es = np.full(1, cs[self.ref] + self.op.ecore, dtype=pyci.c_double)
+        cs[self.ref] = 1
+        return es, cs[None, :]
diff --git a/pyci/include/pyci.h b/pyci/include/pyci.h
index 49159d0..8407672 100644
--- a/pyci/include/pyci.h
+++ b/pyci/include/pyci.h
@@ -41,7 +41,11 @@
 
 #include <parallel_hashmap/phmap.h>
 
+#ifdef _USE_RAPIDHASH
+#include <rapidhash.h>
+#else
 #include <SpookyV2.h>
+#endif
 
 #include <sort_with_arg.h>
 
@@ -118,14 +122,24 @@ inline int Ctz(const unsigned long long t) {
 
 /* Hash function. */
 
+#ifdef _USE_RAPIDHASH
+typedef std::uint64_t Hash;
+
+template<typename T, typename U>
+Hash spookyhash(T length, const U *data) {
+    return rapidhash(reinterpret_cast<const void *>(data), length * sizeof(U));
+}
+#else
 typedef std::pair<ulong, ulong> Hash;
 
 template<typename T, typename U>
 Hash spookyhash(T length, const U *data) {
     Hash h(0x23a23cf5033c3c81UL, 0xb3816f6a2c68e530UL);
-    SpookyHash::Hash128(reinterpret_cast<const void *>(data), length * sizeof(U), &h.first, &h.second);
+    SpookyHash::Hash128(reinterpret_cast<const void *>(data), length * sizeof(U), &h.first,
+                        &h.second);
     return h;
 }
+#endif
 
 /* Vector template types. */
 
@@ -162,7 +176,8 @@ template<typename Scalar>
 using Array = pybind11::array_t<Scalar, pybind11::array::c_style | pybind11::array::forcecast>;
 
 template<typename Scalar>
-using ColMajorArray = pybind11::array_t<Scalar, pybind11::array::f_style | pybind11::array::forcecast>;
+using ColMajorArray =
+    pybind11::array_t<Scalar, pybind11::array::f_style | pybind11::array::forcecast>;
 
 /* Forward-declare classes. */
 
@@ -260,14 +275,14 @@ pybind11::tuple py_compute_rdms_fullci(const FullCIWfn &, const Array<double>);
 
 pybind11::tuple py_compute_rdms_genci(const GenCIWfn &, const Array<double>);
 
-pybind11::tuple py_compute_transition_rdms_doci(const DOCIWfn &, const DOCIWfn &, const Array<double>,
-                                           const Array<double>);
+pybind11::tuple py_compute_transition_rdms_doci(const DOCIWfn &, const DOCIWfn &,
+                                                const Array<double>, const Array<double>);
 
-pybind11::tuple py_compute_transition_rdms_fullci(const FullCIWfn &, const FullCIWfn &, const Array<double>,
-                                           const Array<double>);
+pybind11::tuple py_compute_transition_rdms_fullci(const FullCIWfn &, const FullCIWfn &,
+                                                  const Array<double>, const Array<double>);
 
-pybind11::tuple py_compute_transition_rdms_genci(const GenCIWfn &, const GenCIWfn &, const Array<double>,
-                                           const Array<double>);
+pybind11::tuple py_compute_transition_rdms_genci(const GenCIWfn &, const GenCIWfn &,
+                                                 const Array<double>, const Array<double>);
 
 template<class WfnType>
 double py_compute_overlap(const WfnType &, const WfnType &, const Array<double>,
@@ -710,21 +725,19 @@ class Objective {
     std::vector<double> val_paramcons;
 
 public:
+    Objective(const SparseOp &, const Wfn &, const std::size_t = 0UL, const long * = nullptr,
+              const double * = nullptr, const std::size_t = 0UL, const long * = nullptr,
+              const double * = nullptr);
 
-    Objective(const SparseOp &, const Wfn &,
-              const std::size_t = 0UL, const long * = nullptr, const double * = nullptr,
-              const std::size_t = 0UL, const long * = nullptr, const double * = nullptr);
-
-    Objective(const SparseOp &, const Wfn &,
-              const pybind11::object, const pybind11::object,
+    Objective(const SparseOp &, const Wfn &, const pybind11::object, const pybind11::object,
               const pybind11::object, const pybind11::object);
 
     Objective(const Objective &);
 
     Objective(Objective &&) noexcept;
 
-    void init(const std::size_t, const long *, const double *,
-              const std::size_t, const long *, const double *);
+    void init(const std::size_t, const long *, const double *, const std::size_t, const long *,
+              const double *);
 
     void objective(const SparseOp &, const double *, double *);
 
@@ -758,13 +771,12 @@ class AP1roGObjective : public Objective<DOCIWfn> {
     std::vector<std::size_t> part_list;
 
 public:
-    AP1roGObjective(const SparseOp &, const DOCIWfn &,
-                    const std::size_t = 0UL, const long * = nullptr, const double * = nullptr,
-                    const std::size_t = 0UL, const long * = nullptr, const double * = nullptr);
+    AP1roGObjective(const SparseOp &, const DOCIWfn &, const std::size_t = 0UL,
+                    const long * = nullptr, const double * = nullptr, const std::size_t = 0UL,
+                    const long * = nullptr, const double * = nullptr);
 
-    AP1roGObjective(const SparseOp &, const DOCIWfn &,
-                    const pybind11::object, const pybind11::object,
-                    const pybind11::object, const pybind11::object);
+    AP1roGObjective(const SparseOp &, const DOCIWfn &, const pybind11::object,
+                    const pybind11::object, const pybind11::object, const pybind11::object);
 
     AP1roGObjective(const AP1roGObjective &);
 
@@ -790,13 +802,12 @@ class APIGObjective : public Objective<DOCIWfn> {
     std::vector<std::size_t> part_list;
 
 public:
-    APIGObjective(const SparseOp &, const DOCIWfn &,
-                    const std::size_t = 0UL, const long * = nullptr, const double * = nullptr,
-                    const std::size_t = 0UL, const long * = nullptr, const double * = nullptr);
+    APIGObjective(const SparseOp &, const DOCIWfn &, const std::size_t = 0UL,
+                  const long * = nullptr, const double * = nullptr, const std::size_t = 0UL,
+                  const long * = nullptr, const double * = nullptr);
 
-    APIGObjective(const SparseOp &, const DOCIWfn &,
-                    const pybind11::object, const pybind11::object,
-                    const pybind11::object, const pybind11::object);
+    APIGObjective(const SparseOp &, const DOCIWfn &, const pybind11::object, const pybind11::object,
+                  const pybind11::object, const pybind11::object);
 
     APIGObjective(const APIGObjective &);
 
diff --git a/pyci/include/sort_with_arg.h b/pyci/include/sort_with_arg.h
index 30abeae..f666d54 100644
--- a/pyci/include/sort_with_arg.h
+++ b/pyci/include/sort_with_arg.h
@@ -15,7 +15,11 @@
 
 #pragma once
 
+#include <cstddef>
+
 #include <algorithm>
+#include <iterator>
+#include <utility>
 
 /* see <http://stackoverflow.com/a/46370189> */
 
@@ -63,8 +67,13 @@ struct value_reference_t {
 };
 
 template<typename _Data, typename _Order>
-struct value_iterator_t : iterator<random_access_iterator_tag, value_t<_Data, _Order>, ptrdiff_t,
-                                   value_t<_Data, _Order> *, value_reference_t<_Data, _Order>> {
+struct value_iterator_t {
+    using iterator_category = std::random_access_iterator_tag;
+    using difference_type = ptrdiff_t;
+    using value_type = value_t<_Data, _Order>;
+    using pointer = value_t<_Data, _Order> *;
+    using reference = value_reference_t<_Data, _Order>;
+
     _Data *itData;
     _Order *itVal;
     value_iterator_t(_Data *_itData, _Order *_itVal) : itData(_itData), itVal(_itVal) {