diff --git a/src/ruptures/__init__.py b/src/ruptures/__init__.py
index eaf530de..d76e15e4 100644
--- a/src/ruptures/__init__.py
+++ b/src/ruptures/__init__.py
@@ -1,7 +1,16 @@
 """Offline change point detection for Python."""
 
 from .datasets import pw_constant, pw_linear, pw_normal, pw_wavy
-from .detection import Binseg, BottomUp, Dynp, KernelCPD, L1Potts, Pelt, Window
+from .detection import (
+    Binseg,
+    BottomUp,
+    Dynp,
+    EDivisive,
+    KernelCPD,
+    L1Potts,
+    Pelt,
+    Window,
+)
 from .exceptions import NotEnoughPoints
 from .show import display
 
diff --git a/src/ruptures/detection/__init__.py b/src/ruptures/detection/__init__.py
index 649612dc..ecc07c92 100644
--- a/src/ruptures/detection/__init__.py
+++ b/src/ruptures/detection/__init__.py
@@ -3,6 +3,7 @@
 from .binseg import Binseg
 from .bottomup import BottomUp
 from .dynp import Dynp
+from .edivisive import EDivisive
 from .kernelcpd import KernelCPD
 from .l1potts import L1Potts
 from .pelt import Pelt
diff --git a/src/ruptures/detection/edivisive.py b/src/ruptures/detection/edivisive.py
new file mode 100644
index 00000000..73244755
--- /dev/null
+++ b/src/ruptures/detection/edivisive.py
@@ -0,0 +1,244 @@
+r"""E-Divisive change point detection (Matteson & James, 2014)."""
+
+from typing import Optional
+
+import numpy as np
+from numpy.typing import NDArray
+from scipy.spatial.distance import cdist
+from typing_extensions import Self
+
+from ruptures.base import BaseEstimator
+from ruptures.exceptions import BadSegmentationParameters
+
+
+class EDivisive(BaseEstimator):
+    r"""E-Divisive: energy-statistic-based change point detection.
+
+    Detects multiple change points in multivariate time series by
+    recursively splitting segments that maximise the empirical energy
+    divergence.  No distributional assumption is required.
+
+    At each iteration the algorithm selects the split point
+    :math:`\tau^*` and segment :math:`[a, b)` that maximise the
+    scaled energy statistic
+
+    .. math::
+
+        Q(\tau; a, b) = \frac{mn}{m+n}
+        \left(
+          \frac{2}{mn}\sum_{i=a}^{\tau-1}\sum_{j=\tau}^{b-1}
+              \|\mathbf{x}_i - \mathbf{x}_j\|^\alpha
+          - \frac{1}{m^2}\sum_{i,i'=a}^{\tau-1}
+              \|\mathbf{x}_i - \mathbf{x}_{i'}\|^\alpha
+          - \frac{1}{n^2}\sum_{j,j'=\tau}^{b-1}
+              \|\mathbf{x}_j - \mathbf{x}_{j'}\|^\alpha
+        \right),
+
+    where :math:`m = \tau - a` and :math:`n = b - \tau`.
+
+    A permutation test decides whether the best split is significant.
+    Alternatively, pass ``n_bkps`` to fix the number of breakpoints and
+    skip significance testing.
+
+    References:
+        Matteson, D. S. & James, N. A. (2014). A Nonparametric Approach
+        for Multiple Change Point Analysis of Multivariate Data.
+        *Journal of the American Statistical Association*, 109(505), 334-345.
+
+    Example:
+        >>> import numpy as np
+        >>> from ruptures import EDivisive
+        >>> signal = np.concatenate([np.zeros(50), np.ones(50)])
+        >>> EDivisive(n_perms=0).fit_predict(signal, n_bkps=1)
+        [50, 100]
+    """
+
+    def __init__(
+        self,
+        alpha: float = 1.0,
+        min_size: int = 2,
+        n_perms: int = 200,
+        sig_level: float = 0.05,
+    ) -> None:
+        """Initialise EDivisive.
+
+        Args:
+            alpha: exponent applied to pairwise Euclidean distances,
+                must satisfy ``0 < alpha <= 2``. The default ``alpha=1``
+                corresponds to standard Euclidean distances.
+            min_size: minimum number of samples in any segment.
+            n_perms: number of permutations used in the significance
+                test.  Set to ``0`` to disable the test (greedy
+                splitting until ``n_bkps`` is reached).
+            sig_level: significance level for the permutation test.
+                Ignored when ``n_bkps`` is supplied to
+                :meth:`predict`.
+        """
+        if not (0.0 < alpha <= 2.0):
+            raise ValueError("alpha must satisfy 0 < alpha <= 2.")
+        self.alpha = alpha
+        self.min_size = min_size
+        self.n_perms = n_perms
+        self.sig_level = sig_level
+
+        self.signal: Optional[NDArray[np.floating]] = None
+        self.n_samples: Optional[int] = None
+        self._dist: Optional[NDArray[np.floating]] = None
+        self._rng = np.random.default_rng()
+
+    # ------------------------------------------------------------------
+    # Private helpers
+    # ------------------------------------------------------------------
+
+    def _q_stat(self, idx_a: NDArray[np.intp], idx_b: NDArray[np.intp]) -> float:
+        """Return the scaled energy statistic between two index groups."""
+        m, n = idx_a.shape[0], idx_b.shape[0]
+        assert self._dist is not None  # guaranteed after fit()
+        cross = float(self._dist[np.ix_(idx_a, idx_b)].sum())
+        w_a = float(self._dist[np.ix_(idx_a, idx_a)].sum())
+        w_b = float(self._dist[np.ix_(idx_b, idx_b)].sum())
+        e = 2.0 * cross / (m * n) - w_a / (m * m) - w_b / (n * n)
+        return e * m * n / (m + n)
+
+    def _best_split(self, start: int, end: int) -> tuple[Optional[int], float]:
+        """Return ``(t*, Q*)`` — the optimal split and its Q value."""
+        best_t: Optional[int] = None
+        best_q = -np.inf
+        seg = np.arange(start, end, dtype=np.intp)
+        for t in range(start + self.min_size, end - self.min_size + 1):
+            m = t - start
+            q = self._q_stat(seg[:m], seg[m:])
+            if q > best_q:
+                best_q, best_t = q, t
+        return best_t, best_q
+
+    def _perm_test(
+        self,
+        start: int,
+        t_star: int,
+        end: int,
+        q_star: float,
+    ) -> bool:
+        """Return ``True`` if the split at ``t_star`` is significant.
+
+        Uses a fixed-split permutation test: the indices in ``[start, end)``
+        are shuffled and the energy statistic is recomputed at the same
+        split size as ``t_star``, yielding the null distribution.
+        """
+        if self.n_perms == 0:
+            return True
+        m = t_star - start
+        indices = np.arange(start, end, dtype=np.intp)
+        exceed = 0
+        for _ in range(self.n_perms):
+            perm = self._rng.permutation(indices)
+            if self._q_stat(perm[:m], perm[m:]) >= q_star:
+                exceed += 1
+        return (exceed / self.n_perms) < self.sig_level
+
+    # ------------------------------------------------------------------
+    # Public API
+    # ------------------------------------------------------------------
+
+    def fit(self, signal: NDArray[np.number]) -> Self:
+        """Pre-compute the pairwise distance matrix.
+
+        Args:
+            signal: array of shape ``(n_samples,)`` or
+                ``(n_samples, n_features)``.
+
+        Returns:
+            self
+        """
+        if signal.ndim == 1:
+            self.signal = signal.reshape(-1, 1).astype(float)
+        else:
+            self.signal = signal.astype(float)
+        self.n_samples = self.signal.shape[0]
+        self._dist = cdist(self.signal, self.signal, metric="euclidean") ** self.alpha
+        self._rng = np.random.default_rng()
+        return self
+
+    def predict(self, n_bkps: Optional[int] = None) -> list[int]:
+        """Return the list of change points.
+
+        Must be called after :meth:`fit`.
+
+        Args:
+            n_bkps: target number of change points.  When supplied,
+                splitting stops as soon as ``n_bkps`` breakpoints have
+                been found and no permutation test is performed.  When
+                ``None``, the permutation test (controlled by
+                ``n_perms`` and ``sig_level``) governs when to stop.
+
+        Returns:
+            Sorted list of breakpoint indices, by convention the last
+            element equals ``n_samples``.
+
+        Raises:
+            BadSegmentationParameters: if neither ``n_bkps`` nor the
+                permutation test can produce a valid segmentation.
+        """
+        if self.signal is None:
+            raise BadSegmentationParameters("Call fit() before predict().")
+
+        assert self.n_samples is not None
+
+        if n_bkps is not None and n_bkps < 0:
+            raise BadSegmentationParameters("n_bkps must be a non-negative integer.")
+
+        bkps: set[int] = set()
+        # Segments still eligible for splitting: list of (start, end) pairs.
+        pending: list[tuple[int, int]] = [(0, self.n_samples)]
+
+        while pending:
+            # Find the globally best split across all pending segments.
+            best_candidate: Optional[tuple[float, int, int, int]] = None
+            for start, end in pending:
+                if end - start < 2 * self.min_size:
+                    continue
+                t, q = self._best_split(start, end)
+                if t is not None:
+                    if best_candidate is None or q > best_candidate[0]:
+                        best_candidate = (q, t, start, end)
+
+            if best_candidate is None:
+                break
+
+            q_best, t_best, s_best, e_best = best_candidate
+
+            # Stopping criterion.
+            if n_bkps is not None:
+                if len(bkps) >= n_bkps:
+                    break
+            else:
+                if not self._perm_test(s_best, t_best, e_best, q_best):
+                    break
+
+            # Accept the split.
+            bkps.add(t_best)
+            pending.remove((s_best, e_best))
+            for child in ((s_best, t_best), (t_best, e_best)):
+                if child[1] - child[0] >= 2 * self.min_size:
+                    pending.append(child)
+
+        bkps.add(self.n_samples)
+        return sorted(bkps)
+
+    def fit_predict(
+        self,
+        signal: NDArray[np.number],
+        n_bkps: Optional[int] = None,
+    ) -> list[int]:
+        """Fit and predict in one step.
+
+        Args:
+            signal: array of shape ``(n_samples,)`` or
+                ``(n_samples, n_features)``.
+            n_bkps: see :meth:`predict`.
+
+        Returns:
+            see :meth:`predict`.
+        """
+        self.fit(signal)
+        return self.predict(n_bkps=n_bkps)
diff --git a/tests/test_edivisive.py b/tests/test_edivisive.py
new file mode 100644
index 00000000..8f675334
--- /dev/null
+++ b/tests/test_edivisive.py
@@ -0,0 +1,253 @@
+"""Tests for the EDivisive change point detection algorithm.
+
+Reference: Matteson & James (2014) "A Nonparametric Approach for
+Multiple Change Point Analysis of Multivariate Data", JASA
+109(505):334-345.
+"""
+
+import numpy as np
+import pytest
+from numpy.testing import assert_allclose
+
+from ruptures import EDivisive
+from ruptures.exceptions import BadSegmentationParameters
+
+
+# ---------------------------------------------------------------------------
+# Fixtures
+# ---------------------------------------------------------------------------
+
+
+@pytest.fixture(scope="module")
+def rng():
+    return np.random.default_rng(0)
+
+
+@pytest.fixture(scope="module")
+def signal_1d_two_bkps(rng):
+    """Piecewise-constant 1-D signal with two clear change points."""
+    s = np.concatenate(
+        [
+            rng.normal(0.0, 0.5, 100),
+            rng.normal(5.0, 0.5, 100),
+            rng.normal(0.0, 0.5, 100),
+        ]
+    )
+    return s, [100, 200, 300]
+
+
+@pytest.fixture(scope="module")
+def signal_5d_one_bkp(rng):
+    """Piecewise-constant 5-D signal with one clear change point."""
+    a = rng.normal(0.0, 0.5, (150, 5))
+    b = rng.normal(3.0, 0.5, (150, 5))
+    return np.vstack([a, b]), [150, 300]
+
+
+@pytest.fixture(scope="module")
+def constant_signal():
+    """Signal with no change point."""
+    return np.zeros(200)
+
+
+# ---------------------------------------------------------------------------
+# Interface / attribute tests
+# ---------------------------------------------------------------------------
+
+
+def test_fit_returns_self(signal_1d_two_bkps):
+    signal, _ = signal_1d_two_bkps
+    algo = EDivisive(n_perms=0)
+    assert algo.fit(signal) is algo
+
+
+def test_predict_last_element_equals_n_samples(signal_1d_two_bkps):
+    signal, _ = signal_1d_two_bkps
+    bkps = EDivisive(n_perms=0).fit_predict(signal, n_bkps=2)
+    assert bkps[-1] == signal.shape[0]
+
+
+def test_predict_n_bkps_length(signal_1d_two_bkps):
+    signal, _ = signal_1d_two_bkps
+    for k in [0, 1, 2]:
+        bkps = EDivisive(n_perms=0).fit_predict(signal, n_bkps=k)
+        assert len(bkps) == k + 1
+
+
+def test_predict_sorted(signal_1d_two_bkps):
+    signal, _ = signal_1d_two_bkps
+    bkps = EDivisive(n_perms=0).fit_predict(signal, n_bkps=2)
+    assert bkps == sorted(bkps)
+
+
+def test_predict_1d_input(signal_1d_two_bkps):
+    signal, _ = signal_1d_two_bkps
+    assert signal.ndim == 1
+    bkps = EDivisive(n_perms=0).fit_predict(signal, n_bkps=1)
+    assert bkps[-1] == signal.shape[0]
+
+
+def test_predict_2d_input(signal_5d_one_bkp):
+    signal, _ = signal_5d_one_bkp
+    bkps = EDivisive(n_perms=0).fit_predict(signal, n_bkps=1)
+    assert bkps[-1] == signal.shape[0]
+    assert len(bkps) == 2
+
+
+def test_fit_predict_equiv(signal_1d_two_bkps):
+    """fit().predict() and fit_predict() must agree."""
+    signal, _ = signal_1d_two_bkps
+    algo = EDivisive(n_perms=0)
+    a = algo.fit(signal).predict(n_bkps=2)
+    b = EDivisive(n_perms=0).fit_predict(signal, n_bkps=2)
+    assert a == b
+
+
+# ---------------------------------------------------------------------------
+# Error handling
+# ---------------------------------------------------------------------------
+
+
+def test_alpha_out_of_range():
+    with pytest.raises(ValueError, match="alpha"):
+        EDivisive(alpha=0.0)
+    with pytest.raises(ValueError, match="alpha"):
+        EDivisive(alpha=2.5)
+
+
+def test_predict_before_fit():
+    with pytest.raises(BadSegmentationParameters):
+        EDivisive().predict(n_bkps=1)
+
+
+def test_negative_n_bkps(signal_1d_two_bkps):
+    signal, _ = signal_1d_two_bkps
+    algo = EDivisive(n_perms=0).fit(signal)
+    with pytest.raises(BadSegmentationParameters):
+        algo.predict(n_bkps=-1)
+
+
+# ---------------------------------------------------------------------------
+# Detection accuracy
+# ---------------------------------------------------------------------------
+
+
+def test_detects_two_change_points_1d(signal_1d_two_bkps):
+    """Both breakpoints must be within ±15 samples of the true positions."""
+    signal, true_bkps = signal_1d_two_bkps
+    bkps = EDivisive(n_perms=0).fit_predict(signal, n_bkps=2)
+    # true_bkps without the final n_samples sentinel
+    detected = bkps[:-1]
+    expected = true_bkps[:-1]
+    for est, true in zip(sorted(detected), sorted(expected)):
+        assert abs(est - true) <= 15, f"Expected ~{true}, got {est}"
+
+
+def test_detects_one_change_point_5d(signal_5d_one_bkp):
+    """Breakpoint in a 5-D signal must be within ±15 samples."""
+    signal, true_bkps = signal_5d_one_bkp
+    bkps = EDivisive(n_perms=0).fit_predict(signal, n_bkps=1)
+    detected = bkps[0]
+    assert abs(detected - 150) <= 15, f"Expected ~150, got {detected}"
+
+
+def test_no_split_constant_signal_n_bkps_0(constant_signal):
+    """Zero breakpoints requested → only the sentinel."""
+    bkps = EDivisive(n_perms=0).fit_predict(constant_signal, n_bkps=0)
+    assert bkps == [200]
+
+
+# ---------------------------------------------------------------------------
+# Statistical properties of the energy statistic
+# ---------------------------------------------------------------------------
+
+
+def test_q_stat_increases_with_separation():
+    """Larger mean shift → larger energy statistic at the true split."""
+    rng = np.random.default_rng(1)
+    n = 100
+    results = {}
+    for delta in [0.5, 2.0, 5.0]:
+        signal = np.concatenate([rng.normal(0.0, 1.0, n), rng.normal(delta, 1.0, n)])
+        algo = EDivisive(n_perms=0).fit(signal)
+        _, q = algo._best_split(0, 2 * n)
+        results[delta] = q
+    assert results[0.5] < results[2.0] < results[5.0]
+
+
+def test_q_stat_zero_identical_segments():
+    """Energy statistic must be zero when both segments are identical."""
+    signal = np.ones(100)
+    algo = EDivisive(n_perms=0).fit(signal)
+    _, q = algo._best_split(0, 100)
+    assert abs(q) < 1e-10
+
+
+def test_best_split_exact_midpoint():
+    """Optimal split must be at the true midpoint for a step signal."""
+    signal = np.concatenate([np.zeros(50), np.full(50, 10.0)])
+    algo = EDivisive(n_perms=0).fit(signal)
+    t, _ = algo._best_split(0, 100)
+    assert t == 50
+
+
+# ---------------------------------------------------------------------------
+# alpha parameter
+# ---------------------------------------------------------------------------
+
+
+def test_alpha_2(signal_1d_two_bkps):
+    """Alpha=2 (squared distances) should still detect the change points."""
+    signal, _ = signal_1d_two_bkps
+    bkps = EDivisive(alpha=2.0, n_perms=0).fit_predict(signal, n_bkps=2)
+    assert len(bkps) == 3
+    assert bkps[-1] == signal.shape[0]
+
+
+def test_alpha_0_5(signal_1d_two_bkps):
+    """alpha=0.5 should still detect the change points."""
+    signal, _ = signal_1d_two_bkps
+    bkps = EDivisive(alpha=0.5, n_perms=0).fit_predict(signal, n_bkps=2)
+    assert len(bkps) == 3
+
+
+# ---------------------------------------------------------------------------
+# Permutation test
+# ---------------------------------------------------------------------------
+
+
+def test_perm_test_rejects_null_for_clear_change(rng):
+    """Permutation test must detect a 10-sigma shift at ~position 100."""
+    signal = np.concatenate(
+        [
+            rng.normal(0.0, 1.0, 100),
+            rng.normal(10.0, 1.0, 100),
+        ]
+    )
+    bkps = EDivisive(n_perms=200, sig_level=0.05).fit_predict(signal)
+    # The true change point at 100 must appear in the result (within ±5).
+    assert any(abs(b - 100) <= 5 for b in bkps), f"Breakpoints: {bkps}"
+    assert bkps[-1] == 200
+
+
+def test_perm_test_no_split_for_pure_noise(rng):
+    """Permutation test must not split a homogeneous noise signal."""
+    signal = rng.normal(0.0, 1.0, 200)
+    bkps = EDivisive(n_perms=200, sig_level=0.05).fit_predict(signal)
+    # May occasionally split due to random chance, but last element must be 200
+    assert bkps[-1] == 200
+
+
+# ---------------------------------------------------------------------------
+# min_size constraint
+# ---------------------------------------------------------------------------
+
+
+def test_min_size_respected(signal_1d_two_bkps):
+    """No segment in the result may be shorter than min_size."""
+    signal, _ = signal_1d_two_bkps
+    min_size = 20
+    bkps = EDivisive(min_size=min_size, n_perms=0).fit_predict(signal, n_bkps=2)
+    segments = list(zip([0] + bkps[:-1], bkps))
+    for start, end in segments:
+        assert end - start >= min_size, f"Segment [{start}, {end}) too short"