diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 5fd25af39..43ceec2ac 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -2,7 +2,7 @@ ci:
     skip: [mypy]
 repos:
 -   repo: https://github.com/pre-commit/pre-commit-hooks
-    rev: v4.6.0
+    rev: v5.0.0
     hooks:
     -   id: trailing-whitespace
     -   id: end-of-file-fixer
@@ -20,8 +20,8 @@ repos:
     -   id: fix-byte-order-marker
     -   id: fix-encoding-pragma
         args: ['--remove']
--   repo: https://github.com/PyCQA/docformatter
-    rev: v1.7.5
+-   repo: https://github.com/s-weigand/docformatter
+    rev: 1ec30b7
     hooks:
     -   id: docformatter
         args:
@@ -66,7 +66,7 @@ repos:
       args:
         - --in-place
 -   repo: https://github.com/pre-commit/mirrors-mypy
-    rev: v1.11.1
+    rev: v1.11.2
     hooks:
     - id: mypy
       exclude: tests/qis/test_pauli.py
diff --git a/bqskit/ir/gates/__init__.py b/bqskit/ir/gates/__init__.py
index 183a6426c..fb9affdb1 100644
--- a/bqskit/ir/gates/__init__.py
+++ b/bqskit/ir/gates/__init__.py
@@ -72,6 +72,8 @@
     CRZGate
     CUGate
     FSIMGate
+    MPRYGate
+    MPRZGate
     PauliGate
     PauliZGate
     PhasedXZGate
diff --git a/bqskit/ir/gates/parameterized/__init__.py b/bqskit/ir/gates/parameterized/__init__.py
index 546520b60..34b51aee5 100644
--- a/bqskit/ir/gates/parameterized/__init__.py
+++ b/bqskit/ir/gates/parameterized/__init__.py
@@ -11,6 +11,8 @@
 from bqskit.ir.gates.parameterized.crz import CRZGate
 from bqskit.ir.gates.parameterized.cu import CUGate
 from bqskit.ir.gates.parameterized.fsim import FSIMGate
+from bqskit.ir.gates.parameterized.mpry import MPRYGate
+from bqskit.ir.gates.parameterized.mprz import MPRZGate
 from bqskit.ir.gates.parameterized.pauli import PauliGate
 from bqskit.ir.gates.parameterized.pauliz import PauliZGate
 from bqskit.ir.gates.parameterized.phasedxz import PhasedXZGate
@@ -41,6 +43,8 @@
     'CRZGate',
     'CUGate',
     'FSIMGate',
+    'MPRYGate',
+    'MPRZGate',
     'PauliGate',
     'PauliZGate',
     'PhasedXZGate',
diff --git a/bqskit/ir/gates/parameterized/mpry.py b/bqskit/ir/gates/parameterized/mpry.py
new file mode 100644
index 000000000..19f9e8648
--- /dev/null
+++ b/bqskit/ir/gates/parameterized/mpry.py
@@ -0,0 +1,176 @@
+"""This module implements the MPRYGate."""
+from __future__ import annotations
+
+import numpy as np
+import numpy.typing as npt
+
+from bqskit.ir.gates.qubitgate import QubitGate
+from bqskit.qis.unitary.differentiable import DifferentiableUnitary
+from bqskit.qis.unitary.optimizable import LocallyOptimizableUnitary
+from bqskit.qis.unitary.unitary import RealVector
+from bqskit.qis.unitary.unitarymatrix import UnitaryMatrix
+from bqskit.utils.cachedclass import CachedClass
+
+
+def get_indices(
+    index: int,
+    target_qudit: int,
+    num_qudits: int,
+) -> tuple[int, int]:
+    """Get indices for the matrix based on the target qubit."""
+    shift_qubit = num_qudits - target_qudit - 1
+    shift = 2 ** shift_qubit
+    # Split into two parts around target qubit
+    # 100 | 111
+    left = index // shift
+    right = index % shift
+
+    # Now, shift left by one spot to
+    # make room for the target qubit
+    left *= (shift * 2)
+    # Now add 0 * new_ind and 1 * new_ind to get indices
+    return left + right, left + shift + right
+
+
+class MPRYGate(
+    QubitGate,
+    DifferentiableUnitary,
+    CachedClass,
+    LocallyOptimizableUnitary,
+):
+    """
+    A gate representing a multiplexed Y rotation. A multiplexed Y rotation
+    uses n - 1 qubits as select qubits and applies a Y rotation to the target.
+    If the target qubit is the last qubit, then the unitary is block diagonal.
+    Each block is a 2x2 RY matrix with parameter theta.
+
+    Since there are n - 1 select qubits, there are 2^(n-1) parameters (thetas).
+
+    We allow the target qubit to be specified to any qubit, and the other qubits
+    maintain their order. Qubit 0 is the most significant qubit.
+
+    See this paper: https://arxiv.org/pdf/quant-ph/0406176
+    """
+
+    _qasm_name = 'mpry'
+
+    def __init__(
+        self,
+        num_qudits: int,
+        target_qubit: int = -1,
+    ) -> None:
+        self._num_qudits = num_qudits
+        # 1 param for each configuration of the selec qubits
+        self._num_params = 2 ** (num_qudits - 1)
+        # By default, the controlled qubit is the last qubit
+        if target_qubit == -1:
+            target_qubit = num_qudits - 1
+        self.target_qubit = target_qubit
+        super().__init__()
+
+    def get_unitary(self, params: RealVector = []) -> UnitaryMatrix:
+        """Return the unitary for this gate, see :class:`Unitary` for more."""
+        self.check_parameters(params)
+
+        matrix = np.zeros(
+            (
+                2 ** self.num_qudits,
+                2 ** self.num_qudits,
+            ), dtype=np.complex128,
+        )
+        for i, param in enumerate(params):
+            cos = np.cos(param / 2)
+            sin = np.sin(param / 2)
+
+            # Now, get indices based on target qubit.
+            # i corresponds to the configuration of the
+            # select qubits (e.g 5 = 101). Now, the
+            # target qubit is 0,1 for both the row and col
+            # indices. So, if i = 5 and the target_qubit is 2
+            # Then the rows/cols are 1001 and 1101
+            x1, x2 = get_indices(i, self.target_qubit, self.num_qudits)
+
+            matrix[x1, x1] = cos
+            matrix[x2, x2] = cos
+            matrix[x2, x1] = sin
+            matrix[x1, x2] = -1 * sin
+
+        return UnitaryMatrix(matrix)
+
+    def get_grad(self, params: RealVector = []) -> npt.NDArray[np.complex128]:
+        """
+        Return the gradient for this gate.
+
+        See :class:`DifferentiableUnitary` for more info.
+        """
+        self.check_parameters(params)
+
+        grad = np.zeros(
+            (
+                len(params), 2 ** self.num_qudits,
+                2 ** self.num_qudits,
+            ), dtype=np.complex128,
+        )
+
+        # For each parameter, calculate the derivative
+        # with respect to that parameter
+        for i, param in enumerate(params):
+            dcos = -np.sin(param / 2) / 2
+            dsin = -1j * np.cos(param / 2) / 2
+
+            # Again, get indices based on target qubit.
+            x1, x2 = get_indices(i, self.target_qubit, self.num_qudits)
+
+            grad[i, x1, x1] = dcos
+            grad[i, x2, x2] = dcos
+            grad[i, x2, x1] = dsin
+            grad[i, x1, x2] = -1 * dsin
+
+        return grad
+
+    def optimize(self, env_matrix: npt.NDArray[np.complex128]) -> list[float]:
+        """
+        Return the optimal parameters with respect to an environment matrix.
+
+        See :class:`LocallyOptimizableUnitary` for more info.
+        """
+        self.check_env_matrix(env_matrix)
+        thetas: list[float] = [0] * self.num_params
+
+        for i in range(self.num_params):
+            x1, x2 = get_indices(i, self.target_qubit, self.num_qudits)
+            a = np.real(env_matrix[x1, x1] + env_matrix[x2, x2])
+            b = np.real(env_matrix[x2, x1] - env_matrix[x1, x2])
+            theta = 2 * np.arccos(a / np.sqrt(a ** 2 + b ** 2))
+            theta *= -1 if b > 0 else 1
+            thetas[i] = theta
+
+        return thetas
+
+    @staticmethod
+    def get_decomposition(params: RealVector = []) -> tuple[
+        RealVector,
+        RealVector,
+    ]:
+        """
+        Get the corresponding parameters for one level of decomposition of a
+        multiplexed gate.
+
+        This is used in the decomposition of both the MPRY and MPRZ gates.
+        """
+        new_num_params = len(params) // 2
+        left_params = np.zeros(new_num_params)
+        right_params = np.zeros(new_num_params)
+        for i in range(len(left_params)):
+            left_param = (params[i] + params[i + new_num_params]) / 2
+            right_param = (params[i] - params[i + new_num_params]) / 2
+            left_params[i] = left_param
+            right_params[i] = right_param
+
+        return left_params, right_params
+
+    @property
+    def name(self) -> str:
+        """The name of this gate, with the number of qudits appended."""
+        base_name = getattr(self, '_name', self.__class__.__name__)
+        return f'{base_name}_{self.num_qudits}'
diff --git a/bqskit/ir/gates/parameterized/mprz.py b/bqskit/ir/gates/parameterized/mprz.py
new file mode 100644
index 000000000..f35e8747a
--- /dev/null
+++ b/bqskit/ir/gates/parameterized/mprz.py
@@ -0,0 +1,147 @@
+"""This module implements the MPRZGate."""
+from __future__ import annotations
+
+import numpy as np
+import numpy.typing as npt
+
+from bqskit.ir.gates.parameterized.mpry import get_indices
+from bqskit.ir.gates.qubitgate import QubitGate
+from bqskit.qis.unitary.differentiable import DifferentiableUnitary
+from bqskit.qis.unitary.optimizable import LocallyOptimizableUnitary
+from bqskit.qis.unitary.unitary import RealVector
+from bqskit.qis.unitary.unitarymatrix import UnitaryMatrix
+from bqskit.utils.cachedclass import CachedClass
+
+
+class MPRZGate(
+    QubitGate,
+    DifferentiableUnitary,
+    CachedClass,
+    LocallyOptimizableUnitary,
+):
+    """
+    A gate representing a multiplexed Z rotation. A multiplexed Z rotation
+    uses n - 1 qubits as select qubits and applies a Z rotation to the target.
+    If the target qubit is the last qubit, then the unitary is block diagonal.
+    Each block is a 2x2 RZ matrix with parameter theta.
+
+    Since there are n - 1 select qubits, there are 2^(n-1) parameters (thetas).
+
+    We allow the target qubit to be specified to any qubit, and the other qubits
+    maintain their order. Qubit 0 is the most significant qubit.
+
+
+    Why is 0 the MSB? Typically, in the QSD diagram, we see the block drawn
+    with qubit 0 at the top and qubit n-1 at the bottom. Then, the decomposition
+    slowly moves from the bottom to the top.
+
+    See this paper: https://arxiv.org/pdf/quant-ph/0406176
+    """
+
+    _qasm_name = 'mprz'
+
+    def __init__(
+        self,
+        num_qudits: int,
+        target_qubit: int = -1,
+    ) -> None:
+        """
+        Create a new MPRZGate with `num_qudits` qubits and `target_qubit` as the
+        target qubit. We then have 2^(n-1) parameters for this gate.
+
+        For Example:
+        `num_qudits` = 3, `target_qubit` = 1
+
+        Then, the select qubits are 0 and 2 with 0 as the MSB.
+
+        If the input vector is |0x0> then the selection is 00, and
+        RZ(theta_0) is applied to the target qubit.
+
+        If the input vector is |1x0> then the selection is 01, and
+        RZ(theta_1) is applied to the target qubit.
+        """
+        self._num_qudits = num_qudits
+        # 1 param for each configuration of the selec qubits
+        self._num_params = 2 ** (num_qudits - 1)
+        # By default, the controlled qubit is the last qubit
+        if target_qubit == -1:
+            target_qubit = num_qudits - 1
+        self.target_qubit = target_qubit
+        super().__init__()
+
+    def get_unitary(self, params: RealVector = []) -> UnitaryMatrix:
+        """Return the unitary for this gate, see :class:`Unitary` for more."""
+        self.check_parameters(params)
+        matrix = np.zeros(
+            (
+                2 ** self.num_qudits,
+                2 ** self.num_qudits,
+            ), dtype=np.complex128,
+        )
+        for i, param in enumerate(params):
+            pos = np.exp(1j * param / 2)
+            neg = np.exp(-1j * param / 2)
+
+            # Get correct indices based on target qubit
+            # See :class:`mcry` for more info
+            x1, x2 = get_indices(i, self.target_qubit, self.num_qudits)
+
+            matrix[x1, x1] = neg
+            matrix[x2, x2] = pos
+
+        return UnitaryMatrix(matrix)
+
+    def get_grad(self, params: RealVector = []) -> npt.NDArray[np.complex128]:
+        """
+        Return the gradient for this gate.
+
+        See :class:`DifferentiableUnitary` for more info.
+        """
+        self.check_parameters(params)
+
+        grad = np.zeros(
+            (
+                len(params), 2 ** self.num_qudits,
+                2 ** self.num_qudits,
+            ), dtype=np.complex128,
+        )
+
+        # For each parameter, calculate the derivative
+        # with respect to that parameter
+        for i, param in enumerate(params):
+            dpos = 1j * np.exp(1j * param / 2) / 2
+            dneg = -1j * np.exp(-1j * param / 2) / 2
+
+            # Again, get indices based on target qubit.
+            x1, x2 = get_indices(i, self.target_qubit, self.num_qudits)
+
+            grad[i, x1, x1] = dpos
+            grad[i, x2, x2] = dneg
+
+        return grad
+
+    def optimize(self, env_matrix: npt.NDArray[np.complex128]) -> list[float]:
+        """
+        Return the optimal parameters with respect to an environment matrix.
+
+        See :class:`LocallyOptimizableUnitary` for more info.
+        """
+        self.check_env_matrix(env_matrix)
+        thetas: list[float] = [0] * self.num_params
+
+        for i in range(self.num_params):
+            x1, x2 = get_indices(i, self.target_qubit, self.num_qudits)
+            # Optimize each RZ independently from indices
+            # Taken from QFACTOR repo
+            a = np.angle(env_matrix[x1, x1])
+            b = np.angle(env_matrix[x2, x2])
+            # print(thetas)
+            thetas[i] = a - b
+
+        return thetas
+
+    @property
+    def name(self) -> str:
+        """The name of this gate, with the number of qudits appended."""
+        base_name = getattr(self, '_name', self.__class__.__name__)
+        return f'{base_name}_{self.num_qudits}'
diff --git a/tests/ir/gates/parameterized/test_mpry.py b/tests/ir/gates/parameterized/test_mpry.py
new file mode 100644
index 000000000..8d4b7aa58
--- /dev/null
+++ b/tests/ir/gates/parameterized/test_mpry.py
@@ -0,0 +1,72 @@
+"""This module tests the MPRYGate class."""
+from __future__ import annotations
+
+import numpy as np
+import scipy.linalg as la
+from hypothesis import given
+from hypothesis.strategies import floats
+from hypothesis.strategies import integers
+from hypothesis.strategies import lists
+
+from bqskit.ir.gates.constant import PermutationGate
+from bqskit.ir.gates.parameterized import MPRYGate
+from bqskit.ir.gates.parameterized import RYGate
+
+
+@given(
+    lists(
+        elements=floats(
+            allow_nan=False,
+            allow_infinity=False,
+            width=32,
+        ), min_size=2, max_size=16,
+    ),
+)
+def test_get_unitary(thetas: list[float]) -> None:
+    """
+    Test the get_unitary method of the MPRYGate class.
+
+    Use the default target qubit.
+    """
+    # Ensure that len(thetas) is a power of 2
+    # There are 2 ** (n - 1) parameters
+    num_qudits = int(np.log2(len(thetas))) + 1
+    thetas = thetas[:2 ** (num_qudits - 1)]
+    MPRy = MPRYGate(num_qudits=num_qudits)
+    block_unitaries = [RYGate().get_unitary([theta]) for theta in thetas]
+    blocked_unitary = la.block_diag(*block_unitaries)
+    dist = MPRy.get_unitary(thetas).get_distance_from(blocked_unitary)
+    assert dist < 1e-7
+
+
+@given(integers(min_value=0, max_value=4))
+def test_get_unitary_target_select(target_qubit: int) -> None:
+    """Test the get_unitary method of the MPRYGate class when the target qubit
+    is set."""
+    # Create an MPRY gate with 6 qubits and random parameters
+    num_qudits = 6
+    MPRy = MPRYGate(num_qudits=num_qudits, target_qubit=target_qubit)
+    thetas = list(np.random.rand(2 ** (num_qudits - 1)) * 2 * np.pi)
+
+    # Create the block diagonal matrix
+    block_unitaries = [RYGate().get_unitary([theta]) for theta in thetas]
+    blocked_unitary = la.block_diag(*block_unitaries)
+
+    # Apply a permutation transformation
+    # to the block diagonal matrix
+    # Swap the target qubit with the last qubit
+    # perm = np.arange(num_qudits)
+    perm = list(range(num_qudits))
+    for i in range(target_qubit, num_qudits):
+        perm[i] = i + 1
+    perm[-1] = target_qubit
+
+    perm_gate = PermutationGate(num_qudits, perm)
+
+    full_utry = (
+        perm_gate.get_unitary().conj().T
+        @ blocked_unitary @ perm_gate.get_unitary()
+    )
+
+    dist = MPRy.get_unitary(thetas).get_distance_from(full_utry)
+    assert dist < 1e-7
diff --git a/tests/ir/gates/parameterized/test_mprz.py b/tests/ir/gates/parameterized/test_mprz.py
new file mode 100644
index 000000000..0ba506e30
--- /dev/null
+++ b/tests/ir/gates/parameterized/test_mprz.py
@@ -0,0 +1,72 @@
+"""This module tests the MPRZGate class."""
+from __future__ import annotations
+
+import numpy as np
+import scipy.linalg as la
+from hypothesis import given
+from hypothesis.strategies import floats
+from hypothesis.strategies import integers
+from hypothesis.strategies import lists
+
+from bqskit.ir.gates.constant import PermutationGate
+from bqskit.ir.gates.parameterized import MPRZGate
+from bqskit.ir.gates.parameterized import RZGate
+
+
+@given(
+    lists(
+        elements=floats(
+            allow_nan=False,
+            allow_infinity=False,
+            width=32,
+        ), min_size=2, max_size=16,
+    ),
+)
+def test_get_unitary(thetas: list[float]) -> None:
+    """
+    Test the get_unitary method of the MPRZGate class.
+
+    Use the default target qubit.
+    """
+    # Ensure that len(thetas) is a power of 2
+    # There are 2 ** (n - 1) parameters
+    num_qudits = int(np.log2(len(thetas))) + 1
+    thetas = thetas[:2 ** (num_qudits - 1)]
+    MPRy = MPRZGate(num_qudits=num_qudits)
+    block_unitaries = [RZGate().get_unitary([theta]) for theta in thetas]
+    blocked_unitary = la.block_diag(*block_unitaries)
+    dist = MPRy.get_unitary(thetas).get_distance_from(blocked_unitary)
+    assert dist < 1e-7
+
+
+@given(integers(min_value=0, max_value=4))
+def test_get_unitary_target_select(target_qubit: int) -> None:
+    """Test the get_unitary method of the MPRZGate class when the target qubit
+    is set."""
+    # Create an MPRZ gate with 6 qubits and random parameters
+    num_qudits = 6
+    MPRy = MPRZGate(num_qudits=num_qudits, target_qubit=target_qubit)
+    thetas = list(np.random.rand(2 ** (num_qudits - 1)) * 2 * np.pi)
+
+    # Create the block diagonal matrix
+    block_unitaries = [RZGate().get_unitary([theta]) for theta in thetas]
+    blocked_unitary = la.block_diag(*block_unitaries)
+
+    # Apply a permutation transformation
+    # to the block diagonal matrix
+    # Swap the target qubit with the last qubit
+    # perm = np.arange(num_qudits)
+    perm = list(range(num_qudits))
+    for i in range(target_qubit, num_qudits):
+        perm[i] = i + 1
+    perm[-1] = target_qubit
+
+    perm_gate = PermutationGate(num_qudits, perm)
+
+    full_utry = (
+        perm_gate.get_unitary().conj().T
+        @ blocked_unitary @ perm_gate.get_unitary()
+    )
+
+    dist = MPRy.get_unitary(thetas).get_distance_from(full_utry)
+    assert dist < 1e-7