arnaucube
arnaucube/phantom-zonePublic Self-hosted repository
phantom-zone/src/rgsw.rs
src/rgsw.rsRaw
   1use std::{
   2    clone,
   3    fmt::Debug,
   4    marker::PhantomData,
   5    ops::{Div, Neg, Sub},
   6};
   7
   8use itertools::{izip, Itertools};
   9use num_traits::{PrimInt, Signed, ToPrimitive, Zero};
  10
  11use crate::{
  12    backend::{ArithmeticOps, GetModulus, Modulus, VectorOps},
  13    decomposer::{self, Decomposer, RlweDecomposer},
  14    ntt::{self, Ntt, NttInit},
  15    random::{
  16        DefaultSecureRng, NewWithSeed, RandomElementInModulus, RandomFill,
  17        RandomFillGaussianInModulus, RandomFillUniformInModulus,
  18    },
  19    utils::{fill_random_ternary_secret_with_hamming_weight, TryConvertFrom1, WithLocal},
  20    Matrix, MatrixEntity, MatrixMut, Row, RowEntity, RowMut, Secret,
  21};
  22
  23pub struct SeededAutoKey<M, S, Mod>
  24where
  25    M: Matrix,
  26{
  27    data: M,
  28    seed: S,
  29    modulus: Mod,
  30}
  31
  32impl<M: Matrix + MatrixEntity, S, Mod: Modulus<Element = M::MatElement>> SeededAutoKey<M, S, Mod> {
  33    fn empty<D: Decomposer>(ring_size: usize, auto_decomposer: &D, seed: S, modulus: Mod) -> Self {
  34        SeededAutoKey {
  35            data: M::zeros(auto_decomposer.decomposition_count(), ring_size),
  36            seed,
  37            modulus,
  38        }
  39    }
  40}
  41
  42pub struct AutoKeyEvaluationDomain<M, R, N> {
  43    data: M,
  44    _phantom: PhantomData<(R, N)>,
  45}
  46
  47impl<
  48        M: MatrixMut + MatrixEntity,
  49        Mod: Modulus<Element = M::MatElement> + Clone,
  50        R: RandomFillUniformInModulus<[M::MatElement], Mod> + NewWithSeed,
  51        N: NttInit<Mod> + Ntt<Element = M::MatElement>,
  52    > From<&SeededAutoKey<M, R::Seed, Mod>> for AutoKeyEvaluationDomain<M, R, N>
  53where
  54    <M as Matrix>::R: RowMut,
  55    M::MatElement: Copy,
  56    R::Seed: Clone,
  57{
  58    fn from(value: &SeededAutoKey<M, R::Seed, Mod>) -> Self {
  59        let (d, ring_size) = value.data.dimension();
  60        let mut data = M::zeros(2 * d, ring_size);
  61
  62        // sample RLWE'_A(-s(X^k))
  63        let mut p_rng = R::new_with_seed(value.seed.clone());
  64        data.iter_rows_mut().take(d).for_each(|r| {
  65            RandomFillUniformInModulus::random_fill(&mut p_rng, &value.modulus, r.as_mut());
  66        });
  67
  68        // copy over RLWE'_B(-s(X^k))
  69        izip!(data.iter_rows_mut().skip(d), value.data.iter_rows()).for_each(|(to_r, from_r)| {
  70            to_r.as_mut().copy_from_slice(from_r.as_ref());
  71        });
  72
  73        // send RLWE'(-s(X^k)) polynomials to evaluation domain
  74        let ntt_op = N::new(&value.modulus, ring_size);
  75        data.iter_rows_mut()
  76            .for_each(|r| ntt_op.forward(r.as_mut()));
  77
  78        AutoKeyEvaluationDomain {
  79            data,
  80            _phantom: PhantomData,
  81        }
  82    }
  83}
  84
  85pub struct RgswCiphertext<M: Matrix, Mod> {
  86    /// Rgsw ciphertext polynomials
  87    pub(crate) data: M,
  88    modulus: Mod,
  89    /// Decomposition for RLWE part A
  90    d_a: usize,
  91    /// Decomposition for RLWE part B
  92    d_b: usize,
  93}
  94
  95impl<M: MatrixEntity, Mod: Modulus<Element = M::MatElement>> RgswCiphertext<M, Mod> {
  96    pub(crate) fn empty<D: RlweDecomposer>(
  97        ring_size: usize,
  98        decomposer: &D,
  99        modulus: Mod,
 100    ) -> RgswCiphertext<M, Mod> {
 101        RgswCiphertext {
 102            data: M::zeros(
 103                decomposer.a().decomposition_count() * 2 + decomposer.b().decomposition_count() * 2,
 104                ring_size,
 105            ),
 106            d_a: decomposer.a().decomposition_count(),
 107            d_b: decomposer.b().decomposition_count(),
 108            modulus,
 109        }
 110    }
 111}
 112
 113pub struct SeededRgswCiphertext<M, S, Mod>
 114where
 115    M: Matrix,
 116{
 117    pub(crate) data: M,
 118    seed: S,
 119    modulus: Mod,
 120    /// Decomposition for RLWE part A
 121    d_a: usize,
 122    /// Decomposition for RLWE part B
 123    d_b: usize,
 124}
 125
 126impl<M: Matrix + MatrixEntity, S, Mod> SeededRgswCiphertext<M, S, Mod> {
 127    pub(crate) fn empty<D: RlweDecomposer>(
 128        ring_size: usize,
 129        decomposer: &D,
 130        seed: S,
 131        modulus: Mod,
 132    ) -> SeededRgswCiphertext<M, S, Mod> {
 133        SeededRgswCiphertext {
 134            data: M::zeros(
 135                decomposer.a().decomposition_count() * 2 + decomposer.b().decomposition_count(),
 136                ring_size,
 137            ),
 138            seed,
 139            modulus,
 140            d_a: decomposer.a().decomposition_count(),
 141            d_b: decomposer.b().decomposition_count(),
 142        }
 143    }
 144}
 145
 146impl<M: Debug + Matrix, S: Debug, Mod: Debug> Debug for SeededRgswCiphertext<M, S, Mod>
 147where
 148    M::MatElement: Debug,
 149{
 150    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
 151        f.debug_struct("SeededRgswCiphertext")
 152            .field("data", &self.data)
 153            .field("seed", &self.seed)
 154            .field("modulus", &self.modulus)
 155            .finish()
 156    }
 157}
 158
 159pub struct RgswCiphertextEvaluationDomain<M, R, N> {
 160    pub(crate) data: M,
 161    _phantom: PhantomData<(R, N)>,
 162}
 163
 164impl<
 165        M: MatrixMut + MatrixEntity,
 166        Mod: Modulus<Element = M::MatElement>,
 167        R: NewWithSeed + RandomFillUniformInModulus<[M::MatElement], Mod>,
 168        N: NttInit<Mod> + Ntt<Element = M::MatElement> + Debug,
 169    > From<&SeededRgswCiphertext<M, R::Seed, Mod>> for RgswCiphertextEvaluationDomain<M, R, N>
 170where
 171    <M as Matrix>::R: RowMut,
 172    M::MatElement: Copy,
 173    R::Seed: Clone,
 174    M: Debug,
 175{
 176    fn from(value: &SeededRgswCiphertext<M, R::Seed, Mod>) -> Self {
 177        let mut data = M::zeros(value.d_a * 2 + value.d_b * 2, value.data.dimension().1);
 178
 179        // copy RLWE'(-sm)
 180        izip!(
 181            data.iter_rows_mut().take(value.d_a * 2),
 182            value.data.iter_rows().take(value.d_a * 2)
 183        )
 184        .for_each(|(to_ri, from_ri)| {
 185            to_ri.as_mut().copy_from_slice(from_ri.as_ref());
 186        });
 187
 188        // sample A polynomials of RLWE'(m) - RLWE'A(m)
 189        // TODO(Jay): Do we want to be generic over RandomGenerator used here? I think
 190        // not.
 191        let mut p_rng = R::new_with_seed(value.seed.clone());
 192        izip!(data.iter_rows_mut().skip(value.d_a * 2).take(value.d_b * 1))
 193            .for_each(|ri| p_rng.random_fill(&value.modulus, ri.as_mut()));
 194
 195        // RLWE'_B(m)
 196        izip!(
 197            data.iter_rows_mut().skip(value.d_a * 2 + value.d_b),
 198            value.data.iter_rows().skip(value.d_a * 2)
 199        )
 200        .for_each(|(to_ri, from_ri)| {
 201            to_ri.as_mut().copy_from_slice(from_ri.as_ref());
 202        });
 203
 204        // Send polynomials to evaluation domain
 205        let ring_size = data.dimension().1;
 206        let nttop = N::new(&value.modulus, ring_size);
 207        data.iter_rows_mut()
 208            .for_each(|ri| nttop.forward(ri.as_mut()));
 209
 210        Self {
 211            data: data,
 212            _phantom: PhantomData,
 213        }
 214    }
 215}
 216
 217impl<
 218        M: MatrixMut + MatrixEntity,
 219        Mod: Modulus<Element = M::MatElement>,
 220        R,
 221        N: NttInit<Mod> + Ntt<Element = M::MatElement>,
 222    > From<&RgswCiphertext<M, Mod>> for RgswCiphertextEvaluationDomain<M, R, N>
 223where
 224    <M as Matrix>::R: RowMut,
 225    M::MatElement: Copy,
 226    M: Debug,
 227{
 228    fn from(value: &RgswCiphertext<M, Mod>) -> Self {
 229        let mut data = M::zeros(value.d_a * 2 + value.d_b * 2, value.data.dimension().1);
 230
 231        // copy RLWE'(-sm)
 232        izip!(
 233            data.iter_rows_mut().take(value.d_a * 2),
 234            value.data.iter_rows().take(value.d_a * 2)
 235        )
 236        .for_each(|(to_ri, from_ri)| {
 237            to_ri.as_mut().copy_from_slice(from_ri.as_ref());
 238        });
 239
 240        // copy RLWE'(m)
 241        izip!(
 242            data.iter_rows_mut().skip(value.d_a * 2),
 243            value.data.iter_rows().skip(value.d_a * 2)
 244        )
 245        .for_each(|(to_ri, from_ri)| {
 246            to_ri.as_mut().copy_from_slice(from_ri.as_ref());
 247        });
 248
 249        // Send polynomials to evaluation domain
 250        let ring_size = data.dimension().1;
 251        let nttop = N::new(&value.modulus, ring_size);
 252        data.iter_rows_mut()
 253            .for_each(|ri| nttop.forward(ri.as_mut()));
 254
 255        Self {
 256            data: data,
 257            _phantom: PhantomData,
 258        }
 259    }
 260}
 261
 262impl<M: Debug, R, N> Debug for RgswCiphertextEvaluationDomain<M, R, N> {
 263    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
 264        f.debug_struct("RgswCiphertextEvaluationDomain")
 265            .field("data", &self.data)
 266            .field("_phantom", &self._phantom)
 267            .finish()
 268    }
 269}
 270
 271impl<M: Matrix, R, N> Matrix for RgswCiphertextEvaluationDomain<M, R, N> {
 272    type MatElement = M::MatElement;
 273    type R = M::R;
 274
 275    fn dimension(&self) -> (usize, usize) {
 276        self.data.dimension()
 277    }
 278
 279    fn fits(&self, row: usize, col: usize) -> bool {
 280        self.data.fits(row, col)
 281    }
 282}
 283
 284impl<M: Matrix, R, N> AsRef<[M::R]> for RgswCiphertextEvaluationDomain<M, R, N> {
 285    fn as_ref(&self) -> &[M::R] {
 286        self.data.as_ref()
 287    }
 288}
 289
 290pub struct SeededRlweCiphertext<R, S, Mod> {
 291    pub(crate) data: R,
 292    pub(crate) seed: S,
 293    pub(crate) modulus: Mod,
 294}
 295
 296impl<R: RowEntity, S, Mod> SeededRlweCiphertext<R, S, Mod> {
 297    pub(crate) fn empty(ring_size: usize, seed: S, modulus: Mod) -> Self {
 298        SeededRlweCiphertext {
 299            data: R::zeros(ring_size),
 300            seed,
 301            modulus,
 302        }
 303    }
 304}
 305
 306pub struct RlweCiphertext<M, Rng> {
 307    pub(crate) data: M,
 308    pub(crate) is_trivial: bool,
 309    pub(crate) _phatom: PhantomData<Rng>,
 310}
 311
 312impl<M, Rng> RlweCiphertext<M, Rng> {
 313    pub(crate) fn new_trivial(data: M) -> Self {
 314        RlweCiphertext {
 315            data,
 316            is_trivial: true,
 317            _phatom: PhantomData,
 318        }
 319    }
 320}
 321
 322impl<M: Matrix, Rng> Matrix for RlweCiphertext<M, Rng> {
 323    type MatElement = M::MatElement;
 324    type R = M::R;
 325
 326    fn dimension(&self) -> (usize, usize) {
 327        self.data.dimension()
 328    }
 329
 330    fn fits(&self, row: usize, col: usize) -> bool {
 331        self.data.fits(row, col)
 332    }
 333}
 334
 335impl<M: MatrixMut, Rng> MatrixMut for RlweCiphertext<M, Rng> where <M as Matrix>::R: RowMut {}
 336
 337impl<M: Matrix, Rng> AsRef<[<M as Matrix>::R]> for RlweCiphertext<M, Rng> {
 338    fn as_ref(&self) -> &[<M as Matrix>::R] {
 339        self.data.as_ref()
 340    }
 341}
 342
 343impl<M: MatrixMut, Rng> AsMut<[<M as Matrix>::R]> for RlweCiphertext<M, Rng>
 344where
 345    <M as Matrix>::R: RowMut,
 346{
 347    fn as_mut(&mut self) -> &mut [<M as Matrix>::R] {
 348        self.data.as_mut()
 349    }
 350}
 351
 352impl<M, Rng> IsTrivial for RlweCiphertext<M, Rng> {
 353    fn is_trivial(&self) -> bool {
 354        self.is_trivial
 355    }
 356    fn set_not_trivial(&mut self) {
 357        self.is_trivial = false;
 358    }
 359}
 360
 361impl<
 362        R: Row,
 363        M: MatrixEntity<R = R, MatElement = R::Element> + MatrixMut,
 364        Rng: NewWithSeed + RandomFillUniformInModulus<[M::MatElement], Mod>,
 365        Mod: Modulus<Element = R::Element>,
 366    > From<&SeededRlweCiphertext<R, Rng::Seed, Mod>> for RlweCiphertext<M, Rng>
 367where
 368    Rng::Seed: Clone,
 369    <M as Matrix>::R: RowMut,
 370    R::Element: Copy,
 371{
 372    fn from(value: &SeededRlweCiphertext<R, Rng::Seed, Mod>) -> Self {
 373        let mut data = M::zeros(2, value.data.as_ref().len());
 374
 375        // sample a
 376        let mut p_rng = Rng::new_with_seed(value.seed.clone());
 377        RandomFillUniformInModulus::random_fill(&mut p_rng, &value.modulus, data.get_row_mut(0));
 378
 379        data.get_row_mut(1).copy_from_slice(value.data.as_ref());
 380
 381        RlweCiphertext {
 382            data,
 383            is_trivial: false,
 384            _phatom: PhantomData,
 385        }
 386    }
 387}
 388
 389pub trait IsTrivial {
 390    fn is_trivial(&self) -> bool;
 391    fn set_not_trivial(&mut self);
 392}
 393
 394pub struct SeededRlwePublicKey<Ro: Row, S> {
 395    data: Ro,
 396    seed: S,
 397    modulus: Ro::Element,
 398}
 399
 400impl<Ro: RowEntity, S> SeededRlwePublicKey<Ro, S> {
 401    pub(crate) fn empty(ring_size: usize, seed: S, modulus: Ro::Element) -> Self {
 402        Self {
 403            data: Ro::zeros(ring_size),
 404            seed,
 405            modulus,
 406        }
 407    }
 408}
 409
 410pub struct RlwePublicKey<M, R> {
 411    data: M,
 412    _phantom: PhantomData<R>,
 413}
 414
 415impl<
 416        M: MatrixMut + MatrixEntity,
 417        Rng: NewWithSeed + RandomFillUniformInModulus<[M::MatElement], M::MatElement>,
 418    > From<&SeededRlwePublicKey<M::R, Rng::Seed>> for RlwePublicKey<M, Rng>
 419where
 420    <M as Matrix>::R: RowMut,
 421    M::MatElement: Copy,
 422    Rng::Seed: Copy,
 423{
 424    fn from(value: &SeededRlwePublicKey<M::R, Rng::Seed>) -> Self {
 425        let mut data = M::zeros(2, value.data.as_ref().len());
 426
 427        // sample a
 428        let mut p_rng = Rng::new_with_seed(value.seed);
 429        RandomFillUniformInModulus::random_fill(&mut p_rng, &value.modulus, data.get_row_mut(0));
 430
 431        // copy over b
 432        data.get_row_mut(1).copy_from_slice(value.data.as_ref());
 433
 434        Self {
 435            data,
 436            _phantom: PhantomData,
 437        }
 438    }
 439}
 440
 441#[derive(Clone)]
 442pub struct RlweSecret {
 443    pub(crate) values: Vec<i32>,
 444}
 445
 446impl Secret for RlweSecret {
 447    type Element = i32;
 448    fn values(&self) -> &[Self::Element] {
 449        &self.values
 450    }
 451}
 452
 453impl RlweSecret {
 454    pub fn random(hw: usize, n: usize) -> RlweSecret {
 455        DefaultSecureRng::with_local_mut(|rng| {
 456            let mut out = vec![0i32; n];
 457            fill_random_ternary_secret_with_hamming_weight(&mut out, hw, rng);
 458
 459            RlweSecret { values: out }
 460        })
 461    }
 462}
 463
 464pub(crate) fn generate_auto_map(ring_size: usize, k: isize) -> (Vec<usize>, Vec<bool>) {
 465    assert!(k & 1 == 1, "Auto {k} must be odd");
 466
 467    let k = if k < 0 {
 468        // k is -ve, return k%(2*N)
 469        (2 * ring_size) - (k.abs() as usize % (2 * ring_size))
 470    } else {
 471        k as usize
 472    };
 473    let (auto_map_index, auto_sign_index): (Vec<usize>, Vec<bool>) = (0..ring_size)
 474        .into_iter()
 475        .map(|i| {
 476            let mut to_index = (i * k) % (2 * ring_size);
 477            let mut sign = true;
 478
 479            // wrap around. false implies negative
 480            if to_index >= ring_size {
 481                to_index = to_index - ring_size;
 482                sign = false;
 483            }
 484
 485            (to_index, sign)
 486        })
 487        .unzip();
 488    (auto_map_index, auto_sign_index)
 489}
 490
 491pub(crate) fn routine<R: RowMut, ModOp: VectorOps<Element = R::Element>>(
 492    write_to_row: &mut [R::Element],
 493    matrix_a: &[R],
 494    matrix_b: &[R],
 495    mod_op: &ModOp,
 496) {
 497    izip!(matrix_a.iter(), matrix_b.iter()).for_each(|(a, b)| {
 498        mod_op.elwise_fma_mut(write_to_row, a.as_ref(), b.as_ref());
 499    });
 500}
 501
 502/// Decomposes ring polynomial r(X) into d polynomials using decomposer into
 503/// output matrix decomp_r
 504///
 505/// Note that decomposition of r(X) requires decomposition of each of
 506/// coefficients.
 507///
 508/// - decomp_r: must have dimensions d x ring_size. i^th decomposed polynomial
 509///   will be stored at i^th row.
 510pub(crate) fn decompose_r<R: RowMut, D: Decomposer<Element = R::Element>>(
 511    r: &[R::Element],
 512    decomp_r: &mut [R],
 513    decomposer: &D,
 514) where
 515    R::Element: Copy,
 516{
 517    let ring_size = r.len();
 518    let d = decomposer.decomposition_count();
 519
 520    for ri in 0..ring_size {
 521        let el_decomposed = decomposer.decompose_to_vec(&r[ri]);
 522
 523        for j in 0..d {
 524            decomp_r[j].as_mut()[ri] = el_decomposed[j];
 525        }
 526    }
 527}
 528
 529/// Sends RLWE_{s}(X) -> RLWE_{s}(X^k) where k is some galois element
 530///
 531/// - scratch_matrix: must have dimension at-least d+2 x ring_size. d rows to
 532///   store decomposed polynomials and 2 for rlwe
 533pub(crate) fn galois_auto<
 534    MT: Matrix + IsTrivial + MatrixMut,
 535    Mmut: MatrixMut<MatElement = MT::MatElement>,
 536    ModOp: ArithmeticOps<Element = MT::MatElement> + VectorOps<Element = MT::MatElement>,
 537    NttOp: Ntt<Element = MT::MatElement>,
 538    D: Decomposer<Element = MT::MatElement>,
 539>(
 540    rlwe_in: &mut MT,
 541    ksk: &Mmut,
 542    scratch_matrix: &mut Mmut,
 543    auto_map_index: &[usize],
 544    auto_map_sign: &[bool],
 545    mod_op: &ModOp,
 546    ntt_op: &NttOp,
 547    decomposer: &D,
 548) where
 549    <Mmut as Matrix>::R: RowMut,
 550    <MT as Matrix>::R: RowMut,
 551    MT::MatElement: Copy + Zero,
 552{
 553    let d = decomposer.decomposition_count();
 554    let ring_size = rlwe_in.dimension().1;
 555    assert!(rlwe_in.dimension().0 == 2);
 556    assert!(scratch_matrix.fits(d + 2, ring_size));
 557
 558    let (scratch_matrix_d_ring, tmp_rlwe_out) = scratch_matrix.split_at_row_mut(d);
 559
 560    if !rlwe_in.is_trivial() {
 561        tmp_rlwe_out.iter_mut().for_each(|r| {
 562            r.as_mut().fill(Mmut::MatElement::zero());
 563        });
 564
 565        // send a(X) -> a(X^k) and decompose a(X^k)
 566        izip!(
 567            rlwe_in.get_row(0),
 568            auto_map_index.iter(),
 569            auto_map_sign.iter()
 570        )
 571        .for_each(|(el_in, to_index, sign)| {
 572            let el_out = if !*sign { mod_op.neg(el_in) } else { *el_in };
 573
 574            let el_out_decomposed = decomposer.decompose_to_vec(&el_out);
 575            for j in 0..d {
 576                scratch_matrix_d_ring[j].as_mut()[*to_index] = el_out_decomposed[j];
 577            }
 578        });
 579
 580        // transform decomposed a(X^k) to evaluation domain
 581        scratch_matrix_d_ring.iter_mut().for_each(|r| {
 582            ntt_op.forward(r.as_mut());
 583        });
 584
 585        // RLWE(m^k) = a', b'; RLWE(m) = a, b
 586        // key switch: (a * RLWE'(s(X^k)))
 587        let (ksk_a, ksk_b) = ksk.split_at_row(d);
 588        // a' = decomp<a> * RLWE'_A(s(X^k))
 589        routine(
 590            tmp_rlwe_out[0].as_mut(),
 591            scratch_matrix_d_ring,
 592            ksk_a,
 593            mod_op,
 594        );
 595
 596        // b' += decomp<a(X^k)> * RLWE'_B(s(X^k))
 597        routine(
 598            tmp_rlwe_out[1].as_mut(),
 599            scratch_matrix_d_ring,
 600            ksk_b,
 601            mod_op,
 602        );
 603
 604        // transform RLWE(m^k) to coefficient domain
 605        tmp_rlwe_out
 606            .iter_mut()
 607            .for_each(|r| ntt_op.backward(r.as_mut()));
 608
 609        // send b(X) -> b(X^k) and then b'(X) += b(X^k)
 610        izip!(
 611            rlwe_in.get_row(1),
 612            auto_map_index.iter(),
 613            auto_map_sign.iter()
 614        )
 615        .for_each(|(el_in, to_index, sign)| {
 616            let row = tmp_rlwe_out[1].as_mut();
 617            if !*sign {
 618                row[*to_index] = mod_op.sub(&row[*to_index], el_in);
 619            } else {
 620                row[*to_index] = mod_op.add(&row[*to_index], el_in);
 621            }
 622        });
 623
 624        // copy over A; Leave B for later
 625        rlwe_in
 626            .get_row_mut(0)
 627            .copy_from_slice(tmp_rlwe_out[0].as_ref());
 628    } else {
 629        // RLWE is trivial, a(X) is 0.
 630        // send b(X) -> b(X^k)
 631        izip!(
 632            rlwe_in.get_row(1),
 633            auto_map_index.iter(),
 634            auto_map_sign.iter()
 635        )
 636        .for_each(|(el_in, to_index, sign)| {
 637            if !*sign {
 638                tmp_rlwe_out[1].as_mut()[*to_index] = mod_op.neg(el_in);
 639            } else {
 640                tmp_rlwe_out[1].as_mut()[*to_index] = *el_in;
 641            }
 642        });
 643    }
 644
 645    // Copy over B
 646    rlwe_in
 647        .get_row_mut(1)
 648        .copy_from_slice(tmp_rlwe_out[1].as_ref());
 649}
 650
 651/// Returns RLWE(m0m1) = RLWE(m0) x RGSW(m1). Mutates rlwe_in inplace to equal
 652/// RLWE(m0m1)
 653///
 654/// - rlwe_in: is RLWE(m0) with polynomials in coefficient domain
 655/// - rgsw_in: is RGSW(m1) with polynomials in evaluation domain
 656/// - scratch_matrix_d_ring: is a matrix of dimension (d_rgsw, ring_size) used
 657///   as scratch space to store decomposed Ring elements temporarily
 658pub(crate) fn less1_rlwe_by_rgsw<
 659    Mmut: MatrixMut,
 660    MT: Matrix<MatElement = Mmut::MatElement> + MatrixMut<MatElement = Mmut::MatElement> + IsTrivial,
 661    D: Decomposer<Element = Mmut::MatElement>,
 662    ModOp: VectorOps<Element = Mmut::MatElement>,
 663    NttOp: Ntt<Element = Mmut::MatElement>,
 664>(
 665    rlwe_in: &mut MT,
 666    rgsw_in: &Mmut,
 667    scratch_matrix_dplus2_ring: &mut Mmut,
 668    decomposer: &D,
 669    ntt_op: &NttOp,
 670    mod_op: &ModOp,
 671    skip0: usize,
 672    skip1: usize,
 673) where
 674    Mmut::MatElement: Copy + Zero,
 675    <Mmut as Matrix>::R: RowMut,
 676    <MT as Matrix>::R: RowMut,
 677{
 678    let d_rgsw = decomposer.decomposition_count();
 679    assert!(scratch_matrix_dplus2_ring.dimension() == (d_rgsw + 2, rlwe_in.dimension().1));
 680    assert!(rgsw_in.dimension() == (d_rgsw * 4, rlwe_in.dimension().1));
 681
 682    // decomposed RLWE x RGSW
 683    let (rlwe_dash_nsm, rlwe_dash_m) = rgsw_in.split_at_row(d_rgsw * 2);
 684    let (scratch_matrix_d_ring, scratch_rlwe_out) =
 685        scratch_matrix_dplus2_ring.split_at_row_mut(d_rgsw);
 686    scratch_rlwe_out[0].as_mut().fill(Mmut::MatElement::zero());
 687    scratch_rlwe_out[1].as_mut().fill(Mmut::MatElement::zero());
 688    // RLWE_in = a_in, b_in; RLWE_out = a_out, b_out
 689    if !rlwe_in.is_trivial() {
 690        // a_in = 0 when RLWE_in is trivial RLWE ciphertext
 691        // decomp<a_in>
 692        decompose_r(rlwe_in.get_row_slice(0), scratch_matrix_d_ring, decomposer);
 693        scratch_matrix_d_ring
 694            .iter_mut()
 695            .for_each(|r| ntt_op.forward(r.as_mut()));
 696        // a_out += decomp<a_in> \cdot RLWE_A'(-sm)
 697        routine(
 698            scratch_rlwe_out[0].as_mut(),
 699            scratch_matrix_d_ring[skip0..].as_ref(),
 700            &rlwe_dash_nsm[skip0..d_rgsw],
 701            mod_op,
 702        );
 703        // b_out += decomp<a_in> \cdot RLWE_B'(-sm)
 704        routine(
 705            scratch_rlwe_out[1].as_mut(),
 706            scratch_matrix_d_ring[skip0..].as_ref(),
 707            &rlwe_dash_nsm[d_rgsw + skip0..],
 708            mod_op,
 709        );
 710    }
 711    // decomp<b_in>
 712    decompose_r(rlwe_in.get_row_slice(1), scratch_matrix_d_ring, decomposer);
 713    scratch_matrix_d_ring
 714        .iter_mut()
 715        .for_each(|r| ntt_op.forward(r.as_mut()));
 716    // a_out += decomp<b_in> \cdot RLWE_A'(m)
 717    routine(
 718        scratch_rlwe_out[0].as_mut(),
 719        scratch_matrix_d_ring[skip1..].as_ref(),
 720        &rlwe_dash_m[skip1..d_rgsw],
 721        mod_op,
 722    );
 723    // b_out += decomp<b_in> \cdot RLWE_B'(m)
 724    routine(
 725        scratch_rlwe_out[1].as_mut(),
 726        scratch_matrix_d_ring[skip1..].as_ref(),
 727        &rlwe_dash_m[d_rgsw + skip1..],
 728        mod_op,
 729    );
 730
 731    // transform rlwe_out to coefficient domain
 732    scratch_rlwe_out
 733        .iter_mut()
 734        .for_each(|r| ntt_op.backward(r.as_mut()));
 735
 736    rlwe_in
 737        .get_row_mut(0)
 738        .copy_from_slice(scratch_rlwe_out[0].as_mut());
 739    rlwe_in
 740        .get_row_mut(1)
 741        .copy_from_slice(scratch_rlwe_out[1].as_mut());
 742    rlwe_in.set_not_trivial();
 743}
 744
 745/// Returns RLWE(m0m1) = RLWE(m0) x RGSW(m1). Mutates rlwe_in inplace to equal
 746/// RLWE(m0m1)
 747///
 748/// - rlwe_in: is RLWE(m0) with polynomials in coefficient domain
 749/// - rgsw_in: is RGSW(m1) with polynomials in evaluation domain
 750/// - scratch_matrix_d_ring: is a matrix with atleast max(d_a, d_b) rows and
 751///   ring_size columns. It's used to store decomposed polynomials and out RLWE
 752///   temoporarily
 753pub(crate) fn rlwe_by_rgsw<
 754    Mmut: MatrixMut,
 755    MT: Matrix<MatElement = Mmut::MatElement> + MatrixMut<MatElement = Mmut::MatElement> + IsTrivial,
 756    D: RlweDecomposer<Element = Mmut::MatElement>,
 757    ModOp: VectorOps<Element = Mmut::MatElement>,
 758    NttOp: Ntt<Element = Mmut::MatElement>,
 759>(
 760    rlwe_in: &mut MT,
 761    rgsw_in: &Mmut,
 762    scratch_matrix: &mut Mmut,
 763    decomposer: &D,
 764    ntt_op: &NttOp,
 765    mod_op: &ModOp,
 766) where
 767    Mmut::MatElement: Copy + Zero,
 768    <Mmut as Matrix>::R: RowMut,
 769    <MT as Matrix>::R: RowMut,
 770{
 771    let decomposer_a = decomposer.a();
 772    let decomposer_b = decomposer.b();
 773    let d_a = decomposer_a.decomposition_count();
 774    let d_b = decomposer_b.decomposition_count();
 775    let max_d = std::cmp::max(d_a, d_b);
 776    assert!(scratch_matrix.fits(max_d + 2, rlwe_in.dimension().1));
 777    assert!(rgsw_in.dimension() == (d_a * 2 + d_b * 2, rlwe_in.dimension().1));
 778
 779    // decomposed RLWE x RGSW
 780    let (rlwe_dash_nsm, rlwe_dash_m) = rgsw_in.split_at_row(d_a * 2);
 781    let (scratch_matrix_d_ring, scratch_rlwe_out) = scratch_matrix.split_at_row_mut(max_d);
 782    scratch_rlwe_out[0].as_mut().fill(Mmut::MatElement::zero());
 783    scratch_rlwe_out[1].as_mut().fill(Mmut::MatElement::zero());
 784    // RLWE_in = a_in, b_in; RLWE_out = a_out, b_out
 785    if !rlwe_in.is_trivial() {
 786        // a_in = 0 when RLWE_in is trivial RLWE ciphertext
 787        // decomp<a_in>
 788        decompose_r(
 789            rlwe_in.get_row_slice(0),
 790            scratch_matrix_d_ring,
 791            decomposer_a,
 792        );
 793        scratch_matrix_d_ring
 794            .iter_mut()
 795            .take(d_a)
 796            .for_each(|r| ntt_op.forward(r.as_mut()));
 797        // a_out += decomp<a_in> \cdot RLWE_A'(-sm)
 798        routine(
 799            scratch_rlwe_out[0].as_mut(),
 800            scratch_matrix_d_ring.as_ref(),
 801            &rlwe_dash_nsm[..d_a],
 802            mod_op,
 803        );
 804        // b_out += decomp<a_in> \cdot RLWE_B'(-sm)
 805        routine(
 806            scratch_rlwe_out[1].as_mut(),
 807            scratch_matrix_d_ring.as_ref(),
 808            &rlwe_dash_nsm[d_a..],
 809            mod_op,
 810        );
 811    }
 812    // decomp<b_in>
 813    decompose_r(
 814        rlwe_in.get_row_slice(1),
 815        scratch_matrix_d_ring,
 816        decomposer_b,
 817    );
 818    scratch_matrix_d_ring
 819        .iter_mut()
 820        .take(d_b)
 821        .for_each(|r| ntt_op.forward(r.as_mut()));
 822    // a_out += decomp<b_in> \cdot RLWE_A'(m)
 823    routine(
 824        scratch_rlwe_out[0].as_mut(),
 825        scratch_matrix_d_ring.as_ref(),
 826        &rlwe_dash_m[..d_b],
 827        mod_op,
 828    );
 829    // b_out += decomp<b_in> \cdot RLWE_B'(m)
 830    routine(
 831        scratch_rlwe_out[1].as_mut(),
 832        scratch_matrix_d_ring.as_ref(),
 833        &rlwe_dash_m[d_b..],
 834        mod_op,
 835    );
 836
 837    // transform rlwe_out to coefficient domain
 838    scratch_rlwe_out
 839        .iter_mut()
 840        .for_each(|r| ntt_op.backward(r.as_mut()));
 841
 842    rlwe_in
 843        .get_row_mut(0)
 844        .copy_from_slice(scratch_rlwe_out[0].as_mut());
 845    rlwe_in
 846        .get_row_mut(1)
 847        .copy_from_slice(scratch_rlwe_out[1].as_mut());
 848    rlwe_in.set_not_trivial();
 849}
 850
 851/// Inplace mutates rlwe_0 to equal RGSW(m0m1) = RGSW(m0)xRGSW(m1)
 852/// in evaluation domain
 853///
 854/// Warning -
 855/// Pass a fresh RGSW ciphertext as the second operand, i.e. as `rgsw_1`.
 856/// This is to assure minimal error growth in the resulting RGSW ciphertext.
 857/// RGSWxRGSW boils down to d_rgsw*2 RLWExRGSW multiplications. Hence, the noise
 858/// growth in resulting ciphertext depends on the norm of second RGSW
 859/// ciphertext, not the first. This is useful in cases where one is accumulating
 860/// multiple RGSW ciphertexts into 1. In which case, pass the accumulating RGSW
 861/// ciphertext as rlwe_0 (the one with higher noise) and subsequent RGSW
 862/// ciphertexts, with lower noise, to be accumulated as second
 863/// operand.
 864///
 865/// - rgsw_0: RGSW(m0)
 866/// - rgsw_1_eval: RGSW(m1) in Evaluation domain
 867/// - scratch_matrix_d_plus_rgsw_by_ring: scratch space matrix with rows
 868///   (max(d_a, d_b) + d_a*2+d_b*2) and columns ring_size
 869pub(crate) fn rgsw_by_rgsw_inplace<
 870    Mmut: MatrixMut,
 871    D: RlweDecomposer<Element = Mmut::MatElement>,
 872    ModOp: VectorOps<Element = Mmut::MatElement>,
 873    NttOp: Ntt<Element = Mmut::MatElement>,
 874>(
 875    rgsw_0: &mut Mmut,
 876    rgsw_1_eval: &Mmut,
 877    decomposer: &D,
 878    scratch_matrix: &mut Mmut,
 879    ntt_op: &NttOp,
 880    mod_op: &ModOp,
 881) where
 882    <Mmut as Matrix>::R: RowMut,
 883    Mmut::MatElement: Copy + Zero,
 884{
 885    let decomposer_a = decomposer.a();
 886    let decomposer_b = decomposer.b();
 887    let d_a = decomposer_a.decomposition_count();
 888    let d_b = decomposer_b.decomposition_count();
 889    let max_d = std::cmp::max(d_a, d_b);
 890    let rgsw_rows = d_a * 2 + d_b * 2;
 891    assert!(rgsw_0.dimension().0 == rgsw_rows);
 892    let ring_size = rgsw_0.dimension().1;
 893    assert!(rgsw_1_eval.dimension() == (rgsw_rows, ring_size));
 894    assert!(scratch_matrix.fits(max_d + rgsw_rows, ring_size));
 895
 896    let (decomp_r_space, rgsw_space) = scratch_matrix.split_at_row_mut(max_d);
 897
 898    // zero rgsw_space
 899    rgsw_space
 900        .iter_mut()
 901        .for_each(|ri| ri.as_mut().fill(Mmut::MatElement::zero()));
 902    let (rlwe_dash_space_nsm, rlwe_dash_space_m) = rgsw_space.split_at_mut(d_a * 2);
 903    let (rlwe_dash_space_nsm_parta, rlwe_dash_space_nsm_partb) =
 904        rlwe_dash_space_nsm.split_at_mut(d_a);
 905    let (rlwe_dash_space_m_parta, rlwe_dash_space_m_partb) = rlwe_dash_space_m.split_at_mut(d_b);
 906
 907    let (rgsw0_nsm, rgsw0_m) = rgsw_0.split_at_row(d_a * 2);
 908    let (rgsw1_nsm, rgsw1_m) = rgsw_1_eval.split_at_row(d_a * 2);
 909
 910    // RGSW x RGSW
 911    izip!(
 912        rgsw0_nsm.iter().take(d_a).chain(rgsw0_m.iter().take(d_b)),
 913        rgsw0_nsm.iter().skip(d_a).chain(rgsw0_m.iter().skip(d_b)),
 914        rlwe_dash_space_nsm_parta
 915            .iter_mut()
 916            .chain(rlwe_dash_space_m_parta.iter_mut()),
 917        rlwe_dash_space_nsm_partb
 918            .iter_mut()
 919            .chain(rlwe_dash_space_m_partb.iter_mut()),
 920    )
 921    .for_each(|(rlwe_a, rlwe_b, rlwe_out_a, rlwe_out_b)| {
 922        // Part A
 923        decompose_r(rlwe_a.as_ref(), decomp_r_space.as_mut(), decomposer_a);
 924        decomp_r_space
 925            .iter_mut()
 926            .take(d_a)
 927            .for_each(|ri| ntt_op.forward(ri.as_mut()));
 928        routine(
 929            rlwe_out_a.as_mut(),
 930            &decomp_r_space[..d_a],
 931            &rgsw1_nsm[..d_a],
 932            mod_op,
 933        );
 934        routine(
 935            rlwe_out_b.as_mut(),
 936            &decomp_r_space[..d_a],
 937            &rgsw1_nsm[d_a..],
 938            mod_op,
 939        );
 940
 941        // Part B
 942        decompose_r(rlwe_b.as_ref(), decomp_r_space.as_mut(), decomposer_b);
 943        decomp_r_space
 944            .iter_mut()
 945            .take(d_b)
 946            .for_each(|ri| ntt_op.forward(ri.as_mut()));
 947        routine(
 948            rlwe_out_a.as_mut(),
 949            &decomp_r_space[..d_b],
 950            &rgsw1_m[..d_b],
 951            mod_op,
 952        );
 953        routine(
 954            rlwe_out_b.as_mut(),
 955            &decomp_r_space[..d_b],
 956            &rgsw1_m[d_b..],
 957            mod_op,
 958        );
 959    });
 960
 961    // copy over RGSW(m0m1) into RGSW(m0)
 962    izip!(rgsw_0.iter_rows_mut(), rgsw_space.iter())
 963        .for_each(|(to_ri, from_ri)| to_ri.as_mut().copy_from_slice(from_ri.as_ref()));
 964
 965    // send back to coefficient domain
 966    rgsw_0
 967        .iter_rows_mut()
 968        .for_each(|ri| ntt_op.backward(ri.as_mut()));
 969}
 970
 971/// Encrypts message m as a RGSW ciphertext.
 972///
 973/// - m_eval: is `m` is evaluation domain
 974/// - out_rgsw: RGSW(m) is stored as single matrix of dimension (d_rgsw * 3,
 975///   ring_size). The matrix has the following structure [RLWE'_A(-sm) ||
 976///   RLWE'_B(-sm) || RLWE'_B(m)]^T and RLWE'_A(m) is generated via seed (where
 977///   p_rng is assumed to be seeded with seed)
 978pub(crate) fn secret_key_encrypt_rgsw<
 979    Mmut: MatrixMut + MatrixEntity,
 980    S,
 981    R: RandomFillGaussianInModulus<[Mmut::MatElement], ModOp::M>
 982        + RandomFillUniformInModulus<[Mmut::MatElement], ModOp::M>,
 983    PR: RandomFillUniformInModulus<[Mmut::MatElement], ModOp::M>,
 984    ModOp: VectorOps<Element = Mmut::MatElement> + GetModulus<Element = Mmut::MatElement>,
 985    NttOp: Ntt<Element = Mmut::MatElement>,
 986>(
 987    out_rgsw: &mut Mmut,
 988    m: &[Mmut::MatElement],
 989    gadget_a: &[Mmut::MatElement],
 990    gadget_b: &[Mmut::MatElement],
 991    s: &[S],
 992    mod_op: &ModOp,
 993    ntt_op: &NttOp,
 994    p_rng: &mut PR,
 995    rng: &mut R,
 996) where
 997    <Mmut as Matrix>::R: RowMut + RowEntity + TryConvertFrom1<[S], ModOp::M> + Debug,
 998    Mmut::MatElement: Copy + Debug,
 999{
1000    let d_a = gadget_a.len();
1001    let d_b = gadget_b.len();
1002    let q = mod_op.modulus();
1003    let ring_size = s.len();
1004    assert!(out_rgsw.dimension() == (d_a * 2 + d_b, ring_size));
1005    assert!(m.as_ref().len() == ring_size);
1006
1007    // RLWE(-sm), RLWE(m)
1008    let (rlwe_dash_nsm, b_rlwe_dash_m) = out_rgsw.split_at_row_mut(d_a * 2);
1009
1010    let mut s_eval = Mmut::R::try_convert_from(s, &q);
1011    ntt_op.forward(s_eval.as_mut());
1012
1013    let mut scratch_space = Mmut::R::zeros(ring_size);
1014
1015    // RLWE'(-sm)
1016    let (a_rlwe_dash_nsm, b_rlwe_dash_nsm) = rlwe_dash_nsm.split_at_mut(d_a);
1017    izip!(
1018        a_rlwe_dash_nsm.iter_mut(),
1019        b_rlwe_dash_nsm.iter_mut(),
1020        gadget_a.iter()
1021    )
1022    .for_each(|(ai, bi, beta_i)| {
1023        // Sample a_i
1024        RandomFillUniformInModulus::random_fill(rng, &q, ai.as_mut());
1025
1026        // a_i * s
1027        scratch_space.as_mut().copy_from_slice(ai.as_ref());
1028        ntt_op.forward(scratch_space.as_mut());
1029        mod_op.elwise_mul_mut(scratch_space.as_mut(), s_eval.as_ref());
1030        ntt_op.backward(scratch_space.as_mut());
1031
1032        // b_i = e_i + a_i * s
1033        RandomFillGaussianInModulus::random_fill(rng, &q, bi.as_mut());
1034        mod_op.elwise_add_mut(bi.as_mut(), scratch_space.as_ref());
1035
1036        // a_i + \beta_i * m
1037        mod_op.elwise_scalar_mul(scratch_space.as_mut(), m.as_ref(), beta_i);
1038        mod_op.elwise_add_mut(ai.as_mut(), scratch_space.as_ref());
1039    });
1040
1041    // RLWE(m)
1042    let mut a_rlwe_dash_m = {
1043        // polynomials of part A of RLWE'(m) are sampled from seed
1044        let mut a = Mmut::zeros(d_b, ring_size);
1045        a.iter_rows_mut()
1046            .for_each(|ai| RandomFillUniformInModulus::random_fill(p_rng, &q, ai.as_mut()));
1047        a
1048    };
1049
1050    izip!(
1051        a_rlwe_dash_m.iter_rows_mut(),
1052        b_rlwe_dash_m.iter_mut(),
1053        gadget_b.iter()
1054    )
1055    .for_each(|(ai, bi, beta_i)| {
1056        // ai * s
1057        ntt_op.forward(ai.as_mut());
1058        mod_op.elwise_mul_mut(ai.as_mut(), s_eval.as_ref());
1059        ntt_op.backward(ai.as_mut());
1060
1061        // beta_i * m
1062        mod_op.elwise_scalar_mul(scratch_space.as_mut(), m.as_ref(), beta_i);
1063
1064        // Sample e_i
1065        RandomFillGaussianInModulus::random_fill(rng, &q, bi.as_mut());
1066        // e_i + beta_i * m + ai*s
1067        mod_op.elwise_add_mut(bi.as_mut(), scratch_space.as_ref());
1068        mod_op.elwise_add_mut(bi.as_mut(), ai.as_ref());
1069    });
1070}
1071
1072pub(crate) fn public_key_encrypt_rgsw<
1073    Mmut: MatrixMut + MatrixEntity,
1074    M: Matrix<MatElement = Mmut::MatElement>,
1075    R: RandomFillGaussianInModulus<[Mmut::MatElement], ModOp::M>
1076        + RandomFill<[u8]>
1077        + RandomElementInModulus<usize, usize>,
1078    ModOp: VectorOps<Element = Mmut::MatElement> + GetModulus<Element = Mmut::MatElement>,
1079    NttOp: Ntt<Element = Mmut::MatElement>,
1080>(
1081    out_rgsw: &mut Mmut,
1082    m: &[M::MatElement],
1083    public_key: &M,
1084    gadget_a: &[Mmut::MatElement],
1085    gadget_b: &[Mmut::MatElement],
1086    mod_op: &ModOp,
1087    ntt_op: &NttOp,
1088    rng: &mut R,
1089) where
1090    <Mmut as Matrix>::R: RowMut + RowEntity + TryConvertFrom1<[i32], ModOp::M>,
1091    Mmut::MatElement: Copy,
1092{
1093    let ring_size = public_key.dimension().1;
1094    let d_a = gadget_a.len();
1095    let d_b = gadget_b.len();
1096    assert!(public_key.dimension().0 == 2);
1097    assert!(out_rgsw.dimension() == (d_a * 2 + d_b * 2, ring_size));
1098
1099    let mut pk_eval = Mmut::zeros(2, ring_size);
1100    izip!(pk_eval.iter_rows_mut(), public_key.iter_rows()).for_each(|(to_i, from_i)| {
1101        to_i.as_mut().copy_from_slice(from_i.as_ref());
1102        ntt_op.forward(to_i.as_mut());
1103    });
1104    let p0 = pk_eval.get_row_slice(0);
1105    let p1 = pk_eval.get_row_slice(1);
1106
1107    let q = mod_op.modulus();
1108
1109    // RGSW(m) = RLWE'(-sm), RLWE(m)
1110    let (rlwe_dash_nsm, rlwe_dash_m) = out_rgsw.split_at_row_mut(d_a * 2);
1111
1112    // RLWE(-sm)
1113    let (rlwe_dash_nsm_parta, rlwe_dash_nsm_partb) = rlwe_dash_nsm.split_at_mut(d_a);
1114    izip!(
1115        rlwe_dash_nsm_parta.iter_mut(),
1116        rlwe_dash_nsm_partb.iter_mut(),
1117        gadget_a.iter()
1118    )
1119    .for_each(|(ai, bi, beta_i)| {
1120        // sample ephemeral secret u_i
1121        let mut u = vec![0i32; ring_size];
1122        fill_random_ternary_secret_with_hamming_weight(u.as_mut(), ring_size >> 1, rng);
1123        let mut u_eval = Mmut::R::try_convert_from(u.as_ref(), &q);
1124        ntt_op.forward(u_eval.as_mut());
1125
1126        let mut u_eval_copy = Mmut::R::zeros(ring_size);
1127        u_eval_copy.as_mut().copy_from_slice(u_eval.as_ref());
1128
1129        // p0 * u
1130        mod_op.elwise_mul_mut(u_eval.as_mut(), p0.as_ref());
1131        // p1 * u
1132        mod_op.elwise_mul_mut(u_eval_copy.as_mut(), p1.as_ref());
1133        ntt_op.backward(u_eval.as_mut());
1134        ntt_op.backward(u_eval_copy.as_mut());
1135
1136        // sample error
1137        RandomFillGaussianInModulus::random_fill(rng, &q, ai.as_mut());
1138        RandomFillGaussianInModulus::random_fill(rng, &q, bi.as_mut());
1139
1140        // a = p0*u+e0
1141        mod_op.elwise_add_mut(ai.as_mut(), u_eval.as_ref());
1142        // b = p1*u+e1
1143        mod_op.elwise_add_mut(bi.as_mut(), u_eval_copy.as_ref());
1144
1145        // a = p0*u + e0 + \beta*m
1146        // use u_eval as scratch
1147        mod_op.elwise_scalar_mul(u_eval.as_mut(), m.as_ref(), beta_i);
1148        mod_op.elwise_add_mut(ai.as_mut(), u_eval.as_ref());
1149    });
1150
1151    // RLWE(m)
1152    let (rlwe_dash_m_parta, rlwe_dash_m_partb) = rlwe_dash_m.split_at_mut(d_b);
1153    izip!(
1154        rlwe_dash_m_parta.iter_mut(),
1155        rlwe_dash_m_partb.iter_mut(),
1156        gadget_b.iter()
1157    )
1158    .for_each(|(ai, bi, beta_i)| {
1159        // sample ephemeral secret u_i
1160        let mut u = vec![0i32; ring_size];
1161        fill_random_ternary_secret_with_hamming_weight(u.as_mut(), ring_size >> 1, rng);
1162        let mut u_eval = Mmut::R::try_convert_from(u.as_ref(), &q);
1163        ntt_op.forward(u_eval.as_mut());
1164
1165        let mut u_eval_copy = Mmut::R::zeros(ring_size);
1166        u_eval_copy.as_mut().copy_from_slice(u_eval.as_ref());
1167
1168        // p0 * u
1169        mod_op.elwise_mul_mut(u_eval.as_mut(), p0.as_ref());
1170        // p1 * u
1171        mod_op.elwise_mul_mut(u_eval_copy.as_mut(), p1.as_ref());
1172        ntt_op.backward(u_eval.as_mut());
1173        ntt_op.backward(u_eval_copy.as_mut());
1174
1175        // sample error
1176        RandomFillGaussianInModulus::random_fill(rng, &q, ai.as_mut());
1177        RandomFillGaussianInModulus::random_fill(rng, &q, bi.as_mut());
1178
1179        // a = p0*u+e0
1180        mod_op.elwise_add_mut(ai.as_mut(), u_eval.as_ref());
1181        // b = p1*u+e1
1182        mod_op.elwise_add_mut(bi.as_mut(), u_eval_copy.as_ref());
1183
1184        // b = p1*u + e0 + \beta*m
1185        // use u_eval as scratch
1186        mod_op.elwise_scalar_mul(u_eval.as_mut(), m.as_ref(), beta_i);
1187        mod_op.elwise_add_mut(bi.as_mut(), u_eval.as_ref());
1188    });
1189}
1190
1191/// Generates RLWE Key switching key to key switch ciphertext RLWE_{from_s}(m)
1192/// to RLWE_{to_s}(m).
1193///
1194/// Key switching equals
1195///     \sum decompose(c_1)_i * RLWE_{to_s}(\beta^i -from_s)
1196/// Hence, key switchin key equals RLWE'(-from_s) = RLWE(-from_s), RLWE(beta^1
1197/// -from_s), ..., RLWE(beta^{d-1} -from_s).
1198///
1199/// - ksk_out: Output Key switching key. Key switching key stores only part B
1200///   polynomials of ksk RLWE ciphertexts (i.e. RLWE'_B(-from_s)) in coefficient
1201///   domain
1202/// - neg_from_s: Negative of secret polynomial to key switch from
1203/// - to_s: secret polynomial to key switch to.
1204pub(crate) fn rlwe_ksk_gen<
1205    Mmut: MatrixMut + MatrixEntity,
1206    ModOp: ArithmeticOps<Element = Mmut::MatElement>
1207        + VectorOps<Element = Mmut::MatElement>
1208        + GetModulus<Element = Mmut::MatElement>,
1209    NttOp: Ntt<Element = Mmut::MatElement>,
1210    R: RandomFillGaussianInModulus<[Mmut::MatElement], ModOp::M>,
1211    PR: RandomFillUniformInModulus<[Mmut::MatElement], ModOp::M>,
1212>(
1213    ksk_out: &mut Mmut,
1214    neg_from_s: Mmut::R,
1215    mut to_s: Mmut::R,
1216    gadget_vector: &[Mmut::MatElement],
1217    mod_op: &ModOp,
1218    ntt_op: &NttOp,
1219    p_rng: &mut PR,
1220    rng: &mut R,
1221) where
1222    <Mmut as Matrix>::R: RowMut,
1223{
1224    let ring_size = neg_from_s.as_ref().len();
1225    let d = gadget_vector.len();
1226    assert!(ksk_out.dimension() == (d, ring_size));
1227
1228    let q = mod_op.modulus();
1229
1230    ntt_op.forward(to_s.as_mut());
1231
1232    // RLWE'_{to_s}(-from_s)
1233    let mut part_a = {
1234        let mut a = Mmut::zeros(d, ring_size);
1235        a.iter_rows_mut()
1236            .for_each(|ai| RandomFillUniformInModulus::random_fill(p_rng, q, ai.as_mut()));
1237        a
1238    };
1239    izip!(
1240        part_a.iter_rows_mut(),
1241        ksk_out.iter_rows_mut(),
1242        gadget_vector.iter(),
1243    )
1244    .for_each(|(ai, bi, beta_i)| {
1245        // si * ai
1246        ntt_op.forward(ai.as_mut());
1247        mod_op.elwise_mul_mut(ai.as_mut(), to_s.as_ref());
1248        ntt_op.backward(ai.as_mut());
1249
1250        // ei + to_s*ai
1251        RandomFillGaussianInModulus::random_fill(rng, &q, bi.as_mut());
1252        mod_op.elwise_add_mut(bi.as_mut(), ai.as_ref());
1253
1254        // beta_i * -from_s
1255        // use ai as scratch space
1256        mod_op.elwise_scalar_mul(ai.as_mut(), neg_from_s.as_ref(), beta_i);
1257
1258        // bi = ei + to_s*ai + beta_i*-from_s
1259        mod_op.elwise_add_mut(bi.as_mut(), ai.as_ref());
1260    });
1261}
1262
1263pub(crate) fn galois_key_gen<
1264    Mmut: MatrixMut + MatrixEntity,
1265    ModOp: ArithmeticOps<Element = Mmut::MatElement>
1266        + VectorOps<Element = Mmut::MatElement>
1267        + GetModulus<Element = Mmut::MatElement>,
1268    NttOp: Ntt<Element = Mmut::MatElement>,
1269    S,
1270    R: RandomFillGaussianInModulus<[Mmut::MatElement], ModOp::M>,
1271    PR: RandomFillUniformInModulus<[Mmut::MatElement], ModOp::M>,
1272>(
1273    ksk_out: &mut Mmut,
1274    s: &[S],
1275    auto_k: isize,
1276    gadget_vector: &[Mmut::MatElement],
1277    mod_op: &ModOp,
1278    ntt_op: &NttOp,
1279    p_rng: &mut PR,
1280    rng: &mut R,
1281) where
1282    <Mmut as Matrix>::R: RowMut,
1283    Mmut::R: TryConvertFrom1<[S], ModOp::M> + RowEntity,
1284    Mmut::MatElement: Copy + Sub<Output = Mmut::MatElement>,
1285{
1286    let ring_size = s.len();
1287    let (auto_map_index, auto_map_sign) = generate_auto_map(ring_size, auto_k);
1288
1289    let q = mod_op.modulus();
1290
1291    // s(X) -> -s(X^k)
1292    let s = Mmut::R::try_convert_from(s, q);
1293    let mut neg_s_auto = Mmut::R::zeros(s.as_ref().len());
1294    izip!(s.as_ref(), auto_map_index.iter(), auto_map_sign.iter()).for_each(
1295        |(el, to_index, sign)| {
1296            // if sign is +ve (true), then negate because we need -s(X) (i.e. do the
1297            // opposite than the usual case)
1298            if *sign {
1299                neg_s_auto.as_mut()[*to_index] = mod_op.neg(el);
1300            } else {
1301                neg_s_auto.as_mut()[*to_index] = *el;
1302            }
1303        },
1304    );
1305
1306    // Ksk from -s(X^k) to s(X)
1307    rlwe_ksk_gen(
1308        ksk_out,
1309        neg_s_auto,
1310        s,
1311        gadget_vector,
1312        mod_op,
1313        ntt_op,
1314        p_rng,
1315        rng,
1316    );
1317}
1318
1319/// Encrypt polynomial m(X) as RLWE ciphertext.
1320///
1321/// - rlwe_out: returned RLWE ciphertext RLWE(m) in coefficient domain. RLWE
1322///   ciphertext is a matirx with first row consiting polynomial `a` and the
1323///   second rows consting polynomial `b`
1324pub(crate) fn secret_key_encrypt_rlwe<
1325    Ro: Row + RowMut + RowEntity,
1326    ModOp: VectorOps<Element = Ro::Element> + GetModulus<Element = Ro::Element>,
1327    NttOp: Ntt<Element = Ro::Element>,
1328    S,
1329    R: RandomFillGaussianInModulus<[Ro::Element], ModOp::M>,
1330    PR: RandomFillUniformInModulus<[Ro::Element], ModOp::M>,
1331>(
1332    m: &[Ro::Element],
1333    b_rlwe_out: &mut Ro,
1334    s: &[S],
1335    mod_op: &ModOp,
1336    ntt_op: &NttOp,
1337    p_rng: &mut PR,
1338    rng: &mut R,
1339) where
1340    Ro: TryConvertFrom1<[S], ModOp::M> + Debug,
1341{
1342    let ring_size = s.len();
1343    assert!(m.as_ref().len() == ring_size);
1344    assert!(b_rlwe_out.as_ref().len() == ring_size);
1345
1346    let q = mod_op.modulus();
1347
1348    // sample a
1349    let mut a = {
1350        let mut a = Ro::zeros(ring_size);
1351        RandomFillUniformInModulus::random_fill(p_rng, q, a.as_mut());
1352        a
1353    };
1354
1355    // s * a
1356    let mut sa = Ro::try_convert_from(s, q);
1357    ntt_op.forward(sa.as_mut());
1358    ntt_op.forward(a.as_mut());
1359    mod_op.elwise_mul_mut(sa.as_mut(), a.as_ref());
1360    ntt_op.backward(sa.as_mut());
1361
1362    // sample e
1363    RandomFillGaussianInModulus::random_fill(rng, q, b_rlwe_out.as_mut());
1364    mod_op.elwise_add_mut(b_rlwe_out.as_mut(), m.as_ref());
1365    mod_op.elwise_add_mut(b_rlwe_out.as_mut(), sa.as_ref());
1366}
1367
1368pub(crate) fn public_key_encrypt_rlwe<
1369    M: Matrix,
1370    Mmut: MatrixMut<MatElement = M::MatElement>,
1371    ModOp: VectorOps<Element = M::MatElement> + GetModulus<Element = M::MatElement>,
1372    NttOp: Ntt<Element = M::MatElement>,
1373    S,
1374    R: RandomFillGaussianInModulus<[M::MatElement], ModOp::M>
1375        + RandomFillUniformInModulus<[M::MatElement], ModOp::M>
1376        + RandomFill<[u8]>
1377        + RandomElementInModulus<usize, usize>,
1378>(
1379    rlwe_out: &mut Mmut,
1380    pk: &M,
1381    m: &[M::MatElement],
1382    mod_op: &ModOp,
1383    ntt_op: &NttOp,
1384    rng: &mut R,
1385) where
1386    <Mmut as Matrix>::R: RowMut + TryConvertFrom1<[S], ModOp::M> + RowEntity,
1387    M::MatElement: Copy,
1388    S: Zero + Signed + Copy,
1389{
1390    let ring_size = m.len();
1391    assert!(rlwe_out.dimension() == (2, ring_size));
1392
1393    let q = mod_op.modulus();
1394
1395    let mut u = vec![S::zero(); ring_size];
1396    fill_random_ternary_secret_with_hamming_weight(u.as_mut(), ring_size >> 1, rng);
1397    let mut u = Mmut::R::try_convert_from(&u, q);
1398    ntt_op.forward(u.as_mut());
1399
1400    let mut ua = Mmut::R::zeros(ring_size);
1401    ua.as_mut().copy_from_slice(pk.get_row_slice(0));
1402    let mut ub = Mmut::R::zeros(ring_size);
1403    ub.as_mut().copy_from_slice(pk.get_row_slice(1));
1404
1405    // a*u
1406    ntt_op.forward(ua.as_mut());
1407    mod_op.elwise_mul_mut(ua.as_mut(), u.as_ref());
1408    ntt_op.backward(ua.as_mut());
1409
1410    // b*u
1411    ntt_op.forward(ub.as_mut());
1412    mod_op.elwise_mul_mut(ub.as_mut(), u.as_ref());
1413    ntt_op.backward(ub.as_mut());
1414
1415    // sample error
1416    rlwe_out.iter_rows_mut().for_each(|ri| {
1417        RandomFillGaussianInModulus::random_fill(rng, &q, ri.as_mut());
1418    });
1419
1420    // a*u + e0
1421    mod_op.elwise_add_mut(rlwe_out.get_row_mut(0), ua.as_ref());
1422    // b*u + e1
1423    mod_op.elwise_add_mut(rlwe_out.get_row_mut(1), ub.as_ref());
1424
1425    // b*u + e1 + m
1426    mod_op.elwise_add_mut(rlwe_out.get_row_mut(1), m);
1427}
1428
1429/// Generates RLWE public key
1430pub(crate) fn gen_rlwe_public_key<
1431    Ro: RowMut + RowEntity,
1432    S,
1433    ModOp: VectorOps<Element = Ro::Element> + GetModulus<Element = Ro::Element>,
1434    NttOp: Ntt<Element = Ro::Element>,
1435    PRng: RandomFillUniformInModulus<[Ro::Element], ModOp::M>,
1436    Rng: RandomFillGaussianInModulus<[Ro::Element], ModOp::M>,
1437>(
1438    part_b_out: &mut Ro,
1439    s: &[S],
1440    ntt_op: &NttOp,
1441    mod_op: &ModOp,
1442    p_rng: &mut PRng,
1443    rng: &mut Rng,
1444) where
1445    Ro: TryConvertFrom1<[S], ModOp::M>,
1446{
1447    let ring_size = s.len();
1448    assert!(part_b_out.as_ref().len() == ring_size);
1449
1450    let q = mod_op.modulus();
1451
1452    // sample a
1453    let mut a = {
1454        let mut tmp = Ro::zeros(ring_size);
1455        RandomFillUniformInModulus::random_fill(p_rng, &q, tmp.as_mut());
1456        tmp
1457    };
1458    ntt_op.forward(a.as_mut());
1459
1460    // s*a
1461    let mut sa = Ro::try_convert_from(s, &q);
1462    ntt_op.forward(sa.as_mut());
1463    mod_op.elwise_mul_mut(sa.as_mut(), a.as_ref());
1464    ntt_op.backward(sa.as_mut());
1465
1466    // s*a + e
1467    RandomFillGaussianInModulus::random_fill(rng, &q, part_b_out.as_mut());
1468    mod_op.elwise_add_mut(part_b_out.as_mut(), sa.as_ref());
1469}
1470
1471/// Decrypts degree 1 RLWE ciphertext RLWE(m) and returns m
1472///
1473/// - rlwe_ct: input degree 1 ciphertext RLWE(m).
1474pub(crate) fn decrypt_rlwe<
1475    R: RowMut,
1476    M: Matrix<MatElement = R::Element>,
1477    ModOp: VectorOps<Element = R::Element> + GetModulus<Element = R::Element>,
1478    NttOp: Ntt<Element = R::Element>,
1479    S,
1480>(
1481    rlwe_ct: &M,
1482    s: &[S],
1483    m_out: &mut R,
1484    ntt_op: &NttOp,
1485    mod_op: &ModOp,
1486) where
1487    R: TryConvertFrom1<[S], ModOp::M>,
1488    R::Element: Copy,
1489{
1490    let ring_size = s.len();
1491    assert!(rlwe_ct.dimension() == (2, ring_size));
1492    assert!(m_out.as_ref().len() == ring_size);
1493
1494    // transform a to evluation form
1495    m_out.as_mut().copy_from_slice(rlwe_ct.get_row_slice(0));
1496    ntt_op.forward(m_out.as_mut());
1497
1498    // -s*a
1499    let mut s = R::try_convert_from(&s, mod_op.modulus());
1500    ntt_op.forward(s.as_mut());
1501    mod_op.elwise_mul_mut(m_out.as_mut(), s.as_ref());
1502    mod_op.elwise_neg_mut(m_out.as_mut());
1503    ntt_op.backward(m_out.as_mut());
1504
1505    // m+e = b - s*a
1506    mod_op.elwise_add_mut(m_out.as_mut(), rlwe_ct.get_row_slice(1));
1507}
1508
1509// Measures noise in degree 1 RLWE ciphertext against encoded ideal message
1510// encoded_m
1511pub(crate) fn measure_noise<
1512    Mmut: MatrixMut + Matrix,
1513    ModOp: VectorOps<Element = Mmut::MatElement> + GetModulus<Element = Mmut::MatElement>,
1514    NttOp: Ntt<Element = Mmut::MatElement>,
1515    S,
1516>(
1517    rlwe_ct: &Mmut,
1518    encoded_m_ideal: &Mmut::R,
1519    ntt_op: &NttOp,
1520    mod_op: &ModOp,
1521    s: &[S],
1522) -> f64
1523where
1524    <Mmut as Matrix>::R: RowMut,
1525    Mmut::R: RowEntity + TryConvertFrom1<[S], ModOp::M>,
1526    Mmut::MatElement: PrimInt + ToPrimitive + Debug,
1527{
1528    let ring_size = s.len();
1529    assert!(rlwe_ct.dimension() == (2, ring_size));
1530    assert!(encoded_m_ideal.as_ref().len() == ring_size);
1531
1532    // -(s * a)
1533    let q = mod_op.modulus();
1534    let mut s = Mmut::R::try_convert_from(s, &q);
1535    ntt_op.forward(s.as_mut());
1536    let mut a = Mmut::R::zeros(ring_size);
1537    a.as_mut().copy_from_slice(rlwe_ct.get_row_slice(0));
1538    ntt_op.forward(a.as_mut());
1539    mod_op.elwise_mul_mut(s.as_mut(), a.as_ref());
1540    mod_op.elwise_neg_mut(s.as_mut());
1541    ntt_op.backward(s.as_mut());
1542
1543    // m+e = b - s*a
1544    let mut m_plus_e = s;
1545    mod_op.elwise_add_mut(m_plus_e.as_mut(), rlwe_ct.get_row_slice(1));
1546
1547    // difference
1548    mod_op.elwise_sub_mut(m_plus_e.as_mut(), encoded_m_ideal.as_ref());
1549
1550    let mut max_diff_bits = f64::MIN;
1551    m_plus_e.as_ref().iter().for_each(|v| {
1552        let bits = (q.map_element_to_i64(v).to_f64().unwrap()).log2();
1553
1554        if max_diff_bits < bits {
1555            max_diff_bits = bits;
1556        }
1557    });
1558
1559    return max_diff_bits;
1560}
1561
1562#[cfg(test)]
1563pub(crate) mod tests {
1564    use std::{marker::PhantomData, ops::Mul, vec};
1565
1566    use itertools::{izip, Itertools};
1567    use rand::{thread_rng, Rng};
1568
1569    use crate::{
1570        backend::{GetModulus, ModInit, ModularOpsU64, Modulus, VectorOps},
1571        decomposer::{Decomposer, DefaultDecomposer, RlweDecomposer},
1572        ntt::{self, Ntt, NttBackendU64, NttInit},
1573        random::{DefaultSecureRng, NewWithSeed, RandomFillUniformInModulus},
1574        rgsw::{
1575            gen_rlwe_public_key, measure_noise, public_key_encrypt_rgsw, AutoKeyEvaluationDomain,
1576            RgswCiphertext, RgswCiphertextEvaluationDomain, RlweCiphertext, RlwePublicKey,
1577            SeededAutoKey, SeededRgswCiphertext, SeededRlweCiphertext, SeededRlwePublicKey,
1578        },
1579        utils::{generate_prime, negacyclic_mul, Stats, TryConvertFrom1},
1580        Matrix, Secret,
1581    };
1582
1583    use super::{
1584        decrypt_rlwe, galois_auto, galois_key_gen, generate_auto_map, public_key_encrypt_rlwe,
1585        rgsw_by_rgsw_inplace, rlwe_by_rgsw, secret_key_encrypt_rgsw, secret_key_encrypt_rlwe,
1586        RlweSecret,
1587    };
1588
1589    pub(crate) fn _sk_encrypt_rlwe<T: Modulus<Element = u64> + Clone>(
1590        m: &[u64],
1591        s: &[i32],
1592        ntt_op: &NttBackendU64,
1593        mod_op: &ModularOpsU64<T>,
1594    ) -> RlweCiphertext<Vec<Vec<u64>>, DefaultSecureRng> {
1595        let ring_size = m.len();
1596        let q = mod_op.modulus();
1597        assert!(s.len() == ring_size);
1598
1599        let mut rng = DefaultSecureRng::new();
1600        let mut rlwe_seed = [0u8; 32];
1601        rng.fill_bytes(&mut rlwe_seed);
1602        let mut seeded_rlwe_ct =
1603            SeededRlweCiphertext::<_, [u8; 32], _>::empty(ring_size as usize, rlwe_seed, q.clone());
1604        let mut p_rng = DefaultSecureRng::new_seeded(rlwe_seed);
1605        secret_key_encrypt_rlwe(
1606            &m,
1607            &mut seeded_rlwe_ct.data,
1608            s,
1609            mod_op,
1610            ntt_op,
1611            &mut p_rng,
1612            &mut rng,
1613        );
1614
1615        RlweCiphertext::<Vec<Vec<u64>>, DefaultSecureRng>::from(&seeded_rlwe_ct)
1616    }
1617
1618    // Encrypt m as RGSW ciphertext RGSW(m) using supplied public key
1619    pub(crate) fn _pk_encrypt_rgsw<T: Modulus<Element = u64> + Clone>(
1620        m: &[u64],
1621        public_key: &RlwePublicKey<Vec<Vec<u64>>, DefaultSecureRng>,
1622        decomposer: &(DefaultDecomposer<u64>, DefaultDecomposer<u64>),
1623        mod_op: &ModularOpsU64<T>,
1624        ntt_op: &NttBackendU64,
1625    ) -> RgswCiphertext<Vec<Vec<u64>>, T> {
1626        let (_, ring_size) = Matrix::dimension(&public_key.data);
1627        let gadget_vector_a = decomposer.a().gadget_vector();
1628        let gadget_vector_b = decomposer.b().gadget_vector();
1629
1630        let mut rng = DefaultSecureRng::new();
1631
1632        assert!(m.len() == ring_size);
1633
1634        // public key encrypt RGSW(m1)
1635        let mut rgsw_ct = RgswCiphertext::empty(ring_size, decomposer, mod_op.modulus().clone());
1636        public_key_encrypt_rgsw(
1637            &mut rgsw_ct.data,
1638            m,
1639            &public_key.data,
1640            &gadget_vector_a,
1641            &gadget_vector_b,
1642            mod_op,
1643            ntt_op,
1644            &mut rng,
1645        );
1646
1647        rgsw_ct
1648    }
1649
1650    /// Encrypts m as RGSW ciphertext RGSW(m) using supplied secret key. Returns
1651    /// unseeded RGSW ciphertext in coefficient domain
1652    pub(crate) fn _sk_encrypt_rgsw<T: Modulus<Element = u64> + Clone>(
1653        m: &[u64],
1654        s: &[i32],
1655        decomposer: &(DefaultDecomposer<u64>, DefaultDecomposer<u64>),
1656        mod_op: &ModularOpsU64<T>,
1657        ntt_op: &NttBackendU64,
1658    ) -> SeededRgswCiphertext<Vec<Vec<u64>>, [u8; 32], T> {
1659        let ring_size = s.len();
1660        assert!(m.len() == s.len());
1661
1662        let q = mod_op.modulus();
1663
1664        let gadget_vector_a = decomposer.a().gadget_vector();
1665        let gadget_vector_b = decomposer.b().gadget_vector();
1666
1667        let mut rng = DefaultSecureRng::new();
1668        let mut rgsw_seed = [0u8; 32];
1669        rng.fill_bytes(&mut rgsw_seed);
1670        let mut seeded_rgsw_ct = SeededRgswCiphertext::<Vec<Vec<u64>>, [u8; 32], T>::empty(
1671            ring_size as usize,
1672            decomposer,
1673            rgsw_seed,
1674            q.clone(),
1675        );
1676        let mut p_rng = DefaultSecureRng::new_seeded(rgsw_seed);
1677        secret_key_encrypt_rgsw(
1678            &mut seeded_rgsw_ct.data,
1679            m,
1680            &gadget_vector_a,
1681            &gadget_vector_b,
1682            s,
1683            mod_op,
1684            ntt_op,
1685            &mut p_rng,
1686            &mut rng,
1687        );
1688        seeded_rgsw_ct
1689    }
1690
1691    /// Prints noise in RGSW ciphertext RGSW(m).
1692    ///
1693    /// - rgsw_ct: RGSW ciphertext in coefficient domain
1694    pub(crate) fn _measure_noise_rgsw<T: Modulus<Element = u64> + Clone>(
1695        rgsw_ct: &[Vec<u64>],
1696        m: &[u64],
1697        s: &[i32],
1698        decomposer: &(DefaultDecomposer<u64>, DefaultDecomposer<u64>),
1699        q: &T,
1700    ) {
1701        let gadget_vector_a = decomposer.a().gadget_vector();
1702        let gadget_vector_b = decomposer.b().gadget_vector();
1703        let d_a = gadget_vector_a.len();
1704        let d_b = gadget_vector_b.len();
1705        let ring_size = s.len();
1706        assert!(Matrix::dimension(&rgsw_ct) == (d_a * 2 + d_b * 2, ring_size));
1707        assert!(m.len() == ring_size);
1708
1709        let mod_op = ModularOpsU64::new(q.clone());
1710        let ntt_op = NttBackendU64::new(q, ring_size);
1711
1712        let mul_mod =
1713            |a: &u64, b: &u64| ((*a as u128 * *b as u128) % q.q().unwrap() as u128) as u64;
1714        let s_poly = Vec::<u64>::try_convert_from(s, q);
1715        let mut neg_s = s_poly.clone();
1716        mod_op.elwise_neg_mut(neg_s.as_mut());
1717        let neg_sm0m1 = negacyclic_mul(&neg_s, &m, mul_mod, q.q().unwrap());
1718
1719        // RLWE(\beta^j -s * m)
1720        for j in 0..d_a {
1721            let ideal_m = {
1722                // RLWE(\beta^j -s * m)
1723                let mut beta_neg_sm0m1 = vec![0u64; ring_size as usize];
1724                mod_op.elwise_scalar_mul(beta_neg_sm0m1.as_mut(), &neg_sm0m1, &gadget_vector_a[j]);
1725                beta_neg_sm0m1
1726            };
1727
1728            let mut rlwe = vec![vec![0u64; ring_size as usize]; 2];
1729            rlwe[0].copy_from_slice(rgsw_ct.get_row_slice(j));
1730            rlwe[1].copy_from_slice(rgsw_ct.get_row_slice(d_a + j));
1731            let noise = measure_noise(&rlwe, &ideal_m, &ntt_op, &mod_op, s);
1732
1733            println!(r"Noise RLWE(\beta^{j} -sm0m1): {noise}");
1734        }
1735
1736        // RLWE(\beta^j  m)
1737        for j in 0..d_b {
1738            let ideal_m = {
1739                // RLWE(\beta^j  m)
1740                let mut beta_m0m1 = vec![0u64; ring_size as usize];
1741                mod_op.elwise_scalar_mul(beta_m0m1.as_mut(), &m, &gadget_vector_b[j]);
1742                beta_m0m1
1743            };
1744
1745            let mut rlwe = vec![vec![0u64; ring_size as usize]; 2];
1746            rlwe[0].copy_from_slice(rgsw_ct.get_row_slice(d_a * 2 + j));
1747            rlwe[1].copy_from_slice(rgsw_ct.get_row_slice(d_a * 2 + d_b + j));
1748            let noise = measure_noise(&rlwe, &ideal_m, &ntt_op, &mod_op, s);
1749
1750            println!(r"Noise RLWE(\beta^{j} m0m1): {noise}");
1751        }
1752    }
1753
1754    #[test]
1755    fn rlwe_encrypt_decryption() {
1756        let logq = 50;
1757        let logp = 2;
1758        let ring_size = 1 << 4;
1759        let q = generate_prime(logq, ring_size, 1u64 << logq).unwrap();
1760        let p = 1u64 << logp;
1761
1762        let mut rng = DefaultSecureRng::new();
1763
1764        let s = RlweSecret::random((ring_size >> 1) as usize, ring_size as usize);
1765
1766        // sample m0
1767        let mut m0 = vec![0u64; ring_size as usize];
1768        RandomFillUniformInModulus::<[u64], u64>::random_fill(
1769            &mut rng,
1770            &(1u64 << logp),
1771            m0.as_mut_slice(),
1772        );
1773
1774        let ntt_op = NttBackendU64::new(&q, ring_size as usize);
1775        let mod_op = ModularOpsU64::new(q);
1776
1777        // encrypt m0
1778        let encoded_m = m0
1779            .iter()
1780            .map(|v| (((*v as f64) * q as f64) / (p as f64)).round() as u64)
1781            .collect_vec();
1782        let rlwe_in_ct = _sk_encrypt_rlwe(&encoded_m, s.values(), &ntt_op, &mod_op);
1783
1784        let mut encoded_m_back = vec![0u64; ring_size as usize];
1785        decrypt_rlwe(
1786            &rlwe_in_ct,
1787            s.values(),
1788            &mut encoded_m_back,
1789            &ntt_op,
1790            &mod_op,
1791        );
1792        let m_back = encoded_m_back
1793            .iter()
1794            .map(|v| (((*v as f64 * p as f64) / q as f64).round() as u64) % p)
1795            .collect_vec();
1796        assert_eq!(m0, m_back);
1797
1798        // let noise = measure_noise(&rlwe_in_ct, &encoded_m, &ntt_op, &mod_op,
1799        // s.values()); println!("Noise: {noise}");
1800    }
1801
1802    #[test]
1803    fn rlwe_by_rgsw_works() {
1804        let logq = 50;
1805        let logp = 2;
1806        let ring_size = 1 << 9;
1807        let q = generate_prime(logq, ring_size, 1u64 << logq).unwrap();
1808        let p: u64 = 1u64 << logp;
1809
1810        let mut rng = DefaultSecureRng::new_seeded([0u8; 32]);
1811
1812        let s = RlweSecret::random((ring_size >> 1) as usize, ring_size as usize);
1813
1814        let mut m0 = vec![0u64; ring_size as usize];
1815        RandomFillUniformInModulus::<[u64], _>::random_fill(
1816            &mut rng,
1817            &(1u64 << logp),
1818            m0.as_mut_slice(),
1819        );
1820        let mut m1 = vec![0u64; ring_size as usize];
1821        m1[thread_rng().gen_range(0..ring_size) as usize] = 1;
1822
1823        let ntt_op = NttBackendU64::new(&q, ring_size as usize);
1824        let mod_op = ModularOpsU64::new(q);
1825        let d_rgsw = 10;
1826        let logb = 5;
1827        let decomposer = (
1828            DefaultDecomposer::new(q, logb, d_rgsw),
1829            DefaultDecomposer::new(q, logb, d_rgsw),
1830        );
1831
1832        // create public key
1833        let mut pk_seed = [0u8; 32];
1834        rng.fill_bytes(&mut pk_seed);
1835        let mut pk_prng = DefaultSecureRng::new_seeded(pk_seed);
1836        let mut seeded_pk =
1837            SeededRlwePublicKey::<Vec<u64>, _>::empty(ring_size as usize, pk_seed, q);
1838        gen_rlwe_public_key(
1839            &mut seeded_pk.data,
1840            s.values(),
1841            &ntt_op,
1842            &mod_op,
1843            &mut pk_prng,
1844            &mut rng,
1845        );
1846        let pk = RlwePublicKey::<Vec<Vec<u64>>, DefaultSecureRng>::from(&seeded_pk);
1847
1848        // Encrypt m1 as RGSW(m1)
1849        let rgsw_ct = {
1850            //TODO(Jay): Figure out better way to test secret key and public key variant of
1851            // RGSW ciphertext encryption within the same test
1852
1853            if true {
1854                // Encryption m1 as RGSW(m1) using secret key
1855                let seeded_rgsw_ct =
1856                    _sk_encrypt_rgsw(&m1, s.values(), &decomposer, &mod_op, &ntt_op);
1857                RgswCiphertextEvaluationDomain::<Vec<Vec<u64>>, DefaultSecureRng, NttBackendU64>::from(&seeded_rgsw_ct)
1858            } else {
1859                // Encrypt m1 as RGSW(m1) using public key
1860                let rgsw_ct = _pk_encrypt_rgsw(&m1, &pk, &decomposer, &mod_op, &ntt_op);
1861                RgswCiphertextEvaluationDomain::<_, DefaultSecureRng, NttBackendU64>::from(&rgsw_ct)
1862            }
1863        };
1864
1865        // Encrypt m0 as RLWE(m0)
1866        let mut rlwe_in_ct = {
1867            let encoded_m = m0
1868                .iter()
1869                .map(|v| (((*v as f64) * q as f64) / (p as f64)).round() as u64)
1870                .collect_vec();
1871
1872            _sk_encrypt_rlwe(&encoded_m, s.values(), &ntt_op, &mod_op)
1873        };
1874
1875        // RLWE(m0m1) = RLWE(m0) x RGSW(m1)
1876        let mut scratch_space = vec![
1877            vec![0u64; ring_size as usize];
1878            std::cmp::max(
1879                decomposer.a().decomposition_count(),
1880                decomposer.b().decomposition_count()
1881            ) + 2
1882        ];
1883        rlwe_by_rgsw(
1884            &mut rlwe_in_ct,
1885            &rgsw_ct.data,
1886            &mut scratch_space,
1887            &decomposer,
1888            &ntt_op,
1889            &mod_op,
1890        );
1891
1892        // Decrypt RLWE(m0m1)
1893        let mut encoded_m0m1_back = vec![0u64; ring_size as usize];
1894        decrypt_rlwe(
1895            &rlwe_in_ct,
1896            s.values(),
1897            &mut encoded_m0m1_back,
1898            &ntt_op,
1899            &mod_op,
1900        );
1901        let m0m1_back = encoded_m0m1_back
1902            .iter()
1903            .map(|v| (((*v as f64 * p as f64) / (q as f64)).round() as u64) % p)
1904            .collect_vec();
1905
1906        let mul_mod = |v0: &u64, v1: &u64| (v0 * v1) % p;
1907        let m0m1 = negacyclic_mul(&m0, &m1, mul_mod, p);
1908
1909        // {
1910        //     // measure noise
1911        //     let encoded_m_ideal = m0m1
1912        //         .iter()
1913        //         .map(|v| (((*v as f64) * q as f64) / (p as f64)).round() as u64)
1914        //         .collect_vec();
1915
1916        //     let noise = measure_noise(&rlwe_in_ct, &encoded_m_ideal, &ntt_op,
1917        // &mod_op, s.values());     println!("Noise RLWE(m0m1)(=
1918        // RLWE(m0)xRGSW(m1)) : {noise}"); }
1919
1920        assert!(
1921            m0m1 == m0m1_back,
1922            "Expected {:?} \n Got {:?}",
1923            m0m1,
1924            m0m1_back
1925        );
1926    }
1927
1928    #[test]
1929    fn galois_auto_works() {
1930        let logq = 55;
1931        let ring_size = 1 << 11;
1932        let q = generate_prime(logq, 2 * ring_size, 1u64 << logq).unwrap();
1933        let logp = 3;
1934        let p = 1u64 << logp;
1935        let d_rgsw = 5;
1936        let logb = 11;
1937
1938        let mut rng = DefaultSecureRng::new();
1939        let s = RlweSecret::random((ring_size >> 1) as usize, ring_size as usize);
1940
1941        let mut m = vec![0u64; ring_size as usize];
1942        RandomFillUniformInModulus::random_fill(&mut rng, &p, m.as_mut_slice());
1943        let encoded_m = m
1944            .iter()
1945            .map(|v| (((*v as f64 * q as f64) / (p as f64)).round() as u64))
1946            .collect_vec();
1947
1948        let ntt_op = NttBackendU64::new(&q, ring_size as usize);
1949        let mod_op = ModularOpsU64::new(q);
1950
1951        // RLWE_{s}(m)
1952        let mut seed_rlwe = [0u8; 32];
1953        rng.fill_bytes(&mut seed_rlwe);
1954        let mut seeded_rlwe_m = SeededRlweCiphertext::empty(ring_size as usize, seed_rlwe, q);
1955        let mut p_rng = DefaultSecureRng::new_seeded(seed_rlwe);
1956        secret_key_encrypt_rlwe(
1957            &encoded_m,
1958            &mut seeded_rlwe_m.data,
1959            s.values(),
1960            &mod_op,
1961            &ntt_op,
1962            &mut p_rng,
1963            &mut rng,
1964        );
1965        let mut rlwe_m = RlweCiphertext::<Vec<Vec<u64>>, DefaultSecureRng>::from(&seeded_rlwe_m);
1966
1967        let auto_k = -125;
1968
1969        // Generate galois key to key switch from s^k to s
1970        let decomposer = DefaultDecomposer::new(q, logb, d_rgsw);
1971        let mut seed_auto = [0u8; 32];
1972        rng.fill_bytes(&mut seed_auto);
1973        let mut seeded_auto_key =
1974            SeededAutoKey::empty(ring_size as usize, &decomposer, seed_auto, q);
1975        let mut p_rng = DefaultSecureRng::new_seeded(seed_auto);
1976        let gadget_vector = decomposer.gadget_vector();
1977        galois_key_gen(
1978            &mut seeded_auto_key.data,
1979            s.values(),
1980            auto_k,
1981            &gadget_vector,
1982            &mod_op,
1983            &ntt_op,
1984            &mut p_rng,
1985            &mut rng,
1986        );
1987        let auto_key =
1988            AutoKeyEvaluationDomain::<Vec<Vec<u64>>, DefaultSecureRng, NttBackendU64>::from(
1989                &seeded_auto_key,
1990            );
1991
1992        // Send RLWE_{s}(m) -> RLWE_{s}(m^k)
1993        let mut scratch_space = vec![vec![0u64; ring_size as usize]; d_rgsw + 2];
1994        let (auto_map_index, auto_map_sign) = generate_auto_map(ring_size as usize, auto_k);
1995        galois_auto(
1996            &mut rlwe_m,
1997            &auto_key.data,
1998            &mut scratch_space,
1999            &auto_map_index,
2000            &auto_map_sign,
2001            &mod_op,
2002            &ntt_op,
2003            &decomposer,
2004        );
2005
2006        let rlwe_m_k = rlwe_m;
2007
2008        // Decrypt RLWE_{s}(m^k) and check
2009        let mut encoded_m_k_back = vec![0u64; ring_size as usize];
2010        decrypt_rlwe(
2011            &rlwe_m_k,
2012            s.values(),
2013            &mut encoded_m_k_back,
2014            &ntt_op,
2015            &mod_op,
2016        );
2017        let m_k_back = encoded_m_k_back
2018            .iter()
2019            .map(|v| (((*v as f64 * p as f64) / q as f64).round() as u64) % p)
2020            .collect_vec();
2021
2022        let mut m_k = vec![0u64; ring_size as usize];
2023        // Send \delta m -> \delta m^k
2024        izip!(m.iter(), auto_map_index.iter(), auto_map_sign.iter()).for_each(
2025            |(v, to_index, sign)| {
2026                if !*sign {
2027                    m_k[*to_index] = (p - *v) % p;
2028                } else {
2029                    m_k[*to_index] = *v;
2030                }
2031            },
2032        );
2033
2034        {
2035            let encoded_m_k = m_k
2036                .iter()
2037                .map(|v| ((*v as f64 * q as f64) / p as f64).round() as u64)
2038                .collect_vec();
2039
2040            let noise = measure_noise(&rlwe_m_k, &encoded_m_k, &ntt_op, &mod_op, s.values());
2041            println!("Ksk noise: {noise}");
2042        }
2043
2044        assert_eq!(m_k_back, m_k);
2045    }
2046
2047    #[test]
2048    fn sk_rgsw_by_rgsw() {
2049        let logq = 60;
2050        let logp = 2;
2051        let ring_size = 1 << 11;
2052        let q = generate_prime(logq, ring_size, 1u64 << logq).unwrap();
2053        let p = 1u64 << logp;
2054        let d_rgsw = 5;
2055        let logb = 12;
2056
2057        let s = RlweSecret::random((ring_size >> 1) as usize, ring_size as usize);
2058
2059        let mut rng = DefaultSecureRng::new();
2060        let ntt_op = NttBackendU64::new(&q, ring_size as usize);
2061        let mod_op = ModularOpsU64::new(q);
2062        let decomposer = (
2063            DefaultDecomposer::new(q, logb, d_rgsw),
2064            DefaultDecomposer::new(q, logb, d_rgsw),
2065        );
2066
2067        let mul_mod = |a: &u64, b: &u64| ((*a as u128 * *b as u128) % q as u128) as u64;
2068
2069        let mut carry_m = vec![0u64; ring_size as usize];
2070        carry_m[thread_rng().gen_range(0..ring_size) as usize] = 1;
2071
2072        // RGSW(carry_m)
2073        let mut rgsw_carrym = {
2074            let seeded_rgsw = _sk_encrypt_rgsw(&carry_m, s.values(), &decomposer, &mod_op, &ntt_op);
2075            let mut rgsw_eval =
2076                RgswCiphertextEvaluationDomain::<_, DefaultSecureRng, NttBackendU64>::from(
2077                    &seeded_rgsw,
2078                );
2079            rgsw_eval
2080                .data
2081                .iter_mut()
2082                .for_each(|ri| ntt_op.backward(ri.as_mut()));
2083            rgsw_eval.data
2084        };
2085
2086        let mut scratch_matrix = vec![
2087            vec![0u64; ring_size as usize];
2088            decomposer.a().decomposition_count() * 2
2089                + decomposer.b().decomposition_count() * 2
2090                + std::cmp::max(
2091                    decomposer.a().decomposition_count(),
2092                    decomposer.b().decomposition_count()
2093                )
2094        ];
2095
2096        _measure_noise_rgsw(&rgsw_carrym, &carry_m, s.values(), &decomposer, &q);
2097
2098        for i in 0..2 {
2099            let mut m = vec![0u64; ring_size as usize];
2100            m[thread_rng().gen_range(0..ring_size) as usize] = if (i & 1) == 1 { q - 1 } else { 1 };
2101            let rgsw_m = RgswCiphertextEvaluationDomain::<_, DefaultSecureRng, NttBackendU64>::from(
2102                &_sk_encrypt_rgsw(&m, s.values(), &decomposer, &mod_op, &ntt_op),
2103            );
2104            rgsw_by_rgsw_inplace(
2105                &mut rgsw_carrym,
2106                &rgsw_m.data,
2107                &decomposer,
2108                &mut scratch_matrix,
2109                &ntt_op,
2110                &mod_op,
2111            );
2112
2113            // measure noise
2114            carry_m = negacyclic_mul(&carry_m, &m, mul_mod, q);
2115            println!("########### Noise RGSW(carrym) in {i}^th loop ###########");
2116            _measure_noise_rgsw(&rgsw_carrym, &carry_m, s.values(), &decomposer, &q);
2117        }
2118        {
2119            // RLWE(m) x RGSW(carry_m)
2120            let mut m = vec![0u64; ring_size as usize];
2121            RandomFillUniformInModulus::random_fill(&mut rng, &q, m.as_mut_slice());
2122            let mut rlwe_ct = _sk_encrypt_rlwe(&m, s.values(), &ntt_op, &mod_op);
2123
2124            // send rgsw to evaluation domain
2125            rgsw_carrym
2126                .iter_mut()
2127                .for_each(|ri| ntt_op.forward(ri.as_mut_slice()));
2128
2129            rlwe_by_rgsw(
2130                &mut rlwe_ct,
2131                &rgsw_carrym,
2132                &mut scratch_matrix,
2133                &decomposer,
2134                &ntt_op,
2135                &mod_op,
2136            );
2137            let m_expected = negacyclic_mul(&carry_m, &m, mul_mod, q);
2138            let noise = measure_noise(&rlwe_ct, &m_expected, &ntt_op, &mod_op, s.values());
2139            println!("RLWE(m) x RGSW(carry_m): {noise}");
2140        }
2141    }
2142
2143    #[test]
2144    fn some_work() {
2145        let logq = 50;
2146        let ring_size = 1 << 11;
2147        let q = generate_prime(logq, ring_size, 1u64 << logq).unwrap();
2148        let d_rgsw = 10;
2149        let logb = 5;
2150        let decomposer = (
2151            DefaultDecomposer::new(q, logb, d_rgsw),
2152            DefaultDecomposer::new(q, logb, d_rgsw),
2153        );
2154
2155        let ntt_op = NttBackendU64::new(&q, ring_size as usize);
2156        let mod_op = ModularOpsU64::new(q);
2157        let mut rng = DefaultSecureRng::new_seeded([0u8; 32]);
2158
2159        let s = RlweSecret::random((ring_size >> 1) as usize, ring_size as usize);
2160
2161        let mut check = Stats { samples: vec![] };
2162
2163        for _ in 0..100 {
2164            let mut m0 = vec![0u64; ring_size as usize];
2165            m0[thread_rng().gen_range(0..ring_size) as usize] = 1;
2166            let mut m1 = vec![0u64; ring_size as usize];
2167            m1[thread_rng().gen_range(0..ring_size) as usize] = 1;
2168
2169            let mut rgsw_ct0 = {
2170                let seeded_rgsw_ct =
2171                    _sk_encrypt_rgsw(&m0, s.values(), &decomposer, &mod_op, &ntt_op);
2172                RgswCiphertextEvaluationDomain::<Vec<Vec<u64>>, DefaultSecureRng, NttBackendU64>::from(&seeded_rgsw_ct)
2173            };
2174            let rgsw_ct1 = {
2175                let seeded_rgsw_ct =
2176                    _sk_encrypt_rgsw(&m1, s.values(), &decomposer, &mod_op, &ntt_op);
2177                RgswCiphertextEvaluationDomain::<Vec<Vec<u64>>, DefaultSecureRng, NttBackendU64>::from(&seeded_rgsw_ct)
2178            };
2179
2180            // RGSW x RGSW
2181            // send RGSW(m0) to coefficient domain
2182            rgsw_ct0
2183                .data
2184                .iter_mut()
2185                .for_each(|r| ntt_op.backward(r.as_mut_slice()));
2186            let mut scratch_matrix = vec![
2187                vec![0u64; ring_size as usize];
2188                std::cmp::max(
2189                    decomposer.a().decomposition_count(),
2190                    decomposer.b().decomposition_count()
2191                ) + decomposer.a().decomposition_count() * 2
2192                    + decomposer.b().decomposition_count() * 2
2193            ];
2194            rgsw_by_rgsw_inplace(
2195                &mut rgsw_ct0.data,
2196                &rgsw_ct1.data,
2197                &decomposer,
2198                &mut scratch_matrix,
2199                &ntt_op,
2200                &mod_op,
2201            );
2202            let mut rgsw_m0m1 = rgsw_ct0;
2203            // Back to Evaluation for RLWExRGSW
2204            rgsw_m0m1
2205                .data
2206                .iter_mut()
2207                .for_each(|r| ntt_op.forward(r.as_mut_slice()));
2208
2209            // Sample m2, encrypt it as RLWE(m2) and multiply RLWE(m2)xRGSW(m0m1)
2210            let mut m2 = vec![0u64; ring_size as usize];
2211            RandomFillUniformInModulus::random_fill(&mut rng, &q, m2.as_mut_slice());
2212            let mut rlwe_in_ct = { _sk_encrypt_rlwe(&m2, s.values(), &ntt_op, &mod_op) };
2213            let mut scratch_space = vec![
2214                vec![0u64; ring_size as usize];
2215                std::cmp::max(
2216                    decomposer.a().decomposition_count(),
2217                    decomposer.b().decomposition_count()
2218                ) + 2
2219            ];
2220            rlwe_by_rgsw(
2221                &mut rlwe_in_ct,
2222                &rgsw_m0m1.data,
2223                &mut scratch_space,
2224                &decomposer,
2225                &ntt_op,
2226                &mod_op,
2227            );
2228
2229            // Decrypt RLWE(m0m1m2)
2230            let mut m0m1m2_back = vec![0u64; ring_size as usize];
2231            decrypt_rlwe(&rlwe_in_ct, s.values(), &mut m0m1m2_back, &ntt_op, &mod_op);
2232
2233            // Calculate m0m1m2
2234            let mul_mod = |v0: &u64, v1: &u64| ((*v0 as u128 * *v1 as u128) % q as u128) as u64;
2235            let m0m1 = negacyclic_mul(&m0, &m1, mul_mod, q);
2236            let m0m1m2 = negacyclic_mul(&m2, &m0m1, mul_mod, q);
2237
2238            // diff
2239            mod_op.elwise_sub_mut(m0m1m2_back.as_mut_slice(), m0m1m2.as_ref());
2240
2241            check.add_more(&Vec::<i64>::try_convert_from(&m0m1m2_back, &q));
2242        }
2243
2244        println!("Std: {}", check.std_dev().abs().log2());
2245    }
2246}