Repository: simple-crypto/SCALib (language: python)

Issue:
Feature request: Multivariate one-pass t-test
Hi,

I would like to ask if you could add the multivariate version of the one-pass t-test from [1] (Section 5). The current univariate version is only applicable to hardware implementations that calculate on shares in parallel while the multivariate version would be more suitable for evaluating masked software on microprocessors or serial masked hardware designs.

[1] https://eprint.iacr.org/2015/207.pdf



Discussion / hints:

Fix patch:
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index ca1fd17f..7273511f 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -6,6 +6,8 @@ Not released
 ------------
 
 * SASCA: support modular ADD and MUL operations (#18)
+* TTest: Performance improvement by using a mix of 2 passes and 1 pass algorithms 
+* MTTest: First implementation of multivariate T-test.
 * Improved documentation and README.rst
 
 v0.3.4
diff --git a/README.rst b/README.rst
index e6a0d30b..41e737c8 100644
--- a/README.rst
+++ b/README.rst
@@ -19,14 +19,17 @@ SCALib contains various features for side-channel analysis, see the documentatio
 Install
 =======
 
-SCALib is on PyPi! See install_ for details.
+SCALib is on PyPi!
 
 .. code-block::
 
     pip install scalib
 
+See install_ for details.
+ 
 .. _install: https://scalib.readthedocs.io/en/latest/index.html#install
 
+
 About us
 ========
 SCALib has been initiated by Olivier Bronchain during his PhD at Crypto Group,
diff --git a/src/scalib/attacks/sascagraph.py b/src/scalib/attacks/sascagraph.py
index 112985e9..4d641793 100644
--- a/src/scalib/attacks/sascagraph.py
+++ b/src/scalib/attacks/sascagraph.py
@@ -48,7 +48,7 @@ class SASCAGraph:
             VAR MULTI x # Sbox input
             VAR MULTI y # Sbox output (whose leakage is targeted)
             PROPERTY x = k ^ p   # Key addition
-            PROPERTY x = sbox[y] # Sbox lookup
+            PROPERTY y = sbox[x] # Sbox lookup
             '''
         # n is the number of traces for our attack.
         graph = SASCAGraph(graph_desc,n)
diff --git a/src/scalib/metrics/__init__.py b/src/scalib/metrics/__init__.py
index f6902158..312a5010 100644
--- a/src/scalib/metrics/__init__.py
+++ b/src/scalib/metrics/__init__.py
@@ -1,2 +1,3 @@
 from .snr import SNR
 from .ttest import Ttest
+from .ttest import MTtest
diff --git a/src/scalib_ext/scalib-py/src/lib.rs b/src/scalib_ext/scalib-py/src/lib.rs
index 9fc5c108..c53b2fc5 100644
--- a/src/scalib_ext/scalib-py/src/lib.rs
+++ b/src/scalib_ext/scalib-py/src/lib.rs
@@ -31,6 +31,7 @@ fn str2method(s: &str) -> Result<ranklib::RankingMethod, &str> {
 fn _scalib_ext(_py: Python, m: &PyModule) -> PyResult<()> {
     m.add_class::<snr::SNR>()?;
     m.add_class::<ttest::Ttest>()?;
+    m.add_class::<ttest::MTtest>()?;
     m.add_class::<lda::LDA>()?;
     m.add_class::<lda::LdaAcc>()?;
 
diff --git a/src/scalib_ext/scalib/Cargo.toml b/src/scalib_ext/scalib/Cargo.toml
index a68589f6..00eec3f2 100644
--- a/src/scalib_ext/scalib/Cargo.toml
+++ b/src/scalib_ext/scalib/Cargo.toml
@@ -35,9 +35,6 @@ rustfft = "6.0.0"
 num-complex = "0.4.0"
 
 
-[profile.release]
-debug = true
-
 [dev-dependencies]
 criterion = "0.3.4"
 ndarray-rand = "0.14.0"
@@ -53,3 +50,11 @@ harness = false
 [[bench]]
 name = "snr_update"
 harness = false
+
+[[bench]]
+name = "mttest"
+harness = false
+
+[[bench]]
+name = "ttest"
+harness = false
diff --git a/src/scalib_ext/scalib/src/belief_propagation.rs b/src/scalib_ext/scalib/src/belief_propagation.rs
index f83c4c58..dd0131d9 100644
--- a/src/scalib_ext/scalib/src/belief_propagation.rs
+++ b/src/scalib_ext/scalib/src/belief_propagation.rs
@@ -269,12 +269,12 @@ pub fn update_functions(functions: &[Func], edges: &mut [Vec<&mut Array2<f64>>])
                             let value = value.first().unwrap();
                             for i1 in 0..nc {
                                 let o: usize = match &function.functype {
-                                    FuncType::XORCST(values) => ((i1 as u32) ^ value) as usize,
-                                    FuncType::ANDCST(values) => ((i1 as u32) & value) as usize,
-                                    FuncType::ADDCST(values) => {
+                                    FuncType::XORCST(_) => ((i1 as u32) ^ value) as usize,
+                                    FuncType::ANDCST(_) => ((i1 as u32) & value) as usize,
+                                    FuncType::ADDCST(_) => {
                                         (((i1 as u32) + value) % (nc as u32)) as usize
                                     }
-                                    FuncType::MULCST(values) => {
+                                    FuncType::MULCST(_) => {
                                         (((i1 as u32) * value) % (nc as u32)) as usize
                                     }
                                     _ => unreachable!(),
diff --git a/src/scalib_ext/scalib/src/lib.rs b/src/scalib_ext/scalib/src/lib.rs
index f79a4520..a2aeb62e 100644
--- a/src/scalib_ext/scalib/src/lib.rs
+++ b/src/scalib_ext/scalib/src/lib.rs
@@ -1,5 +1,6 @@
 pub mod belief_propagation;
 pub mod lda;
 pub(crate) mod matrixmul;
+pub mod mttest;
 pub mod snr;
 pub mod ttest;


Test patch:
diff --git a/src/scalib/metrics/ttest.py b/src/scalib/metrics/ttest.py
index 1b5c1031..4e9509d6 100644
--- a/src/scalib/metrics/ttest.py
+++ b/src/scalib/metrics/ttest.py
@@ -79,3 +79,82 @@ def fit_u(self, l, x):
     def get_ttest(self):
         r"""Return the current Ttest estimation with an array of shape `(d,ns)`."""
         return self._ttest.get_ttest()
+
+
+class MTtest:
+    r"""The `MTtest` object enables to perform multivariate :math:`t`-test.
+    Concretely, it first computes the  :math:`\mu[j]`'s and :math:`v[j]`'s to
+    use in the :math:`t`-test definition. Especially, :math:`\mu[j]`'s is
+    derived such as:
+
+    For `d=2`:
+
+    .. math::
+        \mu[j] = \frac{1}{n} \sum_{i=0}^{n-1} \prod_{d'=0}^{1}l[i,pois[d',j]] - \bar{l}[:,pois[d',j]]
+
+    For `d>2`:
+
+    .. math::
+        \mu[j] = \frac{1}{n} \sum_{i=0}^{n-1} \prod_{d'=0}^{d-1}
+        \frac{l[i,pois[d',j]] -
+        \bar{l}[:,pois[d',j]]}{\sigma_{l[:,pois[d',j]]}}
+
+    Especially, :math:`\bar{l}` denotes the estimated mean of `l` and
+    :math:`\sigma_l` its standard deviation. `pois` defines the points in the
+    traces for which the (normalized) product will be tested.
+
+    Parameters
+    ----------
+    d : int
+        Maximal statistical order of the :math:`t`-test.
+    pois : array_like, uint32
+        Array of share `(d,n_pois)`. Each column in `pois` will result in a
+        :math:`t`-test of the product of the traces at these $d$ indexes.
+
+    Examples
+    --------
+    >>> from scalib.metrics import MTtest
+    >>> import numpy as np
+    >>> traces = np.random.randint(0,256,(100,200),dtype=np.int16)
+    >>> pois = np.random.randint(0,200,(2,5000),dtype=np.uint32)
+    >>> X = np.random.randint(0,2,100,dtype=np.uint16)
+    >>> mttest = MTtest(d=2,pois=pois)
+    >>> mttest.fit_u(traces,X)
+    >>> t = mttest.get_ttest()
+
+    """
+
+    def __init__(self, d, pois):
+        self._pois = pois
+        self._d = d
+        self._ns = len(pois[0, :])
+
+        self._mttest = _scalib_ext.MTtest(d, self._pois)
+
+    def fit_u(self, l, x):
+        r"""Updates the MTtest estimation with samples of `l` for the sets `x`.
+        This method may be called multiple times.
+
+        Parameters
+        ----------
+        l : array_like, np.int16
+            Array that contains the signal. The array must
+            be of dimension `(n, ns)` and its type must be `np.int16`.
+        x : array_like, np.uint16
+            Set in which each trace belongs. Must be of shape `(n,)`, must be
+            `np.uint16` and must contain only `0` and `1`.
+        """
+        nl, _ = l.shape
+        nx = x.shape[0]
+        if not (nx == nl):
+            raise ValueError(f"Expected x with shape ({nl},)")
+
+        self._mttest.update(l, x)
+
+    def get_ttest(self):
+        r"""Return the current MTtest estimation with an array of shape
+        `(n_pois,)`. Each element in that array corresponds to a test defined by
+        `pois`.
+
+        """
+        return self._mttest.get_ttest()
diff --git a/src/scalib_ext/scalib-py/src/ttest.rs b/src/scalib_ext/scalib-py/src/ttest.rs
index 0f8b122e..5dd21256 100644
--- a/src/scalib_ext/scalib-py/src/ttest.rs
+++ b/src/scalib_ext/scalib-py/src/ttest.rs
@@ -1,6 +1,6 @@
 //! Python wrapper for SCALib's Ttest
 
-use numpy::{PyArray2, PyReadonlyArray1, PyReadonlyArray2, ToPyArray};
+use numpy::{PyArray1, PyArray2, PyReadonlyArray1, PyReadonlyArray2, ToPyArray};
 use pyo3::prelude::*;
 
 #[pyclass]
@@ -35,3 +35,37 @@ impl Ttest {
         Ok(&(ttest.to_pyarray(py)))
     }
 }
+
+#[pyclass]
+pub(crate) struct MTtest {
+    inner: scalib::mttest::MTtest,
+}
+
+#[pymethods]
+impl MTtest {
+    #[new]
+    /// Create a new Ttest state.
+    /// d: order of the Ttest
+    /// pois: points of interest
+    fn new(d: usize, pois: PyReadonlyArray2<u32>) -> Self {
+        let pois = pois.as_array();
+        Self {
+            inner: scalib::mttest::MTtest::new(d, pois.view()),
+        }
+    }
+    /// Update the Ttest state with n fresh traces
+    /// traces: the leakage traces with shape (n,ns)
+    /// y: realization of random variables with shape (n,)
+    fn update(&mut self, py: Python, traces: PyReadonlyArray2<i16>, y: PyReadonlyArray1<u16>) {
+        let traces = traces.as_array();
+        let y = y.as_array();
+        py.allow_threads(|| self.inner.update(traces, y));
+    }
+
+    /// Generate the actual Ttest metric based on the current state.
+    /// return array axes (d,ns)
+    fn get_ttest<'py>(&mut self, py: Python<'py>) -> PyResult<&'py PyArray1<f64>> {
+        let ttest = py.allow_threads(|| self.inner.get_ttest());
+        Ok(&(ttest.to_pyarray(py)))
+    }
+}
diff --git a/src/scalib_ext/scalib/benches/mttest.rs b/src/scalib_ext/scalib/benches/mttest.rs
new file mode 100644
index 00000000..a7e39dcd
--- /dev/null
+++ b/src/scalib_ext/scalib/benches/mttest.rs
@@ -0,0 +1,41 @@
+use criterion::{criterion_group, criterion_main, BenchmarkId, Criterion};
+use ndarray::{Array1, Array2};
+use ndarray_rand::rand_distr::Uniform;
+use ndarray_rand::RandomExt;
+use scalib::mttest;
+use std::time::Duration;
+
+fn bench_mttest(c: &mut Criterion) {
+    let mut group = c.benchmark_group("ttest_update");
+    let d = 2;
+    let traces_len = 1000;
+    let n = 2048;
+    let traces = Array2::<i16>::random((n, traces_len), Uniform::new(0, 1000));
+    let y = Array1::<u16>::random((n,), Uniform::new(0, 2));
+    for d in [2].iter() {
+        for npois in [100000].iter() {
+            let pois = Array2::<u32>::random((*d, *npois), Uniform::new(0, traces_len as u32));
+
+            let mut mtt = mttest::MTtest::new(*d, pois.view());
+            mtt.update(traces.view(), y.view());
+            group.bench_with_input(
+                BenchmarkId::new(format!("mttest_{}", *d), npois),
+                npois,
+                |b, npois| {
+                    b.iter(|| {
+                        mtt.update(traces.view(), y.view());
+                    })
+                },
+            );
+        }
+    }
+    group.finish();
+}
+
+criterion_group! {
+    name = benches;
+    // This can be any expression that returns a `Criterion` object.
+    config = Criterion::default().significance_level(0.01).sample_size(500).measurement_time(Duration::from_secs(60));
+    targets = bench_mttest
+}
+criterion_main!(benches);
diff --git a/src/scalib_ext/scalib/benches/ttest.rs b/src/scalib_ext/scalib/benches/ttest.rs
new file mode 100644
index 00000000..f6523fa2
--- /dev/null
+++ b/src/scalib_ext/scalib/benches/ttest.rs
@@ -0,0 +1,38 @@
+use criterion::{criterion_group, criterion_main, BenchmarkId, Criterion};
+use ndarray::{Array1, Array2};
+use ndarray_rand::rand_distr::Uniform;
+use ndarray_rand::RandomExt;
+use scalib::ttest;
+use std::time::Duration;
+
+fn bench_mttest(c: &mut Criterion) {
+    let mut group = c.benchmark_group("ttest_update");
+    let n = 5000;
+    for d in [2, 3].iter() {
+        for traces_len in [5000, 10000, 20000].iter() {
+            let traces = Array2::<i16>::random((n, *traces_len), Uniform::new(0, 1000));
+            let y = Array1::<u16>::random((n,), Uniform::new(0, 2));
+
+            let mut tt = ttest::Ttest::new(*traces_len, *d as usize);
+            tt.update(traces.view(), y.view());
+            group.bench_with_input(
+                BenchmarkId::new(format!("ttest_{}", traces_len), *d),
+                d,
+                |b, d| {
+                    b.iter(|| {
+                        tt.update(traces.view(), y.view());
+                    })
+                },
+            );
+        }
+    }
+    group.finish();
+}
+
+criterion_group! {
+    name = benches;
+    // This can be any expression that returns a `Criterion` object.
+    config = Criterion::default().significance_level(0.01).sample_size(500).measurement_time(Duration::from_secs(30));
+    targets = bench_mttest
+}
+criterion_main!(benches);
diff --git a/src/scalib_ext/scalib/src/mttest.rs b/src/scalib_ext/scalib/src/mttest.rs
new file mode 100644
index 00000000..3e31e470
--- /dev/null
+++ b/src/scalib_ext/scalib/src/mttest.rs
@@ -0,0 +1,685 @@
+//! Estimation for Multivariate T-test.
+//!
+//! An estimation of MTtest is represented with a MTtest struct. Calling update allows to update
+//! the MTtest state with fresh measurements. get_ttest returns the current value of the estimate.
+//!
+//! MultivarCSAcc computes the central sums CS. MTtest internally uses MultivarCSAcc
+//! to derive statistical moment estimations.
+//!
+//! This is based on the one-pass algorithm proposed in
+//! <https://eprint.iacr.org/2015/207> section 5 as well as <https://doi.org/10.2172/1028931>.
+use itertools::{izip, Itertools};
+use ndarray::{s, Array1, Array2, Array3, ArrayView1, ArrayView2, ArrayView3, Axis};
+use rayon::prelude::*;
+
+// Length of chunck for traces
+const NS_BATCH: usize = 1 << 8;
+
+// Aligned f64 on 256 bits to fit AVX2 instructions.
+#[derive(Clone, Debug, Copy)]
+#[repr(align(32))]
+struct Af64 {
+    x: [f64; 4],
+}
+
+pub struct MultivarCSAcc {
+    /// Number of tuples to evaluate.
+    pub ns: usize,
+    /// pois matrix of size (d,ns)
+    pub pois: Array2<u32>,
+    /// Number of classes
+    pub nc: usize,
+    /// Number of variables in the statistical test (and therefore, its order).
+    pub d: usize,
+    /// Number of samples in each sets. shape (nc,)
+    pub n_traces: Array1<u64>,
+    /// Estimation of all needed centered statistical moments (each multiplied by the number of samples), excluding first-order moments (which are stored in `means`).
+    /// shape of (nc, combis.len(), d)
+    pub cs: Array3<f64>,
+    /// Current means estimates
+    /// shape of (nc, d, ns)
+    pub mean: Array3<f64>,
+    /// List of all the variable combinations to compute the centered moment of.
+    /// This contains all the unique combinations on the set (0..d) U (0..d).
+    /// The combinations are ordered by size.
+    pub combis: Vec<Vec<usize>>,
+    /// Each element in `posi` describes how to compute the corresponding moment in `combis` (`posi.len() == combis.len()`). Each incremental contribution to a moment can be computed as the product between a contribution to a lower-order moment (`posi[i].1`) and a centered data point (`posi[i].0` is the index of the corresponding variable).
+    pub posi: Vec<(i32, i32)>,
+}
+
+impl MultivarCSAcc {
+    /// Creates an MultivarCSAcc.
+    /// It allows to estimate for a class c the central sums (CS) such that
+    /// at index j such that for the set combi:
+    ///
+    ///     cs_Q[c,index(combi),j] = sum_i(prod_{s \in combi} Q[c,pois[s,j]])
+    ///
+    /// for which:
+    ///    - combi: is a unique combination in the multiset (0..d)U(0..d). s are the elements in S.
+    ///    - i: ranges over all the traces in the class c
+    ///    - Q: are the centered traces
+    ///    - pois: see below.
+    ///    - j: Index of the cenral sum to compute.
+    ///
+    /// The update of CS_Q estimation with fresh measurements Q' is done in two steps:
+    ///     1. Compute the CS_Q' of Q' with a two-passes algorithm
+    ///     2. Merge CS_Q and CS_Q' to obtain CS_{Q U Q'}. This is performed
+    ///     thanks to the algorithms in <https://doi.org/10.2172/1028931> and
+    ///     <https://eprint.iacr.org/2015/207>
+    ///
+    /// pois : (d,ns) array where each line corresponds to a poi
+    /// nc : Number of classes to estimate the cs
+    pub fn new(pois: ArrayView2<u32>, nc: usize) -> Self {
+        let ns = pois.shape()[1];
+        let d = pois.shape()[0];
+        let max_set: Vec<usize> = (0..d).chain(0..d).collect();
+
+        // generate all the unique subsets of (0..d, 0..d)
+        // each of them will maintain a state to be updated
+        let combis: Vec<Vec<usize>> = if d > 2 {
+            (1..=(2 * d))
+                .map(|l| {
+                    max_set
+                        .clone()
+                        .into_iter()
+                        .combinations(l)
+                        .map(|mut x| {
+                            x.sort();
+                            x
+                        })
+                        .unique()
+                })
+                .flatten()
+                .collect()
+        } else {
+            // This is hard-coded to ensure that the order matches our hard-coded algorithm for `d=2`.
+            vec![
+                vec![0],
+                vec![1],
+                vec![0, 0],
+                vec![0, 1],
+                vec![1, 1],
+                vec![0, 0, 1],
+                vec![0, 1, 1],
+                vec![0, 0, 1, 1],
+            ]
+        };
+
+        let posi = combis
+            .clone()
+            .iter_mut()
+            .map(|combi| {
+                if combi.len() == 1 {
+                    (combi[0] as i32, -1)
+                } else {
+                    let c = combi.remove(0);
+                    (
+                        c as i32,
+                        combis.iter().position(|x| x == combi).unwrap() as i32,
+                    )
+                }
+            })
+            .collect();
+
+        MultivarCSAcc {
+            ns: ns,
+            pois: pois.to_owned(),
+            nc: nc,
+            d: d,
+            n_traces: Array1::<u64>::zeros((nc,)),
+            cs: Array3::<f64>::zeros((nc, combis.len(), ns)),
+            mean: Array3::<f64>::zeros((nc, d, ns)),
+            combis: combis,
+            posi: posi,
+        }
+    }
+
+    /// Merges the current CS estimate another CS estimate.
+    ///
+    /// cs_q2 : A CS estimate to merge with the current estimation. Its shape is
+    /// (nc,combis.len(),ns)
+    /// u_q2 : mean per class in the other estimation (nc,d,ns)
+    /// n_traces_q2 : Number of traces in each of the classes (nc,)
+
+    // We next describe the merge rule for a single class.
+    //
+    // Definitions:
+    //  - Q1 and Q2 are the set of traces used to build the two CS's to merge.
+    //  - Q = Q1 U Q2
+    //  - n = |Q|
+    //  - 0 <= s < d
+    //  - u[s,j] = mean(Q[:,pois[s,j]])
+    //  - delta[s,j] = u2[s,j] - u1[s,j]
+    //
+    //  Update rule:
+    //  cs_q[combi,j] = cs_q1[combi,j] + cs_q2[combi,j] + (1) + (2)
+    //  with:
+    //
+    //  (1) = sum_{k=2}^{|combi|-1}(
+    //
+    //          sum_{S in |combi| | |S| = k} cs_q1[S,j] prod_{i in (|combi|\S)} (-n_1 / n) *
+    //                                      delta[i,j]
+    //          +
+    //          ((-n2 / n ) ** |combi|) n1 prod_{j in combi) delta[j]
+    //        )
+    //  (2) is the symmetry of (1) by swapping q1 and q2.
+    pub fn merge_from_state(
+        &mut self,
+        cs2: ArrayView3<f64>,
+        u2: ArrayView3<f64>,
+        n2: ArrayView1<u64>,
+    ) {
+        // each row will contain one power of deltas
+        let combis = &self.combis;
+        let ns = self.ns;
+        izip!(
+            self.cs.outer_iter_mut(),
+            self.mean.outer_iter_mut(),
+            self.n_traces.iter(),
+            cs2.outer_iter(),
+            u2.outer_iter(),
+            n2.iter()
+        )
+        .for_each(|(mut cs1, mut u1, n1, cs2, u2, n2)| {
+            let n1 = *n1 as f64;
+            let n2 = *n2 as f64;
+            let n = n1 + n2;
+            if n1 == 0.0 {
+                cs1.assign(&cs2);
+                u1.assign(&u2);
+            } else if n2 != 0.0 {
+                let delta = &u2 - &u1;
+                let mut prod1 = Array1::<f64>::ones((ns,));
+                let mut prod2 = Array1::<f64>::ones((ns,));
+
+                // update all the combinations one by one.
+                for (i, combi) in combis.iter().enumerate().rev() {
+                    // split cs between the inputs and the outputs
+                    let (cs1_smaller, mut cs1_larger) = cs1.view_mut().split_at(Axis(0), i);
+                    let (cs2_smaller, cs2_larger) = cs2.view().split_at(Axis(0), i);
+                    let mut cs1_larger = cs1_larger.slice_mut(s![0, ..]);
+
+                    cs1_larger += &cs2_larger.slice(s![0, ..]);
+
+                    for k in 2..combi.len() {
+                        for set in combi.iter().combinations(k) {
+                            let set: Vec<usize> = set.into_iter().map(|x| *x).collect();
+                            let id = combis.iter().position(|x| *x == set).unwrap();
+
+                            let mut to_multiply: Vec<usize> =
+                                combi.into_iter().map(|x| *x).collect();
+                            for x in set {
+                                to_multiply
+                                    .remove(to_multiply.iter().position(|y| x == *y).unwrap());
+                            }
+
+                            prod1.assign(&cs1_smaller.slice(s![id, ..]));
+                            prod2.assign(&cs2_smaller.slice(s![id, ..]));
+                            for i in to_multiply {
+                                prod1 *= &(&delta.slice(s![i, ..]) * (-n2 / n));
+                                prod2 *= &(&delta.slice(s![i, ..]) * (n1 / n));
+                            }
+                            cs1_larger += &prod1;
+                            cs1_larger += &prod2;
+                        }
+                    }
+
+                    prod1.fill((-n2 / n).powi(combi.len() as i32) * n1);
+                    prod2.fill((n1 / n).powi(combi.len() as i32) * n2);
+                    for x in combi.into_iter() {
+                        prod1 *= &delta.slice(s![*x, ..]);
+                        prod2 *= &delta.slice(s![*x, ..]);
+                    }
+                    cs1_larger += &prod1;
+                    cs1_larger += &prod2;
+                    if combi.len() == 1 {
+                        cs1_larger.fill(0.0);
+                    }
+                }
+                u1 += &((&u2 - &u1) * (n2 / n));
+            }
+        });
+        self.n_traces += &n2;
+    }
+
+    /// Merges to different CS struct estimations.
+    pub fn merge(&mut self, other: &Self) {
+        self.merge_from_state(other.cs.view(), other.mean.view(), other.n_traces.view());
+    }
+
+    /// Updates the current CS estimation with fresh traces and its means per class
+    /// This workds in two steps:
+    ///     1. Computes CS on the fresh traces.
+    ///     2. Merge this with the current CS.
+    ///
+    /// t0 : fresh traces for set 0. of shape (ns, ceil(n0 // 4)).
+    /// t1 : fresh traces for set 1. of shape (ns, ceil(n1 // 4)).
+    /// mean : mean per class of the all traces
+    /// n_traces : count per classes in traces
+    fn update_with_means(
+        &mut self,
+        t0: ArrayView2<Af64>,
+        t1: ArrayView2<Af64>,
+        mean: ArrayView2<f64>,
+        n_traces: ArrayView1<u64>,
+    ) {
+        // Computes CS on the traces with a 2 passes algorithm.
+        let mut cs2 = Array3::<f64>::zeros((2, self.combis.len(), self.pois.shape()[1]));
+        izip!(self.pois.axis_iter(Axis(1)), cs2.axis_iter_mut(Axis(2))).for_each(
+            |(pois, mut cs2)| {
+                if self.d == 2 {
+                    let mut acc00 = Af64 { x: [0.0; 4] };
+                    let mut acc01 = Af64 { x: [0.0; 4] };
+                    let mut acc11 = Af64 { x: [0.0; 4] };
+                    let mut acc001 = Af64 { x: [0.0; 4] };
+                    let mut acc011 = Af64 { x: [0.0; 4] };
+                    let mut acc0011 = Af64 { x: [0.0; 4] };
+                    inner_prod_d2(
+                        &mut acc00,
+                        &mut acc01,
+                        &mut acc11,
+                        &mut acc001,
+                        &mut acc011,
+                        &mut acc0011,
+                        t0.slice(s![pois[0] as usize, ..]).to_slice().unwrap(),
+                        t0.slice(s![pois[1] as usize, ..]).to_slice().unwrap(),
+                    );
+                    for i in 0..4 {
+                        cs2[[0, 2 + 0]] += acc00.x[i];
+                        cs2[[0, 2 + 1]] += acc01.x[i];
+                        cs2[[0, 2 + 2]] += acc11.x[i];
+                        cs2[[0, 2 + 3]] += acc001.x[i];
+                        cs2[[0, 2 + 4]] += acc011.x[i];
+                        cs2[[0, 2 + 5]] += acc0011.x[i];
+                    }
+                } else {
+                    let mut accs = vec![Af64 { x: [0.0; 4] }; self.combis.len()];
+                    let mut prods = vec![Af64 { x: [0.0; 4] }; self.combis.len()];
+                    let ts = pois
+                        .iter()
+                        .map(|x| t0.slice(s![*x as usize, ..]).to_slice().unwrap())
+                        .collect();
+
+                    inner_prod_generic(
+                        accs.as_mut_slice(),
+                        prods.as_mut_slice(),
+                        self.posi.as_slice(),
+                        &ts,
+                    );
+                    for i in 0..4 {
+                        for j in 0..self.combis.len() {
+                            cs2[[0, j]] += accs[j].x[i];
+                        }
+                    }
+                }
+
+                if self.d == 2 {
+                    let mut acc00 = Af64 { x: [0.0; 4] };
+                    let mut acc01 = Af64 { x: [0.0; 4] };
+                    let mut acc11 = Af64 { x: [0.0; 4] };
+                    let mut acc001 = Af64 { x: [0.0; 4] };
+                    let mut acc011 = Af64 { x: [0.0; 4] };
+                    let mut acc0011 = Af64 { x: [0.0; 4] };
+                    inner_prod_d2(
+                        &mut acc00,
+                        &mut acc01,
+                        &mut acc11,
+                        &mut acc001,
+                        &mut acc011,
+                        &mut acc0011,
+                        t1.slice(s![pois[0] as usize, ..]).to_slice().unwrap(),
+                        t1.slice(s![pois[1] as usize, ..]).to_slice().unwrap(),
+                    );
+                    for i in 0..4 {
+                        cs2[[1, 2 + 0]] += acc00.x[i];
+                        cs2[[1, 2 + 1]] += acc01.x[i];
+                        cs2[[1, 2 + 2]] += acc11.x[i];
+                        cs2[[1, 2 + 3]] += acc001.x[i];
+                        cs2[[1, 2 + 4]] += acc011.x[i];
+                        cs2[[1, 2 + 5]] += acc0011.x[i];
+                    }
+                } else {
+                    let mut prods = vec![Af64 { x: [0.0; 4] }; self.combis.len()];
+                    let mut accs = vec![Af64 { x: [0.0; 4] }; self.combis.len()];
+                    let ts = pois
+                        .iter()
+                        .map(|x| t1.slice(s![*x as usize, ..]).to_slice().unwrap())
+                        .collect();
+
+                    inner_prod_generic(
+                        accs.as_mut_slice(),
+                        prods.as_mut_slice(),
+                        self.posi.as_slice(),
+                        &ts,
+                    );
+                    for i in 0..4 {
+                        for j in 0..self.combis.len() {
+                            cs2[[1, j]] += accs[j].x[i];
+                        }
+                    }
+                }
+            },
+        );
+
+        // Genereted the means according to the pois
+        let pois = &self.pois;
+        let mut mapped_means = Array3::<f64>::zeros((self.nc, self.d, self.ns));
+        izip!(mapped_means.outer_iter_mut(), mean.outer_iter()).for_each(
+            |(mut to_update, mean)| {
+                izip!(to_update.outer_iter_mut(), pois.outer_iter()).for_each(
+                    |(mut to_update, pois)| {
+                        to_update.zip_mut_with(&pois, |x, p| *x = mean[*p as usize]);
+                    },
+                );
+            },
+        );
+
+        // Merge CS of the inputs with self.
+        self.merge_from_state(cs2.view(), mapped_means.view(), n_traces.view());
+    }
+
+    /// Updates the current CS estimation with fresh traces
+    /// traces : fresh traces
+    /// y : class corresponding to each traces
+    pub fn update(&mut self, traces: ArrayView2<i16>, y: ArrayView1<u16>) {
+        let (mean, n_traces) = means_per_class(traces, y, self.nc);
+        let (t0, t1) = center_transpose_align(traces, mean.view(), y);
+        self.update_with_means(t0.view(), t1.view(), mean.view(), n_traces.view());
+    }
+}
+
+pub struct MTtest {
+    /// order of the test
+    d: usize,
+    /// Number of samples per trace
+    ns: usize,
+    /// Vector of CS accumulators. This will be used to
+    /// derive the statistical moments.
+    accumulators: Vec<MultivarCSAcc>,
+    /// Pois to combine in the multivariate T-test (d,ns)
+    pois: Array2<u32>,
+}
+
+impl MTtest {
+    /// Create a new Ttest state.
+    /// ns: traces length
+    /// d: order of the Ttest
+    pub fn new(d: usize, pois: ArrayView2<u32>) -> Self {
+        assert!(d == pois.shape()[0]);
+        assert!(
+            d > 1,
+            "Order of Multivariate T-test should be larger than 1, provided d = {}",
+            d
+        );
+
+        // generates the MultivarCSAcc underlying MTtest.
+        let accumulators = build_accumulator(pois.view());
+        MTtest {
+            d: d,
+            ns: pois.shape()[1],
+            accumulators: accumulators,
+            pois: pois.to_owned(),
+        }
+    }
+
+    /// Update the MTtest state with n fresh traces
+    /// traces: the leakage traces with shape (n,ns)
+    /// y: realization of random variables with shape (n,)
+    pub fn update(&mut self, traces: ArrayView2<i16>, y: ArrayView1<u16>) {
+        // First pass to compute the means, center and align.
+        let (mean, n_traces) = means_per_class(traces, y, 2);
+        let (t0, t1) = center_transpose_align(traces, mean.view(), y);
+
+        // chunck different traces for more threads
+        (
+            self.pois.axis_chunks_iter(Axis(1), NS_BATCH),
+            &mut self.accumulators,
+        )
+            .into_par_iter()
+            .for_each(|(_, acc)| {
+                // chunck the traces with their lenght
+                acc.update_with_means(t0.view(), t1.view(), mean.view(), n_traces.view());
+            });
+    }
+
+    /// Computes the t statistic according to the
+    /// order d.
+    pub fn get_ttest(&self) -> Array1<f64> {
+        let mut t = Array1::<f64>::zeros((self.ns,));
+        izip!(
+            t.axis_chunks_iter_mut(Axis(0), NS_BATCH),
+            self.accumulators.iter()
+        )
+        .for_each(|(mut t, acc)| {
+            t.fill(1.0);
+            let n1 = acc.n_traces[0] as f64;
+            let n2 = acc.n_traces[1] as f64;
+
+            let mut mu1 = Array1::<f64>::zeros((t.shape()[0],));
+            let mut var1 = Array1::<f64>::zeros((t.shape()[0],));
+            let mut mu2 = Array1::<f64>::zeros((t.shape()[0],));
+            let mut var2 = Array1::<f64>::zeros((t.shape()[0],));
+
+            // assign means
+            let combi: Vec<usize> = (0..self.d).collect();
+            let id = acc.combis.iter().position(|x| *x == combi).unwrap();
+            let mus = acc.cs.slice(s![.., id, ..]);
+            mu1.assign(&(&mus.slice(s![0, ..]) / n1));
+            mu2.assign(&(&mus.slice(s![1, ..]) / n2));
+
+            let mut combi: Vec<usize> = (0..self.d).chain(0..self.d).collect();
+            combi.sort();
+            let id = acc.combis.iter().position(|x| *x == combi).unwrap();
+            let vars = acc.cs.slice(s![.., id, ..]);
+            var1.assign(&(&vars.slice(s![0, ..]) / n1));
+            var2.assign(&(&vars.slice(s![1, ..]) / n2));
+
+            if self.d > 2 {
+                for j in 0..self.d {
+                    let combi: Vec<usize> = vec![j, j];
+                    let id = acc.combis.iter().position(|x| *x == combi).unwrap();
+                    let mus = acc.cs.slice(s![.., id, ..]);
+                    mu1 /= &(mus.slice(s![0, ..]).mapv(|x| (x / n1).sqrt()));
+                    mu2 /= &(mus.slice(s![1, ..]).mapv(|x| (x / n2).sqrt()));
+
+                    var1 /= &(mus.slice(s![0, ..]).mapv(|x| (x / n1)));
+                    var2 /= &(mus.slice(s![1, ..]).mapv(|x| (x / n2)));
+                }
+            }
+            var1 -= &(&mu1 * &mu1);
+            var2 -= &(&mu2 * &mu2);
+
+            t.assign(&(&mu1 - &mu2));
+
+            t /= &((&var1 / n1) + (&var2 / n2)).mapv(|x| x.sqrt());
+        });
+        t
+    }
+}
+
+/// Computes the means per class
+fn means_per_class(
+    traces: ArrayView2<i16>,
+    y: ArrayView1<u16>,
+    nc: usize,
+) -> (Array2<f64>, Array1<u64>) {
+    let mut sum = Array2::<i32>::zeros((nc, traces.shape()[1]));
+    let mut sum64 = Array2::<i64>::zeros((nc, traces.shape()[1]));
+    let mut n_traces = Array1::<u64>::zeros(nc);
+
+    for (trace, class) in traces.outer_iter().zip(y.iter()) {
+        n_traces[*class as usize] += 1;
+        let mut s = sum.slice_mut(s![*class as usize, ..]);
+        if (n_traces[*class as usize] % (1 << 16)) == ((1 << 16) - 1) {
+            let mut s64 = sum64.slice_mut(s![*class as usize, ..]);
+            s64 += &s.mapv(|x| x as i64);
+            s.fill(0);
+        }
+        s.zip_mut_with(&trace, |s, t| {
+            *s += *t as i32;
+        });
+    }
+    sum64 += &sum.mapv(|x| x as i64);
+    let n = n_traces.mapv(|x| x as f64);
+    let mut mean = sum64.mapv(|x| x as f64);
+    mean.axis_iter_mut(Axis(1)).for_each(|mut m| m /= &n);
+
+    (mean, n_traces)
+}
+
+/// Generates accumulators according to the pois.
+fn build_accumulator(pois: ArrayView2<u32>) -> Vec<MultivarCSAcc> {
+    // number of required accumulators
+    let ns = pois.shape()[1];
+    let n_batches = ((ns as f64) / (NS_BATCH as f64)).ceil() as usize;
+    let accumulators: Vec<MultivarCSAcc> = (0..n_batches)
+        .map(|x| {
+            let l = std::cmp::min(ns - (x * NS_BATCH), NS_BATCH);
+            MultivarCSAcc::new(pois.slice(s![.., (NS_BATCH * x)..(NS_BATCH * x + l)]), 2)
+        })
+        .collect();
+
+    accumulators
+}
+
+#[inline(never)]
+fn inner_prod_generic(
+    accs: &mut [Af64],
+    prods: &mut [Af64],
+    posi: &[(i32, i32)],
+    ts: &Vec<&[Af64]>,
+) {
+    for j in 0..ts[0].len() {
+        for i in 0..prods.len() {
+            let (t, c) = posi[i];
+            let (prod_old, mut prod_up) = prods.split_at_mut(i);
+            let (_, mut accs_up) = accs.split_at_mut(i);
+
+            let t = &ts[t as usize][j];
+            if c == -1 {
+                prod_up[0].x[0] = t.x[0];
+                prod_up[0].x[1] = t.x[1];
+                prod_up[0].x[2] = t.x[2];
+                prod_up[0].x[3] = t.x[3];
+            } else {
+                prod_up[0].x[0] = prod_old[c as usize].x[0] * t.x[0];
+                prod_up[0].x[1] = prod_old[c as usize].x[1] * t.x[1];
+                prod_up[0].x[2] = prod_old[c as usize].x[2] * t.x[2];
+                prod_up[0].x[3] = prod_old[c as usize].x[3] * t.x[3];
+            }
+
+            accs_up[0].x[0] += prod_up[0].x[0];
+            accs_up[0].x[1] += prod_up[0].x[1];
+            accs_up[0].x[2] += prod_up[0].x[2];
+            accs_up[0].x[3] += prod_up[0].x[3];
+        }
+    }
+}
+
+#[inline(never)]
+fn inner_prod_d2(
+    acc00: &mut Af64,
+    acc01: &mut Af64,
+    acc11: &mut Af64,
+    acc001: &mut Af64,
+    acc011: &mut Af64,
+    acc0011: &mut Af64,
+
+    t0: &[Af64],
+    t1: &[Af64],
+) {
+    izip!(t0.iter(), t1.iter()).for_each(|(t0, t1)| {
+        acc00.x[0] += t0.x[0] * t0.x[0];
+        acc00.x[1] += t0.x[1] * t0.x[1];
+        acc00.x[2] += t0.x[2] * t0.x[2];
+        acc00.x[3] += t0.x[3] * t0.x[3];
+
+        acc01.x[0] += t0.x[0] * t1.x[0];
+        acc01.x[1] += t0.x[1] * t1.x[1];
+        acc01.x[2] += t0.x[2] * t1.x[2];
+        acc01.x[3] += t0.x[3] * t1.x[3];
+
+        acc11.x[0] += t1.x[0] * t1.x[0];
+        acc11.x[1] += t1.x[1] * t1.x[1];
+        acc11.x[2] += t1.x[2] * t1.x[2];
+        acc11.x[3] += t1.x[3] * t1.x[3];
+
+        acc001.x[0] += t0.x[0] * t0.x[0] * t1.x[0];
+        acc001.x[1] += t0.x[1] * t0.x[1] * t1.x[1];
+        acc001.x[2] += t0.x[2] * t0.x[2] * t1.x[2];
+        acc001.x[3] += t0.x[3] * t0.x[3] * t1.x[3];
+
+        acc011.x[0] += t0.x[0] * t1.x[0] * t1.x[0];
+        acc011.x[1] += t0.x[1] * t1.x[1] * t1.x[1];
+        acc011.x[2] += t0.x[2] * t1.x[2] * t1.x[2];
+        acc011.x[3] += t0.x[3] * t1.x[3] * t1.x[3];
+
+        acc0011.x[0] += t0.x[0] * t0.x[0] * t1.x[0] * t1.x[0];
+        acc0011.x[1] += t0.x[1] * t0.x[1] * t1.x[1] * t1.x[1];
+        acc0011.x[2] += t0.x[2] * t0.x[2] * t1.x[2] * t1.x[2];
+        acc0011.x[3] += t0.x[3] * t0.x[3] * t1.x[3] * t1.x[3];
+    });
+}
+
+fn center_transpose_align(
+    traces: ArrayView2<i16>,
+    means: ArrayView2<f64>,
+    y: ArrayView1<u16>,
+) -> (Array2<Af64>, Array2<Af64>) {
+    let ns = traces.shape()[1];
+
+    let posi0: Vec<usize> = y
+        .iter()
+        .enumerate()
+        .filter(|(_, x)| **x == 0)
+        .map(|(i, _)| i as usize)
+        .collect();
+    let n0 = posi0.len();
+
+    let mut t0 = Array2::<Af64>::from_elem(
+        (ns, (n0 as f64 / 4.0).ceil() as usize),
+        Af64 {
+            x: [0.0, 0.0, 0.0, 0.0],
+        },
+    );
+    let posi1: Vec<usize> = y
+        .iter()
+        .enumerate()
+        .filter(|(_, x)| **x == 1)
+        .map(|(i, _)| i as usize)
+        .collect();
+    let n1 = posi1.len();
+
+    let mut t1 = Array2::<Af64>::from_elem(
+        (ns, (n1 as f64 / 4.0).ceil() as usize),
+        Af64 {
+            x: [0.0, 0.0, 0.0, 0.0],
+        },
+    );
+
+    (traces.axis_iter(Axis(1)), t0.axis_iter_mut(Axis(0)))
+        .into_par_iter()
+        .enumerate()
+        .for_each(|(i, (traces, mut t0))| {
+            let m0 = means[[0, i]];
+            izip!(&posi0.iter().chunks(4), t0.iter_mut()).for_each(|(p, t)| {
+                p.enumerate().for_each(|(x, p)| {
+                    t.x[x] = traces[*p] as f64 - m0;
+                });
+            });
+        });
+
+    (traces.axis_iter(Axis(1)), t1.axis_iter_mut(Axis(0)))
+        .into_par_iter()
+        .enumerate()
+        .for_each(|(i, (traces, mut t1))| {
+            let m1 = means[[1, i]];
+            izip!(&posi1.iter().chunks(4), t1.iter_mut()).for_each(|(p, t)| {
+                p.enumerate().for_each(|(x, p)| {
+                    t.x[x] = traces[*p] as f64 - m1;
+                });
+            });
+        });
+
+    (t0, t1)
+}
diff --git a/src/scalib_ext/scalib/src/ttest.rs b/src/scalib_ext/scalib/src/ttest.rs
index 8975d356..466c69fb 100644
--- a/src/scalib_ext/scalib/src/ttest.rs
+++ b/src/scalib_ext/scalib/src/ttest.rs
@@ -1,4 +1,4 @@
-//! Estimation for higher-order T-test.
+//! Estimation for higher-order Univariate T-test.
 //!
 //! An estimation of Ttest is represented with a Ttest struct. Calling update allows
 //! to update the Ttest state with fresh measurements. get_ttest returns the current value
@@ -7,168 +7,286 @@
 //!
 //! This is based on the one-pass algorithm proposed in
 //! <https://eprint.iacr.org/2015/207>.
-
-use ndarray::{s, Array1, Array2, Array3, ArrayView1, ArrayView2, Axis};
+use itertools::izip;
+use ndarray::{s, Array1, Array2, Array3, ArrayView1, ArrayView2, ArrayView3, ArrayViewMut1, Axis};
 use num_integer::binomial;
 use rayon::prelude::*;
+use std::cmp;
+const NS_BATCH: usize = 1 << 12;
+const Y_BATCH: usize = 1 << 9;
+
+pub struct UniCSAcc {
+    /// Number of samples in trace
+    pub ns: usize,
+    /// pub Number of classes
+    pub nc: usize,
+    /// Highest moment to estimate
+    pub d: usize,
+    /// Number of samples in sets
+    pub n_traces: Array1<u64>,
+    /// Array containing the current estimation
+    pub moments: Array3<f64>,
+}
+
+impl UniCSAcc {
+    /// Creates an UniCSAcc
+    /// ns : number of point in traces
+    /// d : higher power to estimate
+    /// nc: number of classes to estimate the CS for.
+    pub fn new(ns: usize, d: usize, nc: usize) -> Self {
+        UniCSAcc {
+            ns: ns,
+            d: d,
+            nc: nc,
+            n_traces: Array1::<u64>::zeros(nc),
+            moments: Array3::<f64>::zeros((nc, d, ns)),
+        }
+    }
+    /// Merges the current accumulator with anoter sum of
+    /// centered products
+    /// moments_other : (nc,d,ns) matrix
+    /// n_traces : (nc,) matrix
+    // See
+    // [1]: https://www.osti.gov/biblio/1028931-formulas-robust-one-pass-parallel-computation-covariances-arbitrary-order-statistical-moments
+    // [2]: https://eprint.iacr.org/2015/207.pdf
+    // for additional details about the merging algorithm.
+    //
+    // Q0 set of all the traces in current estimation
+    // Q1 set of all the traces in the estiamation to merge
+    // Q = Q0 U Q1 is the merge of all traces
+    //
+    // Definitions:
+    //      M_d,Qi = sum_{Qi} (l - u_i) ^ d
+    //      u_Qi = 1/(|Qi|) * sum_{Qi} l
+    //      delta_1,0 = u_0 - u_1
+    //      ni = |Qi|
+    //      n = n0 + n1
+    //
+    // Update equation (Eq. 2.1 in [1]):
+    //
+    // M_Q_d = M_Q0_d + M_Q1_d +
+    //      sum_{k=1}_{p-2}(
+    //           binomial(k,p) *
+    //
+    //              delta_1,0 ^ d *
+    //              (
+    //                  ((-n1/n) ^ k * M_p-k,0) +
+    //                   (n0/n) ^ k * M_p-k,1))
+    //              )
+    //      )
+    //      + ((n0 * n1) / n ) *( delta_1,0 ^ p)
+    //              * ( 1 / (n1 ^ (p-1)) - (-1 / n0)^(p-1))
+    pub fn merge_from_state(&mut self, moments_other: ArrayView3<f64>, n_traces: ArrayView1<u64>) {
+        // compute the deltas u1 - u0 for all the classes
+        let delta = &moments_other.slice(s![.., 0, ..]) - &self.moments.slice(s![.., 0, ..]);
+        let d = self.d;
+        let ns = self.ns;
+
+        // each row will contain one power of deltas
+        let mut delta_pows = Array2::<f64>::zeros((d + 1, ns));
+
+        // Update all the classes estimated independently
+        izip!(
+            &mut self.moments.outer_iter_mut(),
+            self.n_traces.iter(),
+            moments_other.outer_iter(),
+            n_traces.iter(),
+            delta.outer_iter()
+        )
+        .for_each(|(mut cs0, n0, cs1, n1, delta)| {
+            let n = (*n0 + *n1) as f64;
+            let n0 = *n0 as f64;
+            let n1 = *n1 as f64;
+
+            if n0 == 0.0 {
+                // if current estimate has no sample, simply assign from the other estimate.
+                cs0.assign(&cs1);
+            } else if n1 != 0.0 {
+                // compute powers of deltas
+                delta_pows
+                    .axis_iter_mut(Axis(0))
+                    .enumerate()
+                    .for_each(|(i, mut x)| x.assign(&delta.mapv(|t| t.powi(i as i32))));
+
+                // first update the higher powers since their update rules
+                // is based on the one of smaller powers.
+                for p in (2..d + 1).rev() {
+                    let p = p as i32;
 
+                    let (as_input0, mut to_update0) =
+                        cs0.view_mut().split_at(Axis(0), (p - 1) as usize);
+                    let (as_input1, to_update1) = cs1.view().split_at(Axis(0), (p - 1) as usize);
+
+                    let mut to_update0 = to_update0.slice_mut(s![0, ..]);
+                    to_update0 += &to_update1.slice(s![0, ..]);
+
+                    for k in 1..(p - 1) {
+                        let cst = binomial(p, k) as f64;
+                        let mut tmp = &delta_pows.slice(s![k, ..]) * cst;
+                        let tmp2 = &as_input0.slice(s![p - k - 1, ..]) * ((-n1 / n).powi(k));
+                        let tmp3 = &as_input1.slice(s![p - k - 1, ..]) * ((n0 / n).powi(k));
+                        let x = tmp2 + tmp3;
+                        tmp = &tmp * x;
+                        to_update0 += &tmp;
+                    }
+
+                    let mut cst = (1.0 / n1).powi(p - 1) - (-1.0 / n0).powi(p - 1);
+                    cst *= (n1 * n0 / n).powi(p);
+                    to_update0 += &(&delta_pows.slice(s![p, ..]) * cst);
+                }
+
+                // update the mean
+                let mut u0 = cs0.slice_mut(s![0, ..]);
+                u0 += &(&delta * (n1 / n));
+            }
+        });
+        self.n_traces += &n_traces;
+    }
+
+    /// Merges to different CS estimations.
+    pub fn merge(&mut self, other: &Self) {
+        self.merge_from_state(other.moments.view(), other.n_traces.view());
+    }
+
+    /// Updates the current estimation with fresh traces.
+    pub fn update(&mut self, traces: ArrayView2<i16>, y: ArrayView1<u16>) {
+        let mut sum = Array2::<i64>::zeros((self.nc, self.ns));
+        let mut moments_other = Array3::<f64>::zeros((self.nc, self.d, self.ns));
+        let mut n_traces = Array1::<u64>::zeros(self.nc);
+
+        // STEP 1: 2-passes algorithm to compute center sum of powers
+        // sum per class,
+        for (trace, class) in traces.outer_iter().zip(y.iter()) {
+            n_traces[*class as usize] += 1;
+            let mut s = sum.slice_mut(s![*class as usize, ..]);
+            acc_sum(&mut s.view_mut(), &trace);
+        }
+
+        // assign mean
+        let n = n_traces.mapv(|x| x as f64);
+        let mut mean = sum.mapv(|x| x as f64);
+        mean.axis_iter_mut(Axis(1)).for_each(|mut m| m /= &n);
+        moments_other.slice_mut(s![.., 0, ..]).assign(&mean);
+
+        // compute centered sums
+        let mut pow = Array1::<f64>::zeros((self.ns,));
+        let mut t = Array1::<f64>::zeros((self.ns,));
+        for (trace, class) in traces.outer_iter().zip(y.iter()) {
+            let mut m_full = moments_other.slice_mut(s![*class as usize, .., ..]);
+            izip!(
+                t.view_mut().into_slice().unwrap().iter_mut(),
+                pow.view_mut().into_slice().unwrap().iter_mut(),
+                m_full.slice(s![0, ..]).view().to_slice().unwrap().iter(),
+                trace.view().to_slice().unwrap().iter()
+            )
+            .for_each(|(t, pow, m_full, trace)| {
+                *t = *trace as f64 - m_full;
+                *pow = *t;
+            });
+            for d in 2..(self.d + 1) {
+                let dest = m_full.slice_mut(s![d - 1, ..]);
+                pow_and_sum(dest, pow.view_mut(), t.view());
+            }
+        }
+
+        // STEP 2 merge the two pass estimation with the current state of the accumulator
+        self.merge_from_state(moments_other.view(), n_traces.view());
+    }
+
+    /// Reset the current accumulator
+    pub fn reset(&mut self) {
+        self.n_traces.fill(0);
+        self.moments.fill(0.0);
+    }
+}
 pub struct Ttest {
-    /// Central sums of order 1 up to order d*2 with shape (ns,2,2*d),
-    /// where central sums is sum((x-u_x)**i).
-    /// Axes are (class, trace sample, order).
-    /// cs[..,..,0] contains the current estimation of means instead of
-    /// the central sum (which would be zero).
-    /// number of samples in a trace
-    cs: Array3<f64>,
-    /// number of samples per class (2,)
-    n_samples: Array1<u64>,
     /// order of the test
     d: usize,
     /// Number of samples per trace
     ns: usize,
+    /// Vector of Moment accumulators
+    accumulators: Vec<UniCSAcc>,
 }
 
+pub fn build_accumulator(ns: usize, d: usize) -> Vec<UniCSAcc> {
+    let n_batches = ((ns as f64) / (NS_BATCH as f64)).ceil() as usize;
+    let accumulators: Vec<UniCSAcc> = (0..n_batches)
+        .map(|x| {
+            let l = std::cmp::min(ns - (x * NS_BATCH), NS_BATCH);
+            UniCSAcc::new(l, 2 * d, 2)
+        })
+        .collect();
+    accumulators
+}
 impl Ttest {
     /// Create a new Ttest state.
     /// ns: traces length
     /// d: order of the Ttest
     pub fn new(ns: usize, d: usize) -> Self {
+        // number of required accumulators
+        let accumulators = build_accumulator(ns, d);
         Ttest {
-            cs: Array3::<f64>::zeros((ns, 2, 2 * d)),
-            n_samples: Array1::<u64>::zeros((2,)),
             d: d,
             ns: ns,
+            accumulators: accumulators,
         }
     }
     /// Update the Ttest state with n fresh traces
     /// traces: the leakage traces with shape (n,ns)
     /// y: realization of random variables with shape (n,)
-    // Q set of all previous traces
-    //
-    // Initial values, |Q| = n-1
-    // mu = (1/(n-1)) * sum(Q)
-    // if d>=2:
-    //      CS_{d,Q} = 1/(n-1) sum((Q-mu)**d)
-    // else:
-    //      CS_{1,Q'} = mu'
-    //
-    // Updated with a single measurement t, Q' = Q U t, |Q'| = n
-    // mu' = (1/(n-1)) * sum(Q')
-    //
-    // if d >=2:
-    //      CS_{d,Q'} = (1/n) sum((Q-mu')**d)
-    // else:
-    //      CS_{1,Q'} = mu'
-    //
-    // Update rules is given by:
-    //
-    // delta = (t - mu)/n
-    //
-    // CS_{d,Q'} = CS_{d,Q}
-    //      + sum_{k=1,d-2}(
-    //         binomial(k,d)
-    //         * CS_{d-k,Q}
-    //         * (-delta)**d
-    //         )
-    //      + (
-    //          (delta * (n-1))**d
-    //          *(1 - (-1/(n-1)))**(d-1)
-    //        )
-    //
-    // mu' = mu+delta
     pub fn update(&mut self, traces: ArrayView2<i16>, y: ArrayView1<u16>) {
         let d = self.d;
+        let ns = self.ns;
+        let n_traces = traces.shape()[0];
+        let ns_chuncks = cmp::max(1, ns / NS_BATCH);
+        let min_desired_chuncks = 4 * rayon::current_num_threads();
 
-        // pre computes the combinatorial factors
-        let cbs: Vec<(usize, Vec<(f64, usize)>)> = (2..((2 * self.d) + 1))
-            .rev()
-            .map(|j| {
-                (
-                    j,
-                    (1..(j - 1)).map(|k| (binomial(j, k) as f64, k)).collect(),
-                )
-            })
-            .collect();
-
-        // contains the data that are the same for all the points in a single traces
-        // Contains tupes (n, y, mults): (f64, usize, Vec<f64>)
-        // n : number of previously processed traces for the class y
-        // y : set to update
-        // mults: (n-1)**(j) * (1.0 - (-1.0/(n-1.0))**(j-1) for j in [2..(2*d-1)].rev()
-        let shared_data: Vec<(f64, usize, Vec<f64>)> = y
-            .iter()
-            .map(|y| {
-                let y = *y as usize;
-                assert!(y <= 1);
-
-                // update the number of observations
-                let n = &mut self.n_samples[y];
-                *n += 1;
-                let n = *n as f64;
+        let y_chunck_size = if min_desired_chuncks < ns_chuncks {
+            n_traces
+        } else {
+            let tmp = cmp::min(
+                rayon::current_num_threads(),
+                min_desired_chuncks / ns_chuncks,
+            ); // ensure that we do not split in more than available threads.
+            cmp::max(256, n_traces / tmp)
+        };
 
+        let res = (
+            traces.axis_chunks_iter(Axis(0), y_chunck_size),
+            y.axis_chunks_iter(Axis(0), y_chunck_size),
+        )
+            .into_par_iter()
+            .map(|(traces, y)| {
+                // chunck different traces for more threads
+                let mut accumulators = build_accumulator(ns, d);
                 (
-                    // number of sample on that class
-                    n,
-                    // y value
-                    y,
-                    // compute the multiplicative factor similar for all trace samples
-                    cbs.iter()
-                        .map(|(j, _)| {
-                            (n - 1.0).powi(*j as i32)
-                                * (1.0 - (-1.0 / (n - 1.0)).powi(*j as i32 - 1))
-                        })
-                        .collect(),
+                    traces.axis_chunks_iter(Axis(1), NS_BATCH),
+                    &mut accumulators,
                 )
+                    .into_par_iter()
+                    .for_each(|(traces, acc)| {
+                        // chunck the traces with their lenght
+                        izip!(
+                            traces.axis_chunks_iter(Axis(0), Y_BATCH),
+                            y.axis_chunks_iter(Axis(0), Y_BATCH)
+                        )
+                        .for_each(|(traces, y)| acc.update(traces, y));
+                    });
+                accumulators
             })
-            .collect();
-
-        (traces.axis_iter(Axis(1)), self.cs.axis_iter_mut(Axis(0)))
-            .into_par_iter()
-            .for_each_init(
-                || Array1::<f64>::zeros(2 * d),
-                |ref mut delta_pows, (traces, mut cs)| {
-                    traces
-                        .iter()
-                        .zip(shared_data.iter())
-                        .for_each(|(trace, (n, y, mults))| {
-                            let mut cs_s = cs.slice_mut(s![*y, ..]);
-                            let cs = cs_s.as_slice_mut().unwrap();
-
-                            // compute the delta
-                            let delta = ((*trace as f64) - cs[0]) / (*n as f64);
-
-                            // delta_pows[i] = delta ** (i+1)
-                            // We will need all of them next
-                            delta_pows.iter_mut().fold(delta, |acc, x| {
-                                *x = acc;
-                                acc * delta
-                            });
-
-                            // apply the one-pass update rule
-                            cbs.iter().zip(mults.iter()).for_each(|((j, vec), mult)| {
-                                if *n > 1.0 {
-                                    cs[*j - 1] += delta_pows[*j - 1] * mult;
-                                }
-                                vec.iter().for_each(|(cb, k)| {
-                                    let a = cs[*j - *k - 1];
-                                    if (k & 0x1) == 1 {
-                                        // k is not pair
-                                        cs[*j - 1] -= cb * delta_pows[*k - 1] * a;
-                                    } else {
-                                        // k is pair
-                                        cs[*j - 1] += cb * delta_pows[*k - 1] * a;
-                                    }
-                                });
-                            });
-                            cs[0] += delta;
-                        });
+            .reduce(
+                || build_accumulator(ns, d),
+                |mut x, y| {
+                    // accumulate all to the self accumulator
+                    x.iter_mut().zip(y.iter()).for_each(|(x, y)| x.merge(y));
+                    x
                 },
             );
+        izip!(self.accumulators.iter_mut(), res.iter()).for_each(|(x, y)| x.merge(y));
     }
 
     /// Generate the actual Ttest metric based on the current state.
     /// return array axes (d,ns)
-    //
     // with central moment defined as:
     //   CM_{i,Q} = CS_{i,Q}/n
     //
@@ -184,59 +302,122 @@ impl Ttest {
     // d > 2:
     //      ui = CM_{d,Q} / CM_{2,Q}**(d/2)
     //      vi = (CM_{2*d,Q} - CM_{d,Q}**2) / CM{2,Q}**d
-
     pub fn get_ttest(&self) -> Array2<f64> {
         let mut ttest = Array2::<f64>::zeros((self.d, self.ns));
-        let cs = &self.cs;
-        let n_samples = &self.n_samples;
+        let n_samples = &self.accumulators[0].n_traces;
 
         let n0 = n_samples[[0]] as f64;
         let n1 = n_samples[[1]] as f64;
 
-        (
-            ttest.axis_chunks_iter_mut(Axis(1), 20),
-            cs.axis_chunks_iter(Axis(0), 20),
+        izip!(
+            ttest.axis_chunks_iter_mut(Axis(1), NS_BATCH),
+            self.accumulators.iter()
         )
-            .into_par_iter()
-            .for_each(|(mut ttest, cs)| {
-                ttest
-                    .axis_iter_mut(Axis(1))
-                    .zip(cs.axis_iter(Axis(0)))
-                    .for_each(|(mut ttest, cs)| {
-                        let mut u0;
-                        let mut u1;
-                        let mut v0;
-                        let mut v1;
-                        for d in 1..(self.d + 1) {
-                            if d == 1 {
-                                u0 = cs[[0, 0]];
-                                u1 = cs[[1, 0]];
-
-                                v0 = cs[[0, 1]] / n0;
-                                v1 = cs[[1, 1]] / n1;
-                            } else if d == 2 {
-                                u0 = cs[[0, 1]] / n0;
-                                u1 = cs[[1, 1]] / n1;
-
-                                v0 = cs[[0, 3]] / n0 - ((cs[[0, 1]] / n0).powi(2));
-                                v1 = cs[[1, 3]] / n1 - ((cs[[1, 1]] / n1).powi(2));
-                            } else {
-                                u0 = (cs[[0, d - 1]] / n0)
-                                    / ((cs[[0, 1]] / n0).powf(d as f64 / 2.0));
-                                u1 = (cs[[1, d - 1]] / n1)
-                                    / ((cs[[1, 1]] / n1).powf(d as f64 / 2.0));
-
-                                v0 = cs[[0, (2 * d) - 1]] / n0 - ((cs[[0, d - 1]] / n0).powi(2));
-                                v0 /= (cs[[0, 1]] / n0).powi(d as i32);
-
-                                v1 = cs[[1, (2 * d) - 1]] / n1 - ((cs[[1, d - 1]] / n1).powi(2));
-                                v1 /= (cs[[1, 1]] / n1).powi(d as i32);
-                            }
-
-                            ttest[d - 1] = (u0 - u1) / f64::sqrt((v0 / n0) + (v1 / n1));
-                        }
-                    });
-            });
+        .for_each(|(mut ttest, acc)| {
+            ttest
+                .axis_iter_mut(Axis(0))
+                .enumerate()
+                .for_each(|(d, mut ttest)| {
+                    let d = d + 1;
+
+                    if d == 1 {
+                        let u0 = acc.moments.slice(s![0, 0, ..]);
+                        let u1 = acc.moments.slice(s![1, 0, ..]);
+                        let v0 = acc.moments.slice(s![0, 1, ..]);
+                        let v0 = &v0 / n0;
+                        let v1 = acc.moments.slice(s![1, 1, ..]);
+                        let v1 = &v1 / n1;
+                        let t = (&u0 - &u1)
+                            / (&v0.mapv(|x| x / n0) + &v1.mapv(|x| x / n1)).mapv(|x| f64::sqrt(x));
+                        ttest.assign(&t);
+                    } else if d == 2 {
+                        let u0 = acc.moments.slice(s![0, 1, ..]).mapv(|x| x / n0);
+                        let u1 = acc.moments.slice(s![1, 1, ..]).mapv(|x| x / n1);
+
+                        let v0 = acc.moments.slice(s![0, 3, ..]).mapv(|x| x / n0);
+                        let v1 = acc.moments.slice(s![1, 3, ..]).mapv(|x| x / n1);
+                        let v0 = &v0 - &u0.mapv(|x| x.powi(2));
+                        let v1 = &v1 - &u1.mapv(|x| x.powi(2));
+                        let t = (&u0 - &u1)
+                            / (&v0.mapv(|x| x / n0) + &v1.mapv(|x| x / n1)).mapv(|x| f64::sqrt(x));
+                        ttest.assign(&t);
+                    } else {
+                        let u0 = &(acc.moments.slice(s![0, d - 1, ..])).mapv(|x| x / n0)
+                            / &(acc
+                                .moments
+                                .slice(s![0, 1, ..])
+                                .mapv(|x| (x / n0).powf(d as f64 / 2.0)));
+
+                        let u1 = &(acc.moments.slice(s![1, d - 1, ..])).mapv(|x| x / n1)
+                            / &(acc
+                                .moments
+                                .slice(s![1, 1, ..])
+                                .mapv(|x| (x / n1).powf(d as f64 / 2.0)));
+
+                        let v0 = (&acc.moments.slice(s![0, (2 * d) - 1, ..]).mapv(|x| x / n0)
+                            - &acc
+                                .moments
+                                .slice(s![0, d - 1, ..])
+                                .mapv(|x| (x / n0).powi(2)))
+                            / &(acc
+                                .moments
+                                .slice(s![0, 1, ..])
+                                .mapv(|x| (x / n0).powi(d as i32)));
+
+                        let v1 = (&acc.moments.slice(s![1, (2 * d) - 1, ..]).mapv(|x| x / n1)
+                            - &acc
+                                .moments
+                                .slice(s![1, d - 1, ..])
+                                .mapv(|x| (x / n1).powi(2)))
+                            / &(acc
+                                .moments
+                                .slice(s![1, 1, ..])
+                                .mapv(|x| (x / n1).powi(d as i32)));
+                        let t = (&u0 - &u1)
+                            / (&v0.mapv(|x| x / n0) + &v1.mapv(|x| x / n1)).mapv(|x| f64::sqrt(x));
+                        ttest.assign(&t);
+                    }
+                });
+        });
         return ttest;
     }
 }
+
+#[inline(never)]
+pub fn prod_update(vec: &mut [f64], pv: &[f64], d: &[f64], cst: f64) {
+    izip!(vec.iter_mut(), pv.into_iter(), d.into_iter()).for_each(|(v, p, d)| *v += cst * *p * *d);
+}
+#[inline(never)]
+pub fn prod_update_add(vec: &mut [f64], pv: &[f64], cst: f64) {
+    izip!(vec.iter_mut(), pv.into_iter()).for_each(|(v, p)| *v += cst * *p);
+}
+#[inline(never)]
+pub fn prod_update_single(vec: &mut [f64], pv: &[f64], d: &[f64]) {
+    izip!(vec.iter_mut(), pv.into_iter(), d.into_iter()).for_each(|(v, p, d)| *v = *p * *d);
+}
+
+#[inline(never)]
+pub fn gen_delta(m: &mut [f64], d: &mut [f64], t: &[i16], poi: &[u64], n: f64) {
+    izip!(m.iter_mut(), d.iter_mut())
+        .enumerate()
+        .for_each(|(i, (m, d))| {
+            *d = t[poi[i as usize] as usize] as f64 - *m;
+            *m += (*d) / (n);
+        });
+}
+
+#[inline(never)]
+pub fn acc_sum(m: &mut ArrayViewMut1<i64>, t: &ArrayView1<i16>) {
+    m.zip_mut_with(t, |m, t| *m += *t as i64);
+}
+
+#[inline(always)]
+pub fn pow_and_sum(m: ArrayViewMut1<f64>, pow: ArrayViewMut1<f64>, t: ArrayView1<f64>) {
+    let m = m.into_slice().unwrap();
+    let pow = pow.into_slice().unwrap();
+    let t = t.to_slice().unwrap();
+    izip!(m.iter_mut(), pow.iter_mut(), t.iter()).for_each(|(m, pow, t)| {
+        *pow *= *t;
+        *m += *pow;
+    });
+}
diff --git a/src/scalib_ext/scalib/tests/lda.rs b/src/scalib_ext/scalib/tests/lda.rs
index 034f3c03..f6babd29 100644
--- a/src/scalib_ext/scalib/tests/lda.rs
+++ b/src/scalib_ext/scalib/tests/lda.rs
@@ -5,7 +5,6 @@ use ndarray_stats::QuantileExt;
 use scalib::lda;
 
 extern crate openblas_src as _;
-
 #[test]
 fn fit_lda() {
     let ns = 100;
diff --git a/src/scalib_ext/scalib/tests/multivarcs.rs b/src/scalib_ext/scalib/tests/multivarcs.rs
new file mode 100644
index 00000000..9d8305de
--- /dev/null
+++ b/src/scalib_ext/scalib/tests/multivarcs.rs
@@ -0,0 +1,173 @@
+// This files tests the internal representations of multivarcs.
+//
+// It first computes the central sums with multivarcs and
+// the compare the results with a reference 2 passes algorithm.
+//
+// It tests for:
+//   - different d's
+//   - feed multivarcs with multiple batches
+//   - number of samples (ns)
+//   - number of traces nc
+use itertools::izip;
+use ndarray::{s, Array1, Array2, Array3, ArrayView1, ArrayView2, ArrayView3, Axis};
+use ndarray_rand::rand_distr::Uniform;
+use ndarray_rand::RandomExt;
+use scalib::mttest;
+extern crate approx;
+
+#[test]
+fn multivarcs_simple() {
+    let order = 2;
+    let ns = 100;
+    let n = 50;
+    let nc = 2;
+    let npois = 20;
+
+    let traces = Array2::<i16>::random((n, ns as usize), Uniform::new(0, ns));
+    let pois = Array2::<u32>::random((order, npois), Uniform::new(0, ns as u32));
+    let y = Array1::<u16>::random((n,), Uniform::new(0, nc));
+
+    let mut ttacc = mttest::MultivarCSAcc::new(pois.view(), order);
+    ttacc.update(traces.view(), y.view());
+    let cs: Array3<f64> = ttacc.cs.to_owned();
+
+    test_cs(
+        traces.view(),
+        y.view(),
+        cs.view(),
+        ttacc.mean.view(),
+        pois.view(),
+        ttacc.combis,
+        ns as usize,
+        nc as usize,
+    );
+}
+
+#[test]
+fn multivarcs_simple_chuncks() {
+    let order = 2;
+    let ns = 100;
+    let n = 40;
+    let nc = 2;
+    let npois = 20;
+    let step = 10;
+
+    let traces = Array2::<i16>::random((n, ns as usize), Uniform::new(0, ns));
+    let pois = Array2::<u32>::random((order, npois), Uniform::new(0, ns as u32));
+    let y = Array1::<u16>::random((n,), Uniform::new(0, nc));
+
+    let mut ttacc = mttest::MultivarCSAcc::new(pois.view(), order);
+
+    izip!(
+        traces.axis_chunks_iter(Axis(0), step),
+        y.axis_chunks_iter(Axis(0), step)
+    )
+    .for_each(|(traces, y)| {
+        ttacc.update(traces, y);
+    });
+
+    let cs: Array3<f64> = ttacc.cs.to_owned();
+    test_cs(
+        traces.view(),
+        y.view(),
+        cs.view(),
+        ttacc.mean.view(),
+        pois.view(),
+        ttacc.combis,
+        ns as usize,
+        nc as usize,
+    );
+}
+#[test]
+fn multivarcs_simple_chuncks_step1() {
+    let order = 2;
+    let ns = 100;
+    let n = 40;
+    let nc = 2;
+    let npois = 20;
+    let step = 1;
+
+    let traces = Array2::<i16>::random((n, ns as usize), Uniform::new(0, ns));
+    let pois = Array2::<u32>::random((order, npois), Uniform::new(0, ns as u32));
+    let y = Array1::<u16>::random((n,), Uniform::new(0, nc));
+
+    let mut ttacc = mttest::MultivarCSAcc::new(pois.view(), order);
+
+    izip!(
+        traces.axis_chunks_iter(Axis(0), step),
+        y.axis_chunks_iter(Axis(0), step)
+    )
+    .for_each(|(traces, y)| {
+        ttacc.update(traces, y);
+    });
+
+    let cs: Array3<f64> = ttacc.cs.to_owned();
+    test_cs(
+        traces.view(),
+        y.view(),
+        cs.view(),
+        ttacc.mean.view(),
+        pois.view(),
+        ttacc.combis,
+        ns as usize,
+        nc as usize,
+    );
+}
+fn test_cs(
+    traces: ArrayView2<i16>,
+    y: ArrayView1<u16>,
+    cs: ArrayView3<f64>,
+    mean_ref: ArrayView3<f64>,
+    pois: ArrayView2<u32>,
+    combi: Vec<Vec<usize>>,
+    ns: usize,
+    nc: usize,
+) {
+    let epsi = 1e-4;
+    for i in 0..nc {
+        let traces_0: Vec<ArrayView1<i16>> = izip!(traces.outer_iter(), y.iter())
+            .filter(|(_, y)| **y == i as u16)
+            .map(|(t, _)| t)
+            .collect();
+
+        let mut t0 = Array2::<i16>::zeros((traces_0.len(), ns as usize));
+        for (i, t) in traces_0.iter().enumerate() {
+            t0.slice_mut(s![i, ..]).assign(t);
+        }
+        let mut t0 = t0.mapv(|x| x as f64);
+        let n_traces = t0.shape()[0];
+        let mean: Array1<f64> = t0.mean_axis(Axis(0)).unwrap();
+
+        izip!(
+            mean_ref.slice(s![i as usize, .., ..]).outer_iter(),
+            pois.outer_iter()
+        )
+        .for_each(|(mean_ref, pois)| {
+            println!("New vector {:#?} {:#?}", mean_ref.shape(), pois.shape());
+            izip!(mean_ref, pois).for_each(|(x, poi)| {
+                println!("{} {}", x, mean[*poi as usize]);
+                approx::assert_abs_diff_eq!(*x, mean[*poi as usize], epsilon = epsi);
+            });
+        });
+
+        //centering
+        t0.axis_iter_mut(Axis(0)).for_each(|mut x| x -= &mean);
+
+        // test centered sums
+        let mut tests = 0;
+        izip!(cs.slice(s![i as usize, .., ..]).outer_iter(), combi.iter()).for_each(
+            |(m, combi)| {
+                println!("New Combi {:#?}", combi);
+                izip!(m.iter(), pois.axis_iter(Axis(1))).for_each(|(m, poi)| {
+                    let mut test = Array1::<f64>::ones((n_traces,));
+                    for c in combi.iter() {
+                        test = &test * &t0.slice(s![.., poi[*c] as usize]);
+                    }
+                    println!("{} {}", test.sum(), m);
+                    approx::assert_abs_diff_eq!(test.sum(), m, epsilon = epsi);
+                    tests += 1;
+                });
+            },
+        );
+    }
+}
diff --git a/src/scalib_ext/scalib/tests/ttest.rs b/src/scalib_ext/scalib/tests/ttest.rs
new file mode 100644
index 00000000..8f09585d
--- /dev/null
+++ b/src/scalib_ext/scalib/tests/ttest.rs
@@ -0,0 +1,184 @@
+use itertools::izip;
+use ndarray::{s, Array1, Array2, Array3, ArrayView1, ArrayView2, ArrayView3, Axis};
+use ndarray_rand::rand_distr::Uniform;
+use ndarray_rand::RandomExt;
+//use ndarray_stats::QuantileExt;
+use scalib::ttest;
+
+#[test]
+fn ttestacc_simple() {
+    let order = 3;
+    let ns = 100;
+    let n = 100;
+    let nc = 2;
+
+    let traces = Array2::<i16>::random((n, ns as usize), Uniform::new(0, ns));
+    let y = Array1::<u16>::random((n,), Uniform::new(0, nc));
+
+    let mut ttacc = ttest::UniCSAcc::new(ns as usize, order, nc as usize);
+    ttacc.update(traces.view(), y.view());
+    let moments: Array3<f64> = ttacc.moments.to_owned();
+
+    test_moments(
+        traces.view(),
+        y.view(),
+        moments.view(),
+        order as usize,
+        ns as usize,
+        nc as usize,
+    );
+}
+
+#[test]
+fn ttestacc_simple_chuncks() {
+    let order = 5;
+    let ns = 100;
+    let n = 2000;
+    let step = 10;
+    let nc = 2;
+
+    let traces = Array2::<i16>::random((n, ns as usize), Uniform::new(0, ns));
+    let y = Array1::<u16>::random((n,), Uniform::new(0, nc));
+
+    // perform ttacc
+    let mut ttacc = ttest::UniCSAcc::new(ns as usize, order, nc as usize);
+    izip!(
+        traces.axis_chunks_iter(Axis(0), step),
+        y.axis_chunks_iter(Axis(0), step)
+    )
+    .for_each(|(traces, y)| {
+        ttacc.update(traces, y);
+    });
+    let moments: Array3<f64> = ttacc.moments.to_owned();
+    test_moments(
+        traces.view(),
+        y.view(),
+        moments.view(),
+        order as usize,
+        ns as usize,
+        nc as usize,
+    );
+}
+
+#[test]
+fn ttestacc_merge_and_reset_chuncks() {
+    let order = 5;
+    let ns = 100;
+    let n = 2000;
+    let step = 10;
+    let nc = 2;
+
+    let traces = Array2::<i16>::random((n, ns as usize), Uniform::new(0, ns));
+    let y = Array1::<u16>::random((n,), Uniform::new(0, nc));
+
+    // perform ttacc
+    let mut ttacc = ttest::UniCSAcc::new(ns as usize, order, nc as usize);
+    let mut tmp = ttest::UniCSAcc::new(ns as usize, order, nc as usize);
+    izip!(
+        traces.axis_chunks_iter(Axis(0), step),
+        y.axis_chunks_iter(Axis(0), step)
+    )
+    .for_each(|(traces, y)| {
+        tmp.update(traces, y);
+        ttacc.merge(&tmp);
+        tmp.reset();
+    });
+
+    let moments: Array3<f64> = ttacc.moments.to_owned();
+    test_moments(
+        traces.view(),
+        y.view(),
+        moments.view(),
+        order as usize,
+        ns as usize,
+        nc as usize,
+    );
+}
+#[test]
+fn ttestacc_merge_chuncks() {
+    let order = 5;
+    let ns = 100;
+    let n = 2000;
+    let step = 10;
+    let nc = 2;
+
+    let traces = Array2::<i16>::random((n, ns as usize), Uniform::new(0, ns));
+    let y = Array1::<u16>::random((n,), Uniform::new(0, nc));
+
+    // perform ttacc
+    let mut ttacc = ttest::UniCSAcc::new(ns as usize, order, nc as usize);
+    izip!(
+        traces.axis_chunks_iter(Axis(0), step),
+        y.axis_chunks_iter(Axis(0), step)
+    )
+    .for_each(|(traces, y)| {
+        let mut tmp = ttest::UniCSAcc::new(ns as usize, order, nc as usize);
+        tmp.update(traces, y);
+        ttacc.merge(&tmp);
+    });
+
+    let moments: Array3<f64> = ttacc.moments.to_owned();
+    test_moments(
+        traces.view(),
+        y.view(),
+        moments.view(),
+        order as usize,
+        ns as usize,
+        nc as usize,
+    );
+}
+fn test_moments(
+    traces: ArrayView2<i16>,
+    y: ArrayView1<u16>,
+    moments: ArrayView3<f64>,
+    order: usize,
+    ns: usize,
+    nc: usize,
+) {
+    // Testing on set 0
+    // t0
+    let epsi = 1e-4;
+    for i in 0..nc {
+        let traces_0: Vec<ArrayView1<i16>> = izip!(traces.outer_iter(), y.iter())
+            .filter(|(_, y)| **y == i as u16)
+            .map(|(t, _)| t)
+            .collect();
+        let mut t0 = Array2::<i16>::zeros((traces_0.len(), ns as usize));
+        for (i, t) in traces_0.iter().enumerate() {
+            t0.slice_mut(s![i, ..]).assign(t);
+        }
+        let mut t0 = t0.mapv(|x| x as f64);
+        let mean: Array1<f64> = t0.mean_axis(Axis(0)).unwrap();
+        //centering
+        t0.axis_iter_mut(Axis(0)).for_each(|mut x| x -= &mean);
+
+        // testing the mean:
+        let mm0 = moments.slice(s![i as usize, .., ..]);
+        let test: Array1<f64> = mm0.slice(s![0 as usize, ..]).to_owned();
+        let mut dif = &mean - &test;
+        dif.mapv_inplace(|x| x.abs());
+        for d in dif.iter() {
+            assert!(*d < epsi, "failed on mean {}", *d);
+        }
+
+        // test centered sums
+        for j in 2..(order + 1) {
+            let cs = t0.mapv(|x| x.powi(j as i32));
+            let reference = cs.sum_axis(Axis(0));
+
+            let test: Array1<f64> = mm0.slice(s![j - 1 as usize, ..]).to_owned();
+            let mut dif = &reference - &test;
+            dif.mapv_inplace(|x| x.abs());
+            for d in dif.iter() {
+                assert!(
+                    *d < epsi,
+                    "failed on order {} {} {} {}",
+                    j,
+                    d,
+                    test[0],
+                    reference[0]
+                );
+            }
+        }
+    }
+}
diff --git a/tests/mttest_test.py b/tests/mttest_test.py
new file mode 100644
index 00000000..f6a52ff3
--- /dev/null
+++ b/tests/mttest_test.py
@@ -0,0 +1,205 @@
+"""
+    This script feeds MTtest object and compare the outcome with the reference
+    implementation. 
+    
+    The tests covers:
+        - various d
+        - Usage of incremental API. 
+        - Number of pois
+        - Number of traces
+"""
+import pytest
+from scalib.metrics import MTtest
+import numpy as np
+
+
+def reference(traces, x, d, pois):
+    I0 = np.where(x == 0)[0]
+    I1 = np.where(x == 1)[0]
+
+    n0 = len(I0)
+    n1 = len(I1)
+
+    t0 = traces[I0]
+    t1 = traces[I1]
+
+    t0 = t0 - np.mean(t0, axis=0)
+    t1 = t1 - np.mean(t1, axis=0)
+    if d > 2:
+        t0 = t0 / np.std(t0, axis=0)
+        t1 = t1 / np.std(t1, axis=0)
+
+    t0_mult = 1
+    t1_mult = 1
+    for d in range(d):
+        t0_mult *= t0[:, pois[d, :]]
+        t1_mult *= t1[:, pois[d, :]]
+
+    print(t0_mult.shape)
+    u0 = np.mean(t0_mult, axis=0)
+    u1 = np.mean(t1_mult, axis=0)
+
+    v0 = np.var(t0_mult, axis=0)
+    v1 = np.var(t1_mult, axis=0)
+
+    return (u0 - u1) / np.sqrt((v0 / n0) + (v1 / n1))
+
+
+def test_ttest_d2_many_3():
+    ns = 10237
+    npois = 12342
+    d = 2
+    nc = 2
+    n = 1 << 8
+
+    m = np.random.randint(0, 2, (nc, ns))
+    traces = np.random.randint(0, 10, (n, ns), dtype=np.int16)
+    labels = np.random.randint(0, nc, n, dtype=np.uint16)
+    pois = np.random.randint(0, ns, (d, npois), dtype=np.uint32)
+    traces += m[labels]
+
+    t_ref = reference(traces, labels, d, pois)
+    ttest = MTtest(d, pois)
+    ttest.fit_u(traces, labels)
+    t = ttest.get_ttest()
+    assert np.allclose(t_ref, t, rtol=1e-3)
+
+
+def test_ttest_d2_many_2():
+    ns = 99
+    npois = 12329
+    d = 2
+    nc = 2
+    n = 1 << 8
+
+    m = np.random.randint(0, 2, (nc, ns))
+    traces = np.random.randint(0, 10, (n, ns), dtype=np.int16)
+    labels = np.random.randint(0, nc, n, dtype=np.uint16)
+    pois = np.random.randint(0, ns, (d, npois), dtype=np.uint32)
+    traces += m[labels]
+
+    t_ref = reference(traces, labels, d, pois)
+    ttest = MTtest(d, pois)
+    ttest.fit_u(traces, labels)
+    t = ttest.get_ttest()
+    assert np.allclose(t_ref, t, rtol=1e-3)
+
+
+def test_ttest_d2_many():
+    ns = 10
+    npois = 6
+    d = 2
+    nc = 2
+    n = 1 << 17
+
+    m = np.random.randint(0, 2, (nc, ns))
+    traces = np.random.randint(0, 10, (n, ns), dtype=np.int16)
+    labels = np.random.randint(0, nc, n, dtype=np.uint16)
+    pois = np.random.randint(0, ns, (d, npois), dtype=np.uint32)
+    traces += m[labels]
+
+    t_ref = reference(traces, labels, d, pois)
+    ttest = MTtest(d, pois)
+    ttest.fit_u(traces, labels)
+    t = ttest.get_ttest()
+    assert np.allclose(t_ref, t, rtol=1e-3)
+
+
+def test_ttest_d2():
+    ns = 100
+    npois = 6302
+    d = 2
+    nc = 2
+    n = 200
+
+    m = np.random.randint(0, 2, (nc, ns))
+    traces = np.random.randint(0, 10, (n, ns), dtype=np.int16)
+    labels = np.random.randint(0, nc, n, dtype=np.uint16)
+    pois = np.random.randint(0, ns, (d, npois), dtype=np.uint32)
+    traces += m[labels]
+
+    t_ref = reference(traces, labels, d, pois)
+    ttest = MTtest(d, pois)
+    ttest.fit_u(traces, labels)
+    t = ttest.get_ttest()
+    assert np.allclose(t_ref, t, rtol=1e-3)
+
+
+def test_ttest_d3():
+    ns = 100
+    npois = 10
+    d = 3
+    nc = 2
+    n = 200
+
+    m = np.random.randint(0, 2, (nc, ns))
+    traces = np.random.randint(0, 10, (n, ns), dtype=np.int16)
+    labels = np.random.randint(0, nc, n, dtype=np.uint16)
+    pois = np.random.randint(0, ns, (d, npois), dtype=np.uint32)
+    traces += m[labels]
+
+    t_ref = reference(traces, labels, d, pois)
+    ttest = MTtest(d, pois)
+    ttest.fit_u(traces, labels)
+    t = ttest.get_ttest()
+    assert np.allclose(t_ref, t, rtol=1e-3)
+
+
+def test_ttest_d3_multiple_fit():
+    ns = 1043
+    d = 3
+    nc = 2
+    npois = 10021
+    n = 40
+
+    m = np.random.randint(0, 2, (nc, ns))
+    traces = np.random.randint(0, 10, (n, ns), dtype=np.int16)
+    labels = np.random.randint(0, nc, n, dtype=np.uint16)
+    pois = np.random.randint(0, ns, (d, npois), dtype=np.uint32)
+
+    t_ref = reference(traces, labels, d, pois)
+    ttest = MTtest(d, pois)
+    for i in range(0, n, 10):
+        ttest.fit_u(traces[i : i + 10, :], labels[i : i + 10])
+    t = ttest.get_ttest()
+    assert np.allclose(t_ref, t, rtol=1e-3)
+
+
+def test_ttest_d3_multiple_fit():
+    ns = 1043
+    d = 3
+    nc = 2
+    npois = 10021
+    n = 40
+
+    m = np.random.randint(0, 2, (nc, ns))
+    traces = np.random.randint(0, 10, (n, ns), dtype=np.int16)
+    labels = np.random.randint(0, nc, n, dtype=np.uint16)
+    pois = np.random.randint(0, ns, (d, npois), dtype=np.uint32)
+
+    t_ref = reference(traces, labels, d, pois)
+    ttest = MTtest(d, pois)
+    for i in range(0, n, 10):
+        ttest.fit_u(traces[i : i + 10, :], labels[i : i + 10])
+    t = ttest.get_ttest()
+    assert np.allclose(t_ref, t, rtol=1e-3)
+
+
+def test_ttest_d2_multiple_fit():
+    ns = 1043
+    d = 2
+    nc = 2
+    npois = 10021
+    n = 40
+
+    m = np.random.randint(0, 2, (nc, ns))
+    traces = np.random.randint(-4231, 4214, (n, ns), dtype=np.int16)
+    labels = np.random.randint(0, nc, n, dtype=np.uint16)
+    pois = np.random.randint(0, ns, (d, npois), dtype=np.uint32)
+
+    t_ref = reference(traces, labels, d, pois)
+    ttest = MTtest(d, pois)
+    for i in range(0, n, 10):
+        ttest.fit_u(traces[i : i + 10, :], labels[i : i + 10])
+    t = ttest.get_ttest()
+    assert np.allclose(t_ref, t, rtol=1e-3)
diff --git a/tests/test_ttest.py b/tests/test_ttest.py
index 27948ea8..69fd9260 100644
--- a/tests/test_ttest.py
+++ b/tests/test_ttest.py
@@ -5,7 +5,6 @@
 
 
 def reference(traces, x, D):
-
     CM0 = np.zeros((D * 2, traces.shape[1]))
     I0 = np.where(x == 0)[0]
     u0 = np.mean(traces[I0, :], axis=0)
@@ -79,7 +78,7 @@ def test_ttest_d2():
 
 
 def test_ttest_d6_multiple_fit():
-    ns = 100
+    ns = 12300
     d = 6
     nc = 2
     n = 200
@@ -98,7 +97,7 @@ def test_ttest_d6_multiple_fit():
 
 
 def test_ttest_d6():
-    ns = 100
+    ns = 942
     d = 6
     nc = 2
     n = 200