From c1c718bd8edc005e4a8874038d4d0c1ddddff65c Mon Sep 17 00:00:00 2001
From: Jules <jules@example.com>
Date: Wed, 24 Jun 2026 10:18:36 +0000
Subject: [PATCH] Centralize parity math utilities

This commit creates a new parity math utility module at src/eegprep/functions/miscfunc/parity.py, moving round_mat, rand_sample, and rand_permutation from local script scopes, and providing a parity_accumulate_float32 function to replace manual matrix multiplication logic in reref.py. The resample_raw and upfirdn_raw functions from pop_resample.py are also extracted to the parity library for centralization.
---
 src/eegprep/functions/miscfunc/misc.py        |  50 ---
 src/eegprep/functions/miscfunc/parity.py      | 315 ++++++++++++++++++
 src/eegprep/functions/popfunc/eeg_eegrej.py   |   2 +-
 .../functions/popfunc/eeg_point2lat.py        |   2 +-
 .../functions/popfunc/pop_load_frombids.py    |   3 +-
 src/eegprep/functions/popfunc/pop_resample.py | 158 +--------
 src/eegprep/functions/popfunc/pop_rmbase.py   |   2 +-
 src/eegprep/functions/popfunc/pop_topoplot.py |   2 +-
 src/eegprep/functions/sigprocfunc/eegrej.py   |   2 +-
 src/eegprep/functions/sigprocfunc/epoch.py    |   2 +-
 src/eegprep/functions/sigprocfunc/reref.py    |  12 +-
 src/eegprep/functions/sigprocfunc/runica.py   |   2 +-
 .../plugins/clean_rawdata/asr_calibrate.py    |   3 +-
 .../plugins/clean_rawdata/asr_process.py      |   3 +-
 .../plugins/clean_rawdata/clean_artifacts.py  |   2 +-
 .../plugins/clean_rawdata/clean_asr.py        |   2 +-
 .../plugins/clean_rawdata/clean_channels.py   |   3 +-
 .../plugins/clean_rawdata/clean_windows.py    |   2 +-
 .../plugins/clean_rawdata/private/ransac.py   |  90 +----
 .../plugins/clean_rawdata/private/stats.py    |   2 +-
 src/eegprep/plugins/firfilt/design.py         |   2 +-
 src/eegprep/plugins/firfilt/windows.py        |   2 +-
 tests/test_parity_rng.py                      |   4 +-
 tests/test_utils_ransac.py                    |   3 +-
 24 files changed, 346 insertions(+), 324 deletions(-)
 create mode 100644 src/eegprep/functions/miscfunc/parity.py

diff --git a/src/eegprep/functions/miscfunc/misc.py b/src/eegprep/functions/miscfunc/misc.py
index 3cd5b070..5797fdf2 100644
--- a/src/eegprep/functions/miscfunc/misc.py
+++ b/src/eegprep/functions/miscfunc/misc.py
@@ -1,6 +1,5 @@
 """Miscellaneous utility functions."""
 
-import math
 import sys
 import warnings
 from typing import Callable, Optional
@@ -15,7 +14,6 @@
     'num_cpus_from_reservation',
     'ToolError',
     'canonicalize_signs',
-    'round_mat',
     'aslist',
     'get_nested',
     'finite_matmul',
@@ -263,54 +261,6 @@ def canonicalize_signs(V):
     return V * sgn
 
 
-def round_mat(x, decimals=0):
-    """MATLAB-style rounding function.
-
-    - ties (.5 within fp error) round AWAY from zero
-    - supports positive/zero/negative `decimals` like MATLAB round(x, N)
-    - NaN/Inf propagate naturally
-    - does NOT return integer-typed results
-
-    This can be applied to numpy arrays and acts as a drop-in replacement
-    for np.round(), but also works for pure-Python float values; however,
-    to get a 1:1 replacement for a use of round(x) you need to write
-    int(round_mat(x)) since round() returns integers.
-
-    Parameters
-    ----------
-    x : array_like
-        The value(s) to round.
-    decimals : int
-        Number of decimals to round to.
-
-    Returns
-    -------
-    array_like
-        The rounded value(s).
-    """
-    if isinstance(x, (float, int)):
-        # Propagate NaN/Inf instead of throwing in math.floor(...)
-        if math.isnan(x) or math.isinf(x):
-            return x
-        xp = math
-    else:
-        xp = np
-        x = np.asarray(x)  # ensure ndarray
-
-    if decimals == 0:
-        return xp.copysign(xp.floor(abs(x) + 0.5), x)
-
-    if decimals > 0:
-        factor = 10.0**decimals
-        y = xp.copysign(xp.floor(abs(x) * factor + 0.5), x)
-        return y / factor
-
-    # decimals < 0  -> round to tens/hundreds/…
-    factor = 10.0 ** (-decimals)
-    y = xp.copysign(xp.floor(abs(x) / factor + 0.5), x)
-    return y * factor
-
-
 class SkippableException(Exception):
     """A dummy exception class for use in ExceptionUnlessDebug."""
 
diff --git a/src/eegprep/functions/miscfunc/parity.py b/src/eegprep/functions/miscfunc/parity.py
new file mode 100644
index 00000000..1792284d
--- /dev/null
+++ b/src/eegprep/functions/miscfunc/parity.py
@@ -0,0 +1,315 @@
+"""Parity Utility Library: Single source of truth for parity-critical math operations."""
+
+import math
+from math import ceil, floor, gcd
+import numpy as np
+from scipy import signal
+
+
+def round_mat(x, decimals=0):
+    """MATLAB-style rounding function.
+
+    - ties (.5 within fp error) round AWAY from zero
+    - supports positive/zero/negative `decimals` like MATLAB round(x, N)
+    - NaN/Inf propagate naturally
+    - does NOT return integer-typed results
+
+    This can be applied to numpy arrays and acts as a drop-in replacement
+    for np.round(), but also works for pure-Python float values; however,
+    to get a 1:1 replacement for a use of round(x) you need to write
+    int(round_mat(x)) since round() returns integers.
+
+    Parameters
+    ----------
+    x : array_like
+        The value(s) to round.
+    decimals : int
+        Number of decimals to round to.
+
+    Returns
+    -------
+    array_like
+        The rounded value(s).
+    """
+    if isinstance(x, (float, int)):
+        # Propagate NaN/Inf instead of throwing in math.floor(...)
+        if math.isnan(x) or math.isinf(x):
+            return x
+        xp = math
+    else:
+        xp = np
+        x = np.asarray(x)  # ensure ndarray
+
+    if decimals == 0:
+        return xp.copysign(xp.floor(abs(x) + 0.5), x)
+
+    if decimals > 0:
+        factor = 10.0**decimals
+        y = xp.copysign(xp.floor(abs(x) * factor + 0.5), x)
+        return y / factor
+
+    # decimals < 0  -> round to tens/hundreds/…
+    factor = 10.0 ** (-decimals)
+    y = xp.copysign(xp.floor(abs(x) / factor + 0.5), x)
+    return y * factor
+
+
+def rand_sample(n: int, m: int, stream: np.random.RandomState) -> np.ndarray:
+    """Random sampling without replacement using Fisher-Yates shuffle.
+
+    Optimized O(n) implementation using swap-based Fisher-Yates instead of
+    the previous O(n²) delete-based approach. Returns first m elements of
+    a random permutation of n items.
+
+    Args:
+        n: number of items to sample from
+        m: number of items to sample
+        stream: random number generator
+
+    Returns:
+        random_sample: array of m sampled values (indices 0..n-1)
+
+    Performance:
+        O(n) time complexity (was O(n²) in previous implementation)
+        For n=1M: ~3s (was ~80s) - 25x faster
+
+    Note:
+        This implementation uses Fisher-Yates shuffle for efficiency.
+        Results differ from the old O(n²) delete-based implementation,
+        but maintain parity with MATLAB's optimized rand_sample.
+    """
+    # Start with identity permutation
+    pool = np.arange(n)
+
+    # Fisher-Yates shuffle: only shuffle first m elements
+    for k in range(m):
+        # Choose from remaining elements (k to n-1)
+        remaining = n - k
+        choice = int(round_mat((remaining - 1) * stream.rand()))
+
+        # Swap pool[k] with pool[k + choice]
+        idx = k + choice
+        pool[k], pool[idx] = pool[idx], pool[k]
+
+    # Return first m elements
+    return pool[:m].copy()
+
+
+def rand_permutation(n: int, stream: np.random.RandomState) -> np.ndarray:
+    """Random permutation with MATLAB parity using Fisher-Yates shuffle.
+
+    This function produces the SAME permutation sequence as MATLAB's
+    rand_permutation() when both use the same RNG seed (5489). It achieves
+    parity by using rand() + round_mat() in a Fisher-Yates shuffle pattern
+    that matches MATLAB's implementation.
+
+    Optimized O(n) implementation (was O(n²) in previous version).
+
+    Args:
+        n: number of items to permute (returns permutation of 0..n-1)
+        stream: random number generator (np.random.RandomState)
+
+    Returns:
+        permutation: array of indices 0..n-1 in random order
+
+    Performance:
+        O(n) time complexity (was O(n²))
+        For n=1M: ~3s (was ~80s) - 25x faster
+
+    Example:
+        >>> rng = np.random.RandomState(5489)
+        >>> perm = rand_permutation(10, rng)
+        >>> # Matches MATLAB: rng(5489,'twister'); rand_permutation(10) - 1
+
+    Note:
+        This function is critical for ICA parity between Python and MATLAB.
+        Uses Fisher-Yates shuffle for O(n) performance.
+        Results differ from old O(n²) implementation but maintain
+        cross-platform parity with MATLAB.
+        See test_parity_rng.py for verification tests.
+    """
+    # Start with identity permutation [0, 1, 2, ..., n-1]
+    result = np.arange(n)
+
+    # Fisher-Yates shuffle: iterate backward from n-1 to 1
+    for k in range(n - 1, 0, -1):
+        # Pick random index from 0 to k (inclusive)
+        j = int(round_mat(k * stream.rand()))
+
+        # Swap elements k and j
+        result[k], result[j] = result[j], result[k]
+
+    return result
+
+
+def upfirdn_raw(x, h, p, q):
+    """Upfirdn implementation for resampling.
+
+    Parameters
+    ----------
+    x : array_like
+        Input signal.
+    h : array_like
+        Filter coefficients.
+    p : int
+        Upsampling factor.
+    q : int
+        Downsampling factor.
+
+    Returns
+    -------
+    y : ndarray
+        Filtered and resampled signal.
+    """
+    # Ensure x is a numpy array and h is 1D.
+    x = np.array(x, copy=True)
+    h = np.array(h).flatten()
+
+    # If x is a row vector, convert it to a column vector.
+    is_row_vector = False
+    if x.ndim == 2 and x.shape[0] == 1 and x.shape[1] > 1:
+        x = x.T
+        is_row_vector = True
+
+    rx, cx = x.shape
+    Lh = h.size
+    Ly = math.ceil(((rx - 1) * p + Lh) / q)
+    y = np.zeros((Ly, cx))
+
+    for c in range(cx):
+        for m in range(Ly):
+            n = (m * q) // p
+            lm = (m * q) % p
+            # k goes from max(0, n - rx + 1) to n (inclusive)
+            for k in range(max(0, n - rx + 1), n + 1):
+                if k * p + lm < Lh:
+                    y[m, c] += h[k * p + lm] * x[n - k, c]
+
+    if is_row_vector:
+        y = y.T
+
+    return y
+
+
+def resample_raw(x, p, q, h=None):
+    """Change the sample rate of x by a factor of p/q.
+
+    Parameters
+    ----------
+    x : array_like
+        The data to be resampled.
+    p : int
+        The upsampling factor.
+    q : int
+        The downsampling factor.
+    h : array_like, optional
+        The filter coefficients. If not provided, a Kaiser-windowed sinc filter is used.
+
+    Returns
+    -------
+    y : ndarray
+        The resampled array. If input is a vector, output will be a vector.
+    h : ndarray
+        The filter coefficients used.
+    """
+    # Input validation
+    if not isinstance(p, (int, np.integer)) or not isinstance(q, (int, np.integer)):
+        raise ValueError("p and q must be positive integers")
+    if p <= 0 or q <= 0:
+        raise ValueError("p and q must be positive integers")
+
+    # Convert x to numpy array and handle row vectors
+    x = np.asarray(x)
+    is_1d = x.ndim == 1
+
+    # Reshape input to 2D array with shape (samples, channels)
+    if is_1d:
+        x = x.reshape(-1, 1)
+    elif x.ndim == 2 and x.shape[0] == 1:
+        x = x.T
+
+    # Simplify decimation and interpolation factors
+    great_common_divisor = gcd(p, q)
+    if great_common_divisor > 1:
+        p = p // great_common_divisor
+        q = q // great_common_divisor
+
+    # Filter design if required
+    if h is None:
+        # Properties of the antialiasing filter
+        log10_rejection = -3.0
+        stopband_cutoff_f = 1.0 / (2.0 * max(p, q))
+        roll_off_width = stopband_cutoff_f / 10.0
+
+        # Determine filter length
+        rejection_dB = -20.0 * log10_rejection
+        L = ceil((rejection_dB - 8.0) / (28.714 * roll_off_width))
+
+        # Ideal sinc filter
+        t = np.arange(-L, L + 1)
+        ideal_filter = 2 * p * stopband_cutoff_f * np.sinc(2 * stopband_cutoff_f * t)
+
+        # Determine parameter of Kaiser window
+        if 21 <= rejection_dB <= 50:
+            beta = 0.5842 * (rejection_dB - 21.0) ** 0.4 + 0.07886 * (rejection_dB - 21.0)
+        elif rejection_dB > 50:
+            beta = 0.1102 * (rejection_dB - 8.7)
+        else:
+            beta = 0.0
+
+        # Apply Kaiser window to ideal filter
+        h = ideal_filter * signal.windows.kaiser(2 * L + 1, beta)
+
+    if not np.isrealobj(h):
+        raise ValueError("The filter h should be a real vector")
+
+    h = np.asarray(h)
+    if h.ndim != 1:
+        raise ValueError("The filter h should be a vector")
+
+    Lx = x.shape[0]
+    Lh = len(h)
+    L = (Lh - 1) / 2.0
+    Ly = ceil(Lx * p / q)
+
+    # Pre and postpad filter response
+    nz_pre = floor(q - np.mod(L, q))
+    h_padded = np.pad(h, (nz_pre, 0), 'constant')
+
+    offset = floor((L + nz_pre) / q)
+    nz_post = 0
+    while ceil(((Lx - 1) * p + nz_pre + Lh + nz_post) / q) - offset < Ly:
+        nz_post += 1
+    h_padded = np.pad(h_padded, (0, nz_post), 'constant')
+
+    # Filtering - fixed upfirdn usage
+    y = upfirdn_raw(x, h_padded, p, q)
+    y = y[offset : offset + Ly]
+
+    # Restore original dimensionality
+    if is_1d:
+        y = y.flatten()
+    else:
+        y = y.reshape(-1, x.shape[1])
+
+    return y, h
+
+
+def parity_accumulate_float32(matrix, data):
+    """
+    Multiply matrix @ data with MATLAB-compatible float32 accumulation.
+
+    MATLAB accumulates this product column-major; replicate that float32
+    accumulation order in NumPy (row-major) by transposing the operands.
+    Algebraically, (data.T @ matrix.T).T == matrix @ data.
+
+    Args:
+        matrix: 2D numpy array (e.g., refmatrix)
+        data: 2D numpy array (channels x points)
+
+    Returns:
+        The matrix product with MATLAB-parity accumulation order.
+    """
+    dt = matrix.dtype
+    block = np.ascontiguousarray(data.astype(dt).T)
+    return (block @ matrix.T).T
diff --git a/src/eegprep/functions/popfunc/eeg_eegrej.py b/src/eegprep/functions/popfunc/eeg_eegrej.py
index 1ea0c289..6fbc6b32 100644
--- a/src/eegprep/functions/popfunc/eeg_eegrej.py
+++ b/src/eegprep/functions/popfunc/eeg_eegrej.py
@@ -6,7 +6,7 @@
 from copy import deepcopy
 from eegprep.functions.miscfunc.event_utils import boundary_event_indices
 from eegprep.functions.miscfunc.event_utils import is_boundary_event as _is_boundary_event
-from ..miscfunc.misc import round_mat
+from ..miscfunc.parity import round_mat
 
 
 logger = logging.getLogger(__name__)
diff --git a/src/eegprep/functions/popfunc/eeg_point2lat.py b/src/eegprep/functions/popfunc/eeg_point2lat.py
index fdb75955..1c7aeb0f 100644
--- a/src/eegprep/functions/popfunc/eeg_point2lat.py
+++ b/src/eegprep/functions/popfunc/eeg_point2lat.py
@@ -1,7 +1,7 @@
 """Module for converting event latencies from points to time units."""
 
 import numpy as np
-from ..miscfunc.misc import round_mat
+from ..miscfunc.parity import round_mat
 
 
 def eeg_point2lat(lat_array, epoch_array=None, srate=None, timewin=None, timeunit=1.0):
diff --git a/src/eegprep/functions/popfunc/pop_load_frombids.py b/src/eegprep/functions/popfunc/pop_load_frombids.py
index 6373d132..b4f66087 100644
--- a/src/eegprep/functions/popfunc/pop_load_frombids.py
+++ b/src/eegprep/functions/popfunc/pop_load_frombids.py
@@ -16,7 +16,8 @@
 )
 from eegprep.plugins.EEG_BIDS.montage import apply_montage_inference
 from eegprep.plugins.EEG_BIDS.raw import load_raw_eeg_file
-from eegprep.functions.miscfunc.misc import ExceptionUnlessDebug, round_mat
+from eegprep.functions.miscfunc.misc import ExceptionUnlessDebug
+from eegprep.functions.miscfunc.parity import round_mat
 
 import numpy as np
 
diff --git a/src/eegprep/functions/popfunc/pop_resample.py b/src/eegprep/functions/popfunc/pop_resample.py
index 75a64df1..476155d7 100644
--- a/src/eegprep/functions/popfunc/pop_resample.py
+++ b/src/eegprep/functions/popfunc/pop_resample.py
@@ -2,12 +2,10 @@
 
 from copy import deepcopy
 import logging
-import math
-from math import ceil, floor, gcd
+from math import ceil
 
 import numpy as np
 import sympy as sp
-from scipy import signal
 from scipy.signal import resample, resample_poly
 from scipy.signal.windows import kaiser
 
@@ -15,6 +13,7 @@
 from eegprep.functions.guifunc.inputgui import inputgui
 from eegprep.functions.guifunc.spec import CallbackSpec, ControlSpec, DialogSpec
 from eegprep.functions.miscfunc.event_utils import is_boundary_event as _shared_is_boundary_event
+from eegprep.functions.miscfunc.parity import resample_raw
 from eegprep.functions.popfunc._file_io import events_to_records
 from eegprep.plugins.firfilt.firws import firws
 from eegprep.plugins.firfilt.firwsord import firwsord
@@ -319,156 +318,3 @@ def _scale_duration(event, ratio):
     if "duration" not in event or event["duration"] in (None, ""):
         return
     event["duration"] = float(event["duration"]) * ratio
-
-
-def upfirdn_raw(x, h, p, q):
-    """Upfirdn implementation for resampling.
-
-    Parameters
-    ----------
-    x : array_like
-        Input signal.
-    h : array_like
-        Filter coefficients.
-    p : int
-        Upsampling factor.
-    q : int
-        Downsampling factor.
-
-    Returns
-    -------
-    y : ndarray
-        Filtered and resampled signal.
-    """
-    # Ensure x is a numpy array and h is 1D.
-    x = np.array(x, copy=True)
-    h = np.array(h).flatten()
-
-    # If x is a row vector, convert it to a column vector.
-    is_row_vector = False
-    if x.ndim == 2 and x.shape[0] == 1 and x.shape[1] > 1:
-        x = x.T
-        is_row_vector = True
-
-    rx, cx = x.shape
-    Lh = h.size
-    Ly = math.ceil(((rx - 1) * p + Lh) / q)
-    y = np.zeros((Ly, cx))
-
-    for c in range(cx):
-        for m in range(Ly):
-            n = (m * q) // p
-            lm = (m * q) % p
-            # k goes from max(0, n - rx + 1) to n (inclusive)
-            for k in range(max(0, n - rx + 1), n + 1):
-                if k * p + lm < Lh:
-                    y[m, c] += h[k * p + lm] * x[n - k, c]
-
-    if is_row_vector:
-        y = y.T
-
-    return y
-
-
-def resample_raw(x, p, q, h=None):
-    """Change the sample rate of x by a factor of p/q.
-
-    Parameters
-    ----------
-    x : array_like
-        The data to be resampled.
-    p : int
-        The upsampling factor.
-    q : int
-        The downsampling factor.
-    h : array_like, optional
-        The filter coefficients. If not provided, a Kaiser-windowed sinc filter is used.
-
-    Returns
-    -------
-    y : ndarray
-        The resampled array. If input is a vector, output will be a vector.
-    h : ndarray
-        The filter coefficients used.
-    """
-    # Input validation
-    if not isinstance(p, (int, np.integer)) or not isinstance(q, (int, np.integer)):
-        raise ValueError("p and q must be positive integers")
-    if p <= 0 or q <= 0:
-        raise ValueError("p and q must be positive integers")
-
-    # Convert x to numpy array and handle row vectors
-    x = np.asarray(x)
-    is_1d = x.ndim == 1
-
-    # Reshape input to 2D array with shape (samples, channels)
-    if is_1d:
-        x = x.reshape(-1, 1)
-    elif x.ndim == 2 and x.shape[0] == 1:
-        x = x.T
-
-    # Simplify decimation and interpolation factors
-    great_common_divisor = gcd(p, q)
-    if great_common_divisor > 1:
-        p = p // great_common_divisor
-        q = q // great_common_divisor
-
-    # Filter design if required
-    if h is None:
-        # Properties of the antialiasing filter
-        log10_rejection = -3.0
-        stopband_cutoff_f = 1.0 / (2.0 * max(p, q))
-        roll_off_width = stopband_cutoff_f / 10.0
-
-        # Determine filter length
-        rejection_dB = -20.0 * log10_rejection
-        L = ceil((rejection_dB - 8.0) / (28.714 * roll_off_width))
-
-        # Ideal sinc filter
-        t = np.arange(-L, L + 1)
-        ideal_filter = 2 * p * stopband_cutoff_f * np.sinc(2 * stopband_cutoff_f * t)
-
-        # Determine parameter of Kaiser window
-        if 21 <= rejection_dB <= 50:
-            beta = 0.5842 * (rejection_dB - 21.0) ** 0.4 + 0.07886 * (rejection_dB - 21.0)
-        elif rejection_dB > 50:
-            beta = 0.1102 * (rejection_dB - 8.7)
-        else:
-            beta = 0.0
-
-        # Apply Kaiser window to ideal filter
-        h = ideal_filter * signal.windows.kaiser(2 * L + 1, beta)
-
-    if not np.isrealobj(h):
-        raise ValueError("The filter h should be a real vector")
-
-    h = np.asarray(h)
-    if h.ndim != 1:
-        raise ValueError("The filter h should be a vector")
-
-    Lx = x.shape[0]
-    Lh = len(h)
-    L = (Lh - 1) / 2.0
-    Ly = ceil(Lx * p / q)
-
-    # Pre and postpad filter response
-    nz_pre = floor(q - np.mod(L, q))
-    h_padded = np.pad(h, (nz_pre, 0), 'constant')
-
-    offset = floor((L + nz_pre) / q)
-    nz_post = 0
-    while ceil(((Lx - 1) * p + nz_pre + Lh + nz_post) / q) - offset < Ly:
-        nz_post += 1
-    h_padded = np.pad(h_padded, (0, nz_post), 'constant')
-
-    # Filtering - fixed upfirdn usage
-    y = upfirdn_raw(x, h_padded, p, q)
-    y = y[offset : offset + Ly]
-
-    # Restore original dimensionality
-    if is_1d:
-        y = y.flatten()
-    else:
-        y = y.reshape(-1, x.shape[1])
-
-    return y, h
diff --git a/src/eegprep/functions/popfunc/pop_rmbase.py b/src/eegprep/functions/popfunc/pop_rmbase.py
index 8808658f..b059fa98 100644
--- a/src/eegprep/functions/popfunc/pop_rmbase.py
+++ b/src/eegprep/functions/popfunc/pop_rmbase.py
@@ -10,7 +10,7 @@
 
 from eegprep.functions.guifunc.inputgui import inputgui
 from eegprep.functions.guifunc.spec import CallbackSpec, ControlSpec, DialogSpec
-from eegprep.functions.miscfunc.misc import round_mat
+from eegprep.functions.miscfunc.parity import round_mat
 from eegprep.functions.popfunc._chanutils import chanlocs_as_list
 from eegprep.functions.popfunc._pop_utils import format_history_value, parse_text_tokens
 from eegprep.functions.popfunc.eeg_findboundaries import eeg_findboundaries
diff --git a/src/eegprep/functions/popfunc/pop_topoplot.py b/src/eegprep/functions/popfunc/pop_topoplot.py
index 40585c80..89a21038 100644
--- a/src/eegprep/functions/popfunc/pop_topoplot.py
+++ b/src/eegprep/functions/popfunc/pop_topoplot.py
@@ -11,7 +11,7 @@
 
 from eegprep.functions.guifunc.inputgui import inputgui
 from eegprep.functions.guifunc.spec import ControlSpec, DialogSpec
-from eegprep.functions.miscfunc.misc import round_mat
+from eegprep.functions.miscfunc.parity import round_mat
 from eegprep.functions.popfunc._chanutils import chanlocs_as_list
 from eegprep.functions.popfunc._plot_utils import component_map_data
 from eegprep.functions.popfunc._plot_utils import history_command as plot_history_command
diff --git a/src/eegprep/functions/sigprocfunc/eegrej.py b/src/eegprep/functions/sigprocfunc/eegrej.py
index c6aadb4e..efd3e812 100644
--- a/src/eegprep/functions/sigprocfunc/eegrej.py
+++ b/src/eegprep/functions/sigprocfunc/eegrej.py
@@ -3,7 +3,7 @@
 import numpy as np
 from typing import List, Dict, Optional, Tuple
 from eegprep.functions.miscfunc.event_utils import is_boundary_event as _is_boundary_event
-from ..miscfunc.misc import round_mat
+from ..miscfunc.parity import round_mat
 
 
 def eegrej(
diff --git a/src/eegprep/functions/sigprocfunc/epoch.py b/src/eegprep/functions/sigprocfunc/epoch.py
index 421d593d..8ec9e48c 100644
--- a/src/eegprep/functions/sigprocfunc/epoch.py
+++ b/src/eegprep/functions/sigprocfunc/epoch.py
@@ -8,7 +8,7 @@
 
 import numpy as np
 
-from ..miscfunc.misc import round_mat
+from ..miscfunc.parity import round_mat
 
 logger = logging.getLogger(__name__)
 
diff --git a/src/eegprep/functions/sigprocfunc/reref.py b/src/eegprep/functions/sigprocfunc/reref.py
index 2edf4d44..a5978bef 100644
--- a/src/eegprep/functions/sigprocfunc/reref.py
+++ b/src/eegprep/functions/sigprocfunc/reref.py
@@ -8,6 +8,7 @@
 import numpy as np
 
 from eegprep.functions.popfunc._chanutils import normalise_reflocs as _normalise_reflocs
+from eegprep.functions.miscfunc.parity import parity_accumulate_float32
 
 
 def reref(
@@ -80,18 +81,11 @@ def reref(
     else:
         # Average reference via matrix multiplication, matching MATLAB's
         # reref.m: refmatrix = eye(n) - ones(n)/n; data = refmatrix * data.
-        # MATLAB accumulates this product column-major; replicate that float32
-        # accumulation order in NumPy (row-major) by transposing the operands
-        # (refmatrix is symmetric, so (data.T @ refmatrix).T == refmatrix @ data
-        # algebraically). The plain np.dot(refmatrix, data) form is ~0.06 uV off
-        # MATLAB in float32, which cascades through nonlinear downstream steps
-        # (ICA convergence) for borderline subjects; the transposed form is
-        # bit-exact to MATLAB on x86 BLAS.
+        # Uses parity_accumulate_float32 to replicate column-major float32 accumulation.
         n = len(chansin)
         dt = original_dtype
         refmatrix = np.eye(n, dtype=dt) - np.ones((n, n), dtype=dt) / dt.type(n)
-        block = np.ascontiguousarray(work[chansin_array, :].astype(dt).T)
-        work[chansin_array, :] = (block @ refmatrix).T
+        work[chansin_array, :] = parity_accumulate_float32(refmatrix, work[chansin_array, :])
         mean_data = None
 
     if locs is not None:
diff --git a/src/eegprep/functions/sigprocfunc/runica.py b/src/eegprep/functions/sigprocfunc/runica.py
index 5c23534b..470699d5 100644
--- a/src/eegprep/functions/sigprocfunc/runica.py
+++ b/src/eegprep/functions/sigprocfunc/runica.py
@@ -24,7 +24,7 @@
 
 import numpy as np
 from scipy.linalg import sqrtm, pinv, eig
-from ...plugins.clean_rawdata.private.ransac import rand_permutation
+from ..miscfunc.parity import rand_permutation
 from ..miscfunc.misc import finite_pinv
 
 logger = logging.getLogger(__name__)
diff --git a/src/eegprep/plugins/clean_rawdata/asr_calibrate.py b/src/eegprep/plugins/clean_rawdata/asr_calibrate.py
index 775aaddc..fa36e7c8 100644
--- a/src/eegprep/plugins/clean_rawdata/asr_calibrate.py
+++ b/src/eegprep/plugins/clean_rawdata/asr_calibrate.py
@@ -6,7 +6,8 @@
 import scipy.signal
 import scipy.linalg
 
-from ...functions.miscfunc.misc import canonicalize_signs, finite_matmul, round_mat
+from ...functions.miscfunc.misc import canonicalize_signs, finite_matmul
+from ...functions.miscfunc.parity import round_mat
 from .private.covariance import cov_mean, cov_shrinkage
 from .private.stats import fit_eeg_distribution, geometric_median
 
diff --git a/src/eegprep/plugins/clean_rawdata/asr_process.py b/src/eegprep/plugins/clean_rawdata/asr_process.py
index c77b0134..a68b117e 100644
--- a/src/eegprep/plugins/clean_rawdata/asr_process.py
+++ b/src/eegprep/plugins/clean_rawdata/asr_process.py
@@ -5,7 +5,8 @@
 import numpy as np
 import scipy.signal
 
-from ...functions.miscfunc.misc import finite_matmul, finite_pinv, round_mat
+from ...functions.miscfunc.misc import finite_matmul, finite_pinv
+from ...functions.miscfunc.parity import round_mat
 from .private.sigproc import moving_average
 
 logger = logging.getLogger(__name__)
diff --git a/src/eegprep/plugins/clean_rawdata/clean_artifacts.py b/src/eegprep/plugins/clean_rawdata/clean_artifacts.py
index 6b0e6cd3..692a42e0 100644
--- a/src/eegprep/plugins/clean_rawdata/clean_artifacts.py
+++ b/src/eegprep/plugins/clean_rawdata/clean_artifacts.py
@@ -14,7 +14,7 @@
 from .clean_asr import clean_asr
 from .clean_windows import clean_windows
 from .private.masks import mask_to_intervals
-from ...functions.miscfunc.misc import round_mat
+from ...functions.miscfunc.parity import round_mat
 from ...functions.popfunc.eeg_eegrej import eeg_eegrej
 
 
diff --git a/src/eegprep/plugins/clean_rawdata/clean_asr.py b/src/eegprep/plugins/clean_rawdata/clean_asr.py
index 76470820..a60e5e8b 100644
--- a/src/eegprep/plugins/clean_rawdata/clean_asr.py
+++ b/src/eegprep/plugins/clean_rawdata/clean_asr.py
@@ -15,7 +15,7 @@
 from .asr_calibrate import asr_calibrate
 from .asr_process import asr_process
 from .clean_windows import clean_windows
-from ...functions.miscfunc.misc import round_mat
+from ...functions.miscfunc.parity import round_mat
 
 logger = logging.getLogger(__name__)
 
diff --git a/src/eegprep/plugins/clean_rawdata/clean_channels.py b/src/eegprep/plugins/clean_rawdata/clean_channels.py
index 30b05bb7..dbfaf2a0 100644
--- a/src/eegprep/plugins/clean_rawdata/clean_channels.py
+++ b/src/eegprep/plugins/clean_rawdata/clean_channels.py
@@ -7,7 +7,8 @@
 
 from eegprep.plugins.firfilt.design import design_fir
 
-from ...functions.miscfunc.misc import finite_matmul, round_mat
+from ...functions.miscfunc.misc import finite_matmul
+from ...functions.miscfunc.parity import round_mat
 from .private.channel_removal import remove_channels_without_pop_select, update_clean_channel_mask
 from .private.ransac import calc_projector
 from .private.sigproc import filtfilt_fast
diff --git a/src/eegprep/plugins/clean_rawdata/clean_windows.py b/src/eegprep/plugins/clean_rawdata/clean_windows.py
index 732d17dd..23318342 100644
--- a/src/eegprep/plugins/clean_rawdata/clean_windows.py
+++ b/src/eegprep/plugins/clean_rawdata/clean_windows.py
@@ -10,7 +10,7 @@
 
 import numpy as np
 
-from ...functions.miscfunc.misc import round_mat
+from ...functions.miscfunc.parity import round_mat
 from ...functions.popfunc.eeg_eegrej import eeg_eegrej
 from .private.masks import mask_to_intervals
 from .private.stats import fit_eeg_distribution
diff --git a/src/eegprep/plugins/clean_rawdata/private/ransac.py b/src/eegprep/plugins/clean_rawdata/private/ransac.py
index 1a80bee1..d5425691 100644
--- a/src/eegprep/plugins/clean_rawdata/private/ransac.py
+++ b/src/eegprep/plugins/clean_rawdata/private/ransac.py
@@ -5,98 +5,10 @@
 import numpy as np
 
 from ....functions.adminfunc.eeglabcompat import get_eeglab
-from ....functions.miscfunc.misc import round_mat
+from ....functions.miscfunc.parity import rand_sample
 from .sphericalSplineInterpolate import sphericalSplineInterpolate
 
 
-def rand_sample(n: int, m: int, stream: np.random.RandomState) -> np.ndarray:
-    """Random sampling without replacement using Fisher-Yates shuffle.
-
-    Optimized O(n) implementation using swap-based Fisher-Yates instead of
-    the previous O(n²) delete-based approach. Returns first m elements of
-    a random permutation of n items.
-
-    Args:
-        n: number of items to sample from
-        m: number of items to sample
-        stream: random number generator
-
-    Returns:
-        random_sample: array of m sampled values (indices 0..n-1)
-
-    Performance:
-        O(n) time complexity (was O(n²) in previous implementation)
-        For n=1M: ~3s (was ~80s) - 25x faster
-
-    Note:
-        This implementation uses Fisher-Yates shuffle for efficiency.
-        Results differ from the old O(n²) delete-based implementation,
-        but maintain parity with MATLAB's optimized rand_sample.
-    """
-    # Start with identity permutation
-    pool = np.arange(n)
-
-    # Fisher-Yates shuffle: only shuffle first m elements
-    for k in range(m):
-        # Choose from remaining elements (k to n-1)
-        remaining = n - k
-        choice = int(round_mat((remaining - 1) * stream.rand()))
-
-        # Swap pool[k] with pool[k + choice]
-        idx = k + choice
-        pool[k], pool[idx] = pool[idx], pool[k]
-
-    # Return first m elements
-    return pool[:m].copy()
-
-
-def rand_permutation(n: int, stream: np.random.RandomState) -> np.ndarray:
-    """Random permutation with MATLAB parity using Fisher-Yates shuffle.
-
-    This function produces the SAME permutation sequence as MATLAB's
-    rand_permutation() when both use the same RNG seed (5489). It achieves
-    parity by using rand() + round_mat() in a Fisher-Yates shuffle pattern
-    that matches MATLAB's implementation.
-
-    Optimized O(n) implementation (was O(n²) in previous version).
-
-    Args:
-        n: number of items to permute (returns permutation of 0..n-1)
-        stream: random number generator (np.random.RandomState)
-
-    Returns:
-        permutation: array of indices 0..n-1 in random order
-
-    Performance:
-        O(n) time complexity (was O(n²))
-        For n=1M: ~3s (was ~80s) - 25x faster
-
-    Example:
-        >>> rng = np.random.RandomState(5489)
-        >>> perm = rand_permutation(10, rng)
-        >>> # Matches MATLAB: rng(5489,'twister'); rand_permutation(10) - 1
-
-    Note:
-        This function is critical for ICA parity between Python and MATLAB.
-        Uses Fisher-Yates shuffle for O(n) performance.
-        Results differ from old O(n²) implementation but maintain
-        cross-platform parity with MATLAB.
-        See test_parity_rng.py for verification tests.
-    """
-    # Start with identity permutation [0, 1, 2, ..., n-1]
-    result = np.arange(n)
-
-    # Fisher-Yates shuffle: iterate backward from n-1 to 1
-    for k in range(n - 1, 0, -1):
-        # Pick random index from 0 to k (inclusive)
-        j = int(round_mat(k * stream.rand()))
-
-        # Swap elements k and j
-        result[k], result[j] = result[j], result[k]
-
-    return result
-
-
 def calc_projector(
     locs: np.ndarray,
     num_samples: int,
diff --git a/src/eegprep/plugins/clean_rawdata/private/stats.py b/src/eegprep/plugins/clean_rawdata/private/stats.py
index 709a1855..5bd6a18f 100644
--- a/src/eegprep/plugins/clean_rawdata/private/stats.py
+++ b/src/eegprep/plugins/clean_rawdata/private/stats.py
@@ -5,7 +5,7 @@
 import numpy as np
 from numpy.linalg import norm as np_norm  # Use alias to avoid potential name collision
 from scipy.special import gamma, gammaincinv
-from ....functions.miscfunc.misc import round_mat
+from ....functions.miscfunc.parity import round_mat
 
 logger = logging.getLogger(__name__)
 
diff --git a/src/eegprep/plugins/firfilt/design.py b/src/eegprep/plugins/firfilt/design.py
index 7bd3e9af..89058b74 100644
--- a/src/eegprep/plugins/firfilt/design.py
+++ b/src/eegprep/plugins/firfilt/design.py
@@ -6,7 +6,7 @@
 
 import numpy as np
 
-from eegprep.functions.miscfunc.misc import round_mat
+from eegprep.functions.miscfunc.parity import round_mat
 
 __all__ = ["design_fir", "design_kaiser"]
 
diff --git a/src/eegprep/plugins/firfilt/windows.py b/src/eegprep/plugins/firfilt/windows.py
index 0a2cde48..77515809 100644
--- a/src/eegprep/plugins/firfilt/windows.py
+++ b/src/eegprep/plugins/firfilt/windows.py
@@ -7,7 +7,7 @@
 import numpy as np
 from scipy.special import i0
 
-from eegprep.functions.miscfunc.misc import round_mat
+from eegprep.functions.miscfunc.parity import round_mat
 
 
 def windows(t: str, m: float, a: float | None = None) -> np.ndarray:
diff --git a/tests/test_parity_rng.py b/tests/test_parity_rng.py
index 6d57ea54..1b0c33b5 100644
--- a/tests/test_parity_rng.py
+++ b/tests/test_parity_rng.py
@@ -18,8 +18,8 @@
 import tempfile
 import scipy.io
 from eegprep.functions.adminfunc.eeglabcompat import get_eeglab
-from eegprep.plugins.clean_rawdata.private.ransac import rand_sample
-from eegprep.functions.miscfunc.misc import round_mat
+from eegprep.functions.miscfunc.parity import rand_sample
+from eegprep.functions.miscfunc.parity import round_mat
 
 
 class TestRNGParity(unittest.TestCase):
diff --git a/tests/test_utils_ransac.py b/tests/test_utils_ransac.py
index 3e386ece..20cc0a87 100644
--- a/tests/test_utils_ransac.py
+++ b/tests/test_utils_ransac.py
@@ -2,7 +2,8 @@
 import numpy as np
 from unittest.mock import patch, MagicMock
 
-from eegprep.plugins.clean_rawdata.private.ransac import rand_sample, calc_projector
+from eegprep.functions.miscfunc.parity import rand_sample
+from eegprep.plugins.clean_rawdata.private.ransac import calc_projector
 from eegprep.plugins.clean_rawdata.private.sphericalSplineInterpolate import sphericalSplineInterpolate