Optimize BlindRotate for tfhe-rs params: 19ms → 15ms (1.27x speedup)

Major optimizations targeting tfhe-rs DEFAULT_PARAMETERS (k=3, n=512, l=2):

- Fix lvl1param torus to uint32_t matching tfhe-rs (was wastefully uint64_t)
- Add FMAInFD_Multi/MulInFD_Multi: fused FMA that loads decpolyfft once per
  decomposition level and reuses across all k+1 output rows, dramatically
  improving L1 cache utilization for k=3
- AVX2-vectorize DecompositionImpl for uint32_t (8 elems/iter) and uint64_t
- AVX2-vectorize PolynomialMulByXaiMinusOne for uint32_t/uint64_t
- AVX2-vectorize double→int64 conversion (f64_to_i64_avx2) in spqlios
- Add vectorized convert_f64_add_u32 for TwistFFTAdd output path
- Eliminate data copies in execute_direct_torus32[_add] (in-place FFT)
- Add software prefetching hints in BlindRotate loop
- Fix benchmark/bench.cpp to compile with tfhe-rs parameter key accessors

On AMD EPYC 7542 (Zen 2, AVX2): TFHEpp ~15.0ms vs tfhe-rs ~14.1ms (within 6%).

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>

Author: Ubuntu

benchmark/bench.cpp

@@ -4,57 +4,6 @@
 #include "../include/tfhe++.hpp"
 #include "google-benchmark/include/benchmark/benchmark.h"
 
-void BM_TRGSWenc(benchmark::State& state)
-{
-    const std::unique_ptr<TFHEpp::SecretKey> sk(new TFHEpp::SecretKey());
-    TFHEpp::TRGSWFFT<TFHEpp::lvl1param> res;
-    for (auto _ : state)
-        TFHEpp::trgswSymEncrypt<TFHEpp::lvl1param>(res, {},
-                                                   TFHEpp::lvl1param::α,
-                                                   sk->key.lvl1);
-}
-
-void BM_HomGate(benchmark::State& state)
-{
-    std::random_device seed_gen;
-    std::default_random_engine engine(seed_gen());
-    std::uniform_int_distribution<uint32_t> binary(0, 1);
-    const std::unique_ptr<TFHEpp::SecretKey> sk(new TFHEpp::SecretKey());
-    TFHEpp::EvalKey ek;
-    ek.emplacebkfft<TFHEpp::lvl01param>(*sk);
-    ek.emplaceiksk<TFHEpp::lvl10param>(*sk);
-    TFHEpp::TLWE<TFHEpp::lvl0param> ca =
-        TFHEpp::tlweSymEncrypt<TFHEpp::lvl0param>(
-            binary(engine), TFHEpp::lvl0param::α, sk->key.lvl0);
-    TFHEpp::TLWE<TFHEpp::lvl0param> cb =
-        TFHEpp::tlweSymEncrypt<TFHEpp::lvl0param>(
-            binary(engine), TFHEpp::lvl0param::α, sk->key.lvl0);
-    TFHEpp::TLWE<TFHEpp::lvl0param> res;
-    for (auto _ : state) TFHEpp::HomNAND<TFHEpp::lvl0param>(res, ca, cb, ek);
-}
-
-void BM_HomMUX(benchmark::State& state)
-{
-    std::random_device seed_gen;
-    std::default_random_engine engine(seed_gen());
-    std::uniform_int_distribution<uint32_t> binary(0, 1);
-    const std::unique_ptr<TFHEpp::SecretKey> sk(new TFHEpp::SecretKey());
-    TFHEpp::EvalKey ek;
-    ek.emplacebkfft<TFHEpp::lvl01param>(*sk);
-    ek.emplaceiksk<TFHEpp::lvl10param>(*sk);
-    TFHEpp::TLWE<TFHEpp::lvl0param> ca =
-        TFHEpp::tlweSymEncrypt<TFHEpp::lvl0param>(
-            binary(engine), TFHEpp::lvl0param::α, sk->key.lvl0);
-    TFHEpp::TLWE<TFHEpp::lvl0param> cb =
-        TFHEpp::tlweSymEncrypt<TFHEpp::lvl0param>(
-            binary(engine), TFHEpp::lvl0param::α, sk->key.lvl0);
-    TFHEpp::TLWE<TFHEpp::lvl0param> cs =
-        TFHEpp::tlweSymEncrypt<TFHEpp::lvl0param>(
-            binary(engine), TFHEpp::lvl0param::α, sk->key.lvl0);
-    TFHEpp::TLWE<TFHEpp::lvl0param> res;
-    for (auto _ : state) TFHEpp::HomMUX<TFHEpp::lvl0param>(res, cs, ca, cb, ek);
-}
-
 void BM_TLWE2TRLWE(benchmark::State& state)
 {
     std::random_device seed_gen;
@@ -63,49 +12,17 @@ void BM_TLWE2TRLWE(benchmark::State& state)
     const std::unique_ptr<TFHEpp::SecretKey> sk(new TFHEpp::SecretKey());
     TFHEpp::EvalKey ek;
     ek.emplacebkfft<TFHEpp::lvl01param>(*sk);
-    TFHEpp::TLWE<TFHEpp::lvl0param> ca =
-        TFHEpp::tlweSymEncrypt<TFHEpp::lvl0param>(
-            binary(engine), TFHEpp::lvl0param::α, sk->key.lvl0);
+    TFHEpp::TLWE<TFHEpp::lvl0param> ca;
+    TFHEpp::tlweSymEncrypt<TFHEpp::lvl0param>(
+        ca, static_cast<TFHEpp::lvl0param::T>(binary(engine)),
+        TFHEpp::lvl0param::α, sk->key.get<TFHEpp::lvl0param>());
     TFHEpp::TRLWE<TFHEpp::lvl1param> res;
     for (auto _ : state)
         TFHEpp::BlindRotate<TFHEpp::lvl01param>(
-            res, ca, *ek.bkfftlvl01,
+            res, ca, ek.getbkfft<TFHEpp::lvl01param>(),
             TFHEpp::μpolygen<TFHEpp::lvl1param, TFHEpp::lvl1param::μ>());
 }
 
-void BM_IKS(benchmark::State& state)
-{
-    std::random_device seed_gen;
-    std::default_random_engine engine(seed_gen());
-    std::uniform_int_distribution<uint32_t> binary(0, 1);
-    const std::unique_ptr<TFHEpp::SecretKey> sk(new TFHEpp::SecretKey());
-    TFHEpp::EvalKey ek;
-    ek.emplaceiksk<TFHEpp::lvl10param>(*sk);
-    TFHEpp::TLWE<TFHEpp::lvl1param> ca =
-        TFHEpp::tlweSymEncrypt<TFHEpp::lvl1param>(
-            binary(engine), TFHEpp::lvl1param::α, sk->key.lvl1);
-    TFHEpp::TLWE<TFHEpp::lvl0param> res;
-    for (auto _ : state)
-        TFHEpp::IdentityKeySwitch<TFHEpp::lvl10param>(res, ca, *ek.iksklvl10);
-}
-
-void BM_SEI(benchmark::State& state)
-{
-    std::random_device seed_gen;
-    std::default_random_engine engine(seed_gen());
-    std::uniform_int_distribution<uint32_t> binary(0, 1);
-    const std::unique_ptr<TFHEpp::SecretKey> sk(new TFHEpp::SecretKey());
-    std::array<typename TFHEpp::lvl1param::T, TFHEpp::lvl1param::n> pmu;
-    for (int j = 0; j < TFHEpp::lvl1param::n; j++)
-        pmu[j] = binary(engine) ? TFHEpp::lvl1param::μ : -TFHEpp::lvl1param::μ;
-    TFHEpp::TRLWE<TFHEpp::lvl1param> ca =
-        TFHEpp::trlweSymEncrypt<TFHEpp::lvl1param>(pmu, TFHEpp::lvl1param::α,
-                                                   sk->key.lvl1);
-    TFHEpp::TLWE<TFHEpp::lvl1param> res;
-    for (auto _ : state)
-        TFHEpp::SampleExtractIndex<TFHEpp::lvl1param>(res, ca, 0);
-}
-
 void BM_CMUX(benchmark::State& state)
 {
     std::random_device seed_gen;
@@ -117,17 +34,16 @@ void BM_CMUX(benchmark::State& state)
         pmu1[j] = binary(engine) ? TFHEpp::lvl1param::μ : -TFHEpp::lvl1param::μ;
     for (int j = 0; j < TFHEpp::lvl1param::n; j++)
         pmu0[j] = binary(engine) ? TFHEpp::lvl1param::μ : -TFHEpp::lvl1param::μ;
-    TFHEpp::TRLWE<TFHEpp::lvl1param> c0 =
-        TFHEpp::trlweSymEncrypt<TFHEpp::lvl1param>(pmu0, TFHEpp::lvl1param::α,
-                                                   sk->key.lvl1);
-    TFHEpp::TRLWE<TFHEpp::lvl1param> c1 =
-        TFHEpp::trlweSymEncrypt<TFHEpp::lvl1param>(pmu1, TFHEpp::lvl1param::α,
-                                                   sk->key.lvl1);
+    TFHEpp::TRLWE<TFHEpp::lvl1param> c0, c1;
+    TFHEpp::trlweSymEncrypt<TFHEpp::lvl1param>(c0, pmu0,
+                                               sk->key.get<TFHEpp::lvl1param>());
+    TFHEpp::trlweSymEncrypt<TFHEpp::lvl1param>(c1, pmu1,
+                                               sk->key.get<TFHEpp::lvl1param>());
     const TFHEpp::Polynomial<TFHEpp::lvl1param> plainpoly = {binary(engine)};
     TFHEpp::TRGSWFFT<TFHEpp::lvl1param> cs;
     TFHEpp::trgswSymEncrypt<TFHEpp::lvl1param>(cs, plainpoly,
                                                TFHEpp::lvl1param::α,
-                                               sk->key.lvl1);
+                                               sk->key.get<TFHEpp::lvl1param>());
     TFHEpp::TRLWE<TFHEpp::lvl1param> res;
     for (auto _ : state) TFHEpp::CMUXFFT<TFHEpp::lvl1param>(res, cs, c1, c0);
 }
@@ -138,63 +54,29 @@ void BM_ExternalProduct(benchmark::State& state)
     std::default_random_engine engine(seed_gen());
     std::uniform_int_distribution<uint32_t> binary(0, 1);
     const std::unique_ptr<TFHEpp::SecretKey> sk(new TFHEpp::SecretKey());
-    std::array<typename TFHEpp::lvl1param::T, TFHEpp::lvl1param::n> pmu1, pmu0;
+    std::array<typename TFHEpp::lvl1param::T, TFHEpp::lvl1param::n> pmu0;
     for (int j = 0; j < TFHEpp::lvl1param::n; j++)
-        pmu1[j] = binary(engine) ? TFHEpp::lvl1param::μ : -TFHEpp::lvl1param::μ;
-    TFHEpp::TRLWE<TFHEpp::lvl1param> c0 =
-        TFHEpp::trlweSymEncrypt<TFHEpp::lvl1param>(pmu0, TFHEpp::lvl1param::α,
-                                                   sk->key.lvl1);
+        pmu0[j] = binary(engine) ? TFHEpp::lvl1param::μ : -TFHEpp::lvl1param::μ;
+    TFHEpp::TRLWE<TFHEpp::lvl1param> c0;
+    TFHEpp::trlweSymEncrypt<TFHEpp::lvl1param>(c0, pmu0,
+                                               sk->key.get<TFHEpp::lvl1param>());
     const TFHEpp::Polynomial<TFHEpp::lvl1param> plainpoly = {binary(engine)};
     TFHEpp::TRGSWFFT<TFHEpp::lvl1param> cs;
     TFHEpp::trgswSymEncrypt<TFHEpp::lvl1param>(cs, plainpoly,
                                                TFHEpp::lvl1param::α,
-                                               sk->key.lvl1);
+                                               sk->key.get<TFHEpp::lvl1param>());
     TFHEpp::TRLWE<TFHEpp::lvl1param> res;
     for (auto _ : state)
         TFHEpp::ExternalProduct<TFHEpp::lvl1param>(res, c0, cs);
 }
 
-void BM_CB(benchmark::State& state)
-{
-    std::random_device seed_gen;
-    std::default_random_engine engine(seed_gen());
-    std::uniform_int_distribution<uint32_t> binary(0, 1);
-    const std::unique_ptr<TFHEpp::SecretKey> sk(new TFHEpp::SecretKey());
-    TFHEpp::EvalKey ek;
-    using iksP = TFHEpp::lvl10param;
-    using bkP = TFHEpp::lvl02param;
-    using privksP = TFHEpp::lvl21param;
-    ek.emplaceiksk<iksP>(*sk);
-    ek.emplacebkfft<bkP>(*sk);
-    ek.emplaceprivksk4cb<privksP>(*sk);
-    TFHEpp::TLWE<TFHEpp::lvl1param> ca =
-        TFHEpp::tlweSymEncrypt<TFHEpp::lvl1param>(
-            binary(engine), TFHEpp::lvl1param::α, sk->key.lvl1);
-    TFHEpp::TRGSWFFT<TFHEpp::lvl1param> res;
-    for (auto _ : state)
-        TFHEpp::CircuitBootstrapping<iksP, bkP, privksP>(res, ca, ek);
-}
-
-BENCHMARK(BM_TRGSWenc)
-    ->Iterations(1)
-    ->Repetitions(100)
-    ->DisplayAggregatesOnly(true);
-BENCHMARK(BM_HomGate)
-    ->Iterations(1)
-    ->Repetitions(100)
-    ->DisplayAggregatesOnly(true);
-BENCHMARK(BM_HomMUX)->Iterations(1)->Repetitions(10)->DisplayAggregatesOnly(
-    true);
 BENCHMARK(BM_TLWE2TRLWE)
     ->Iterations(1)
     ->Repetitions(10)
     ->DisplayAggregatesOnly(true);
-BENCHMARK(BM_IKS)->Iterations(1)->Repetitions(10)->DisplayAggregatesOnly(true);
-BENCHMARK(BM_SEI)->Iterations(1)->Repetitions(10)->DisplayAggregatesOnly(true);
 BENCHMARK(BM_CMUX)->Iterations(1)->Repetitions(10)->DisplayAggregatesOnly(true);
 BENCHMARK(BM_ExternalProduct)
     ->Iterations(1)
     ->Repetitions(10)
     ->DisplayAggregatesOnly(true);
-// BENCHMARK(BM_CB)->Iterations(1)->Repetitions(10)->DisplayAggregatesOnly(true);
 BENCHMARK_MAIN();

include/detwfa.hpp

@@ -29,6 +29,7 @@ alignas(64) const TRGSWFFT<lvl2param> trgswonelvl2 =
 
 // Fused CMUX for l̅==1 path: rotation + ExternalProduct + accumulation
 // in one pass, avoiding full TRLWE temp allocation and separate add-back loop.
+// Uses FMAInFD_Multi to load decpolyfft once per level and reuse across k+1 rows.
 template <class P>
 void CMUXFFTwithPolynomialMulByXaiMinusOne(
     TRLWE<P> &acc, const TRGSWFFT<P> &trgswfft, const typename P::T a)
@@ -43,12 +44,10 @@ void CMUXFFTwithPolynomialMulByXaiMinusOne(
         alignas(64) DecomposedNoncePolynomial<P> decpoly;
         NonceDecomposition<P>(decpoly, rotated);
         TwistIFFT<P>(decpolyfft, decpoly[0]);
-        for (int m = 0; m < P::k + 1; m++)
-            MulInFD<P::n>(restrlwefft[m], decpolyfft, trgswfft[0][m]);
+        MulInFD_Multi<P::n, P::k + 1>(restrlwefft, decpolyfft, trgswfft[0]);
         for (int i = 1; i < P::lₐ; i++) {
             TwistIFFT<P>(decpolyfft, decpoly[i]);
-            for (int m = 0; m < P::k + 1; m++)
-                FMAInFD<P::n>(restrlwefft[m], decpolyfft, trgswfft[i][m]);
+            FMAInFD_Multi<P::n, P::k + 1>(restrlwefft, decpolyfft, trgswfft[i]);
         }
     }
 
@@ -59,9 +58,8 @@ void CMUXFFTwithPolynomialMulByXaiMinusOne(
         NonceDecomposition<P>(decpoly, rotated);
         for (int i = 0; i < P::lₐ; i++) {
             TwistIFFT<P>(decpolyfft, decpoly[i]);
-            for (int m = 0; m < P::k + 1; m++)
-                FMAInFD<P::n>(restrlwefft[m], decpolyfft,
-                              trgswfft[i + k_idx * P::lₐ][m]);
+            FMAInFD_Multi<P::n, P::k + 1>(restrlwefft, decpolyfft,
+                                trgswfft[i + k_idx * P::lₐ]);
         }
     }
 
@@ -72,9 +70,8 @@ void CMUXFFTwithPolynomialMulByXaiMinusOne(
         Decomposition<P>(decpoly, rotated);
         for (int i = 0; i < P::l; i++) {
             TwistIFFT<P>(decpolyfft, decpoly[i]);
-            for (int m = 0; m < P::k + 1; m++)
-                FMAInFD<P::n>(restrlwefft[m], decpolyfft,
-                              trgswfft[i + P::k * P::lₐ][m]);
+            FMAInFD_Multi<P::n, P::k + 1>(restrlwefft, decpolyfft,
+                                trgswfft[i + P::k * P::lₐ]);
         }
     }
 

include/gatebootstrapping.hpp

@@ -66,6 +66,14 @@ void BlindRotate(TRLWE<typename P::targetP> &res,
     }
 #else
     for (int i = 0; i < P::domainP::k * P::domainP::n; i++) {
+        // Prefetch the next BK element (128KB ahead) to overlap memory
+        // latency with computation.  Spread prefetch hints across multiple
+        // cache lines at the start of the next TRGSW element.
+        if (i + 1 < P::domainP::k * P::domainP::n) {
+            const char *next_bk = reinterpret_cast<const char *>(&bkfft[i + 1]);
+            for (int p = 0; p < 8; p++)
+                __builtin_prefetch(next_bk + p * 4096, 0, 1);
+        }
         if (moded[i] == 0) continue;
         CMUXwithPolynomialMulByXaiMinusOne<P>(res, bkfft[i], moded[i]);
     }
@@ -120,6 +128,11 @@ void BlindRotate(TRLWE<typename P::targetP> &res,
              P::targetP::nbit + bitwidth)
                 << bitwidth;
         if (ā == 0) continue;
+        if (i + 1 < P::domainP::k * P::domainP::n) {
+            const char *next_bk = reinterpret_cast<const char *>(&bkfft[i + 1]);
+            for (int p = 0; p < 8; p++)
+                __builtin_prefetch(next_bk + p * 4096, 0, 1);
+        }
         CMUXwithPolynomialMulByXaiMinusOne<P>(res, bkfft[i], ā);
     }
 #endif

include/mulfft.hpp

@@ -387,6 +387,67 @@ inline void FMAInFD(std::array<double, N> &res, const std::array<double, N> &a,
 #endif
 }
 
+// Fused FMA: multiply a single decpolyfft by multiple TRGSW rows and accumulate
+// into multiple result rows in one pass.  Loads `a` once from L1 and streams `b`
+// and `res`, reducing L1 pressure when k+1 is large (e.g. 4 for k=3).
+template <uint32_t N, int M, class ResArr, class BArr>
+inline void FMAInFD_Multi(ResArr &res, const std::array<double, N> &a,
+                          const BArr &b_row)
+{
+#if defined(__AVX2__) && !defined(__AVX512F__) && !defined(USE_INTERLEAVED_FORMAT)
+    const double *are = a.data(), *aim = a.data() + N / 2;
+    for (uint32_t i = 0; i < N / 2; i += 4) {
+        __m256d va_re = _mm256_load_pd(are + i);
+        __m256d va_im = _mm256_load_pd(aim + i);
+        for (int m = 0; m < M; m++) {
+            const double *bre = b_row[m].data(), *bim = b_row[m].data() + N / 2;
+            double *rre = res[m].data(), *rim = res[m].data() + N / 2;
+            __m256d vb_re = _mm256_load_pd(bre + i);
+            __m256d vb_im = _mm256_load_pd(bim + i);
+            __m256d vr_re = _mm256_load_pd(rre + i);
+            __m256d vr_im = _mm256_load_pd(rim + i);
+            vr_re = _mm256_fmadd_pd(va_re, vb_re, vr_re);
+            vr_re = _mm256_fnmadd_pd(va_im, vb_im, vr_re);
+            vr_im = _mm256_fmadd_pd(va_im, vb_re, vr_im);
+            vr_im = _mm256_fmadd_pd(va_re, vb_im, vr_im);
+            _mm256_store_pd(rre + i, vr_re);
+            _mm256_store_pd(rim + i, vr_im);
+        }
+    }
+#else
+    for (int m = 0; m < M; m++)
+        FMAInFD<N>(res[m], a, b_row[m]);
+#endif
+}
+
+template <uint32_t N, int M, class ResArr, class BArr>
+inline void MulInFD_Multi(ResArr &res, const std::array<double, N> &a,
+                          const BArr &b_row)
+{
+#if defined(__AVX2__) && !defined(__AVX512F__) && !defined(USE_INTERLEAVED_FORMAT)
+    const double *are = a.data(), *aim = a.data() + N / 2;
+    for (uint32_t i = 0; i < N / 2; i += 4) {
+        __m256d va_re = _mm256_load_pd(are + i);
+        __m256d va_im = _mm256_load_pd(aim + i);
+        for (int m = 0; m < M; m++) {
+            const double *bre = b_row[m].data(), *bim = b_row[m].data() + N / 2;
+            double *rre = res[m].data(), *rim = res[m].data() + N / 2;
+            __m256d vb_re = _mm256_load_pd(bre + i);
+            __m256d vb_im = _mm256_load_pd(bim + i);
+            __m256d vr_re = _mm256_mul_pd(va_re, vb_re);
+            vr_re = _mm256_fnmadd_pd(va_im, vb_im, vr_re);
+            __m256d vr_im = _mm256_mul_pd(va_im, vb_re);
+            vr_im = _mm256_fmadd_pd(va_re, vb_im, vr_im);
+            _mm256_store_pd(rre + i, vr_re);
+            _mm256_store_pd(rim + i, vr_im);
+        }
+    }
+#else
+    for (int m = 0; m < M; m++)
+        MulInFD<N>(res[m], a, b_row[m]);
+#endif
+}
+
 template <class P>
 inline void PolyMul(Polynomial<P> &res, const Polynomial<P> &a,
                     const Polynomial<P> &b)

include/params/tfhe-rs.hpp

@@ -94,11 +94,12 @@ struct lvl1param {
         ErrorDistribution::ModularGaussian;
     // StandardDev from tfhe-rs DEFAULT_PARAMETERS glwe_noise_distribution
     static const inline double α = 9.315272083503367e-10;
-    using T = uint64_t;  // Torus representation
-    static constexpr std::make_signed_t<T> μ = 1ULL << 61;
+    using T = uint32_t;  // Torus representation (matches tfhe-rs 32-bit torus)
+    static constexpr std::make_signed_t<T> μ =
+        1U << (std::numeric_limits<T>::digits - 3);  // 1/8 = 1<<29
     static constexpr uint32_t plain_modulus = 2;
     static constexpr double Δ =
-        2 * static_cast<double>(1ULL << (std::numeric_limits<T>::digits - 1)) /
+        static_cast<double>(1ULL << std::numeric_limits<T>::digits) /
         plain_modulus;
     // Double Decomposition (bivariate representation) parameters
     // For now, set to trivial values (no actual second decomposition)