I have found that the implementation of Barrett reduction can give incorrect results in some cases. Here is an example:

use concrete_ntt::prime32::Plan;

fn main() {
    let p: u32 = 0x7fe0_1001;
    let polynomial_size = 1024;
    let plan = Plan::try_new(polynomial_size, p).unwrap();

    let value = 0x6e63593a;
    let mut acc = [0u32; 8];
    let input = [value, 0, 0, 0, 0, 0, 0, 0];

    plan.mul_accumulate(&mut acc, &input, &input);

    let expected = (u64::from(value) * u64::from(value) % u64::from(p)) as u32;
    println!("acc[0] = {}", acc[0]);
    assert_eq!(acc[0], expected);
}

The output is:

acc[0] = 360086499
thread 'main' panicked at examples/reduction.rs:16:5:
assertion `left == right` failed
  left: 360086499
 right: 364272609

Note that 364272609 - 360086499 = 4186110 = 2 * (0x8000_0000 - 0x7fe0_1001). Performing two conditional subtractions at the conclusion of the algorithm, rather than one, would have given the correct result. (Unfortunately, there are not enough bits in the intermediate result to see that two conditional subtractions are necessary.)

The possible need for two conditional subtractions was noted in Barrett's paper ("Our calculations show that the result x so obtained will always be in the range 0 to 3M - 1") and is also captured in https://eprint.iacr.org/2015/785 ("approximates the result $c = d \bmod n$ by a quasi-reduced number $c + \epsilon n$ where $0 \le \epsilon \le 2$"), but is omitted by the presentation in https://arxiv.org/abs/2103.16400.

If step 1 of the algorithm is "d >> X", I calculate that one subtraction will be sufficient if:

$$\frac{d}{2^L} + \frac{2^X}{p} < 1$$

But, for a 31-bit prime, that increases the size of the multiplication beyond 32 bits.