diff --git a/tachyon/base/parallelize.h b/tachyon/base/parallelize.h
index f4c4a5849c..7efc215096 100644
--- a/tachyon/base/parallelize.h
+++ b/tachyon/base/parallelize.h
@@ -233,6 +233,73 @@ auto ParallelizeMap(size_t size, Callable callback,
                                    std::move(callback));
 }
 
+template <typename Container>
+void ParallelizeFill(Container& container, typename Container::value_type value,
+                     std::optional<size_t> threshold = std::nullopt) {
+  Parallelize(
+      container,
+      [&value](absl::Span<typename Container::value_type> chunk) {
+        for (auto& v : chunk) {
+          v = value;
+        }
+      },
+      threshold);
+}
+
+template <typename Container>
+void ParallelizeResize(Container& container, size_t size,
+                       std::optional<size_t> threshold = std::nullopt) {
+  if (container.capacity() > size) {
+    container.resize(size);
+  } else {
+    std::vector<typename Container::value_type> new_container(size);
+    auto copy_span = absl::MakeSpan(new_container).first(container.size());
+    Parallelize(
+        copy_span,
+        [&container](absl::Span<typename Container::value_type> chunk,
+                     size_t chunk_offset, size_t chunk_size) {
+          size_t start = chunk_offset * chunk_size;
+          for (size_t i = 0; i < chunk.size(); ++i) {
+            chunk[i] = std::move(container[start + i]);
+          }
+        },
+        threshold);
+    container = std::move(new_container);
+  }
+}
+
+template <typename Container>
+void ParallelizeResize(Container& container, size_t size,
+                       typename Container::value_type value,
+                       std::optional<size_t> threshold = std::nullopt) {
+  if (container.capacity() > size) {
+    size_t old_size = container.size();
+    container.resize(size);
+    auto init_span = absl::MakeSpan(container).last(size - old_size);
+    Parallelize(
+        init_span,
+        [&value](absl::Span<typename Container::value_type> chunk) {
+          std::fill(chunk.begin(), chunk.end(), value);
+        },
+        threshold);
+  } else {
+    std::vector<typename Container::value_type> new_container(size);
+    Parallelize(
+        new_container,
+        [&container, &value](absl::Span<typename Container::value_type> chunk,
+                             size_t chunk_offset, size_t chunk_size) {
+          size_t start = chunk_offset * chunk_size;
+          for (size_t i = 0; i < chunk.size(); ++i) {
+            chunk[i] = (start + i) < container.size()
+                           ? std::move(container[start + i])
+                           : value;
+          }
+        },
+        threshold);
+    container = std::move(new_container);
+  }
+}
+
 }  // namespace tachyon::base
 
 #endif  // TACHYON_BASE_PARALLELIZE_H_
diff --git a/tachyon/crypto/challenger/multi_field32_conversions.h b/tachyon/crypto/challenger/multi_field32_conversions.h
index 399ae4c560..9a3e44ada2 100644
--- a/tachyon/crypto/challenger/multi_field32_conversions.h
+++ b/tachyon/crypto/challenger/multi_field32_conversions.h
@@ -26,7 +26,7 @@ BigF Reduce(absl::Span<const SmallF> values) {
   using BigInt = typename BigF::BigIntTy;
   CHECK_LT(values.size(), BigInt::kLimbNums * 2);
 
-  BigInt ret;
+  BigInt ret(0);
   for (size_t i = 0; i < values.size(); i += 2) {
     uint32_t value = values[i].value();
     if constexpr (SmallF::Config::kUseMontgomery) {
diff --git a/tachyon/crypto/commitments/fri/two_adic_multiplicative_coset.h b/tachyon/crypto/commitments/fri/two_adic_multiplicative_coset.h
index efbcd219ce..976e352a41 100644
--- a/tachyon/crypto/commitments/fri/two_adic_multiplicative_coset.h
+++ b/tachyon/crypto/commitments/fri/two_adic_multiplicative_coset.h
@@ -114,6 +114,8 @@ class TwoAdicMultiplicativeCoset {
                 domain_->group_gen_inv();
 
     size_t sz = coset.domain()->size();
+    // NOTE(batzor): These vectors are initialized below in the parallel loop so
+    // it is safe to keep it uninitialized here.
     std::vector<F> first_row(sz);
     std::vector<F> last_row(sz);
     std::vector<F> transition(sz);
diff --git a/tachyon/crypto/commitments/kzg/kzg.h b/tachyon/crypto/commitments/kzg/kzg.h
index 61637c45da..5260115eeb 100644
--- a/tachyon/crypto/commitments/kzg/kzg.h
+++ b/tachyon/crypto/commitments/kzg/kzg.h
@@ -114,6 +114,8 @@ class KZG {
 #endif
 
   void ResizeBatchCommitments(size_t size) {
+    // WARN(batzor): When resizing to a larger size, the last values will be
+    // garbage and should be filled with commitment results.
 #if TACHYON_CUDA
     if (msm_gpu_) {
       gpu_batch_commitments_.resize(size);
@@ -124,6 +126,8 @@ class KZG {
   }
 
   std::vector<Commitment> GetBatchCommitments(BatchCommitmentState& state) {
+    // NOTE(batzor): Resizing this vector without initialization is safe since
+    // |BatchNormalize| will overwrite them.
     std::vector<Commitment> batch_commitments;
 #if TACHYON_CUDA
     if (msm_gpu_) {
diff --git a/tachyon/crypto/commitments/pedersen/pedersen.h b/tachyon/crypto/commitments/pedersen/pedersen.h
index ecfec9cd78..136d109575 100644
--- a/tachyon/crypto/commitments/pedersen/pedersen.h
+++ b/tachyon/crypto/commitments/pedersen/pedersen.h
@@ -85,6 +85,8 @@ class Pedersen final
 #endif
 
   void ResizeBatchCommitments() {
+    // WARN(batzor): When resizing to a larger size, the last values will be
+    // garbage and should be filled with commitment results.
     size_t size = this->batch_commitment_state_.batch_count;
 #if TACHYON_CUDA
     if (msm_gpu_) {
@@ -96,6 +98,8 @@ class Pedersen final
   }
 
   std::vector<Commitment> GetBatchCommitments() {
+    // NOTE(batzor): Resizing this vector without initialization is safe since
+    // |BatchNormalize| will overwrite them.
     std::vector<Commitment> batch_commitments;
 #if TACHYON_CUDA
     if (msm_gpu_) {
diff --git a/tachyon/crypto/hashes/sponge/poseidon/poseidon_sponge_base.h b/tachyon/crypto/hashes/sponge/poseidon/poseidon_sponge_base.h
index 787a87e71c..e68f8b3364 100644
--- a/tachyon/crypto/hashes/sponge/poseidon/poseidon_sponge_base.h
+++ b/tachyon/crypto/hashes/sponge/poseidon/poseidon_sponge_base.h
@@ -70,7 +70,7 @@ struct PoseidonSpongeBase : public FieldBasedCryptographicSponge<Derived> {
       bytes.insert(bytes.end(), elem_bytes.begin(), elem_bytes.end());
     }
 
-    bytes.resize(num_bytes);
+    bytes.resize(num_bytes, F::Zero());
     return bytes;
   }
 
@@ -121,6 +121,8 @@ struct PoseidonSpongeBase : public FieldBasedCryptographicSponge<Derived> {
                                             size_t num_elements) const {
     const Derived& derived = static_cast<const Derived&>(*this);
 
+    // NOTE(batzor): |SqueezeInternal| will fill all the garbage values, so it
+    // is safe to have it uninitialized.
     std::vector<F> ret(num_elements);
     switch (state.mode.type) {
       case DuplexSpongeMode::Type::kAbsorbing: {
diff --git a/tachyon/math/base/BUILD.bazel b/tachyon/math/base/BUILD.bazel
index b17fc32b32..04c9a18d05 100644
--- a/tachyon/math/base/BUILD.bazel
+++ b/tachyon/math/base/BUILD.bazel
@@ -66,6 +66,11 @@ tachyon_cc_library(
     ],
 )
 
+tachyon_cc_library(
+    name = "const_init",
+    hdrs = ["const_init.h"],
+)
+
 tachyon_cc_library(
     name = "egcd",
     hdrs = ["egcd.h"],
@@ -88,6 +93,11 @@ tachyon_cc_library(
     ],
 )
 
+tachyon_cc_library(
+    name = "parallelize_threshold",
+    hdrs = ["parallelize_threshold.h"],
+)
+
 tachyon_cc_library(
     name = "rational_field",
     hdrs = ["rational_field.h"],
diff --git a/tachyon/math/base/arithmetics_results.h b/tachyon/math/base/arithmetics_results.h
index 1959d47090..ea66be5a60 100644
--- a/tachyon/math/base/arithmetics_results.h
+++ b/tachyon/math/base/arithmetics_results.h
@@ -9,8 +9,8 @@ namespace tachyon::math {
 
 template <typename T>
 struct AddResult {
-  T result{};
-  T carry{};
+  T result;
+  T carry{0};
 
   constexpr bool operator==(const AddResult& other) const {
     return result == other.result && carry == other.carry;
@@ -26,8 +26,8 @@ struct AddResult {
 
 template <typename T>
 struct SubResult {
-  T result{};
-  T borrow{};
+  T result;
+  T borrow{0};
 
   constexpr bool operator==(const SubResult& other) const {
     return result == other.result && borrow == other.borrow;
@@ -43,8 +43,8 @@ struct SubResult {
 
 template <typename T>
 struct MulResult {
-  T hi{};
-  T lo{};
+  T hi{0};
+  T lo{0};
 
   constexpr bool operator==(const MulResult& other) const {
     return hi == other.hi && lo == other.lo;
@@ -60,8 +60,8 @@ struct MulResult {
 
 template <typename T>
 struct DivResult {
-  T quotient{};
-  T remainder{};
+  T quotient;
+  T remainder = T::Zero();
 
   constexpr bool operator==(const DivResult& other) const {
     return quotient == other.quotient && remainder == other.remainder;
diff --git a/tachyon/math/base/big_int.h b/tachyon/math/base/big_int.h
index 41d8895b3c..881af7a566 100644
--- a/tachyon/math/base/big_int.h
+++ b/tachyon/math/base/big_int.h
@@ -48,9 +48,7 @@ TACHYON_EXPORT std::string LimbsToHexString(const uint64_t* limbs,
 // designed to support a wide range of big integer arithmetic operations.
 template <size_t N>
 struct BigInt {
-  uint64_t limbs[N] = {
-      0,
-  };
+  uint64_t limbs[N];
   constexpr static size_t kLimbNums = N;
   constexpr static size_t kSmallestLimbIdx = SMALLEST_INDEX(N);
   constexpr static size_t kBiggestLimbIdx = BIGGEST_INDEX(N);
@@ -61,16 +59,16 @@ struct BigInt {
 
   constexpr BigInt() = default;
   template <typename T, std::enable_if_t<std::is_signed_v<T>>* = nullptr>
-  constexpr explicit BigInt(T value) {
+  constexpr explicit BigInt(T value) : limbs{0} {
     DCHECK_GE(value, 0);
     limbs[kSmallestLimbIdx] = value;
   }
   template <typename T, std::enable_if_t<std::is_unsigned_v<T>>* = nullptr>
-  constexpr explicit BigInt(T value) {
+  constexpr explicit BigInt(T value) : limbs{0} {
     limbs[kSmallestLimbIdx] = value;
   }
   template <typename T, std::enable_if_t<std::is_signed_v<T>>* = nullptr>
-  constexpr explicit BigInt(std::initializer_list<T> values) {
+  constexpr explicit BigInt(std::initializer_list<T> values) : limbs{0} {
     DCHECK_LE(values.size(), N);
     auto it = values.begin();
     for (size_t i = 0; i < values.size(); ++i, ++it) {
@@ -79,7 +77,7 @@ struct BigInt {
     }
   }
   template <typename T, std::enable_if_t<std::is_unsigned_v<T>>* = nullptr>
-  constexpr explicit BigInt(std::initializer_list<T> values) {
+  constexpr explicit BigInt(std::initializer_list<T> values) : limbs{0} {
     DCHECK_LE(values.size(), N);
     auto it = values.begin();
     for (size_t i = 0; i < values.size(); ++i, ++it) {
@@ -117,14 +115,14 @@ struct BigInt {
 
   // Convert a decimal string to a BigInt.
   static std::optional<BigInt> FromDecString(std::string_view str) {
-    BigInt ret;
+    BigInt ret(0);
     if (!internal::StringToLimbs(str, ret.limbs, N)) return std::nullopt;
     return ret;
   }
 
   // Convert a hexadecimal string to a BigInt.
   static std::optional<BigInt> FromHexString(std::string_view str) {
-    BigInt ret;
+    BigInt ret(0);
     if (!(internal::HexStringToLimbs(str, ret.limbs, N))) return std::nullopt;
     return ret;
   }
@@ -255,8 +253,8 @@ struct BigInt {
   constexpr BigInt<N2> Extend() const {
     static_assert(N2 > N);
     BigInt<N2> ret;
-    for (size_t i = 0; i < N; ++i) {
-      ret[i] = limbs[i];
+    for (size_t i = 0; i < N2; ++i) {
+      ret[i] = i < N ? limbs[i] : 0;
     }
     return ret;
   }
@@ -707,8 +705,8 @@ struct BigInt {
       LOG_IF_NOT_GPU(ERROR) << "Division by zero attempted";
       return false;
     }
-    BigInt quotient;
-    BigInt remainder;
+    BigInt quotient(0);
+    BigInt remainder(0);
     size_t bits = BitTraits<BigInt>::GetNumBits(*this);
     uint64_t carry = 0;
     uint64_t& smallest_bit = remainder.limbs[kSmallestLimbIdx];
diff --git a/tachyon/math/base/const_init.h b/tachyon/math/base/const_init.h
new file mode 100644
index 0000000000..acc22fff6c
--- /dev/null
+++ b/tachyon/math/base/const_init.h
@@ -0,0 +1,24 @@
+#ifndef TACHYON_MATH_BASE_CONST_INIT_H_
+#define TACHYON_MATH_BASE_CONST_INIT_H_
+
+namespace tachyon::math {
+
+enum ZeroInitType {
+  kZeroInit,
+};
+
+enum OneInitType {
+  kOneInit,
+};
+
+enum MinusOneInitType {
+  kMinusOneInit,
+};
+
+enum TwoInvInitType {
+  kTwoInvInit,
+};
+
+}  // namespace tachyon::math
+
+#endif  // TACHYON_MATH_BASE_CONST_INIT_H_
diff --git a/tachyon/math/base/parallelize_threshold.h b/tachyon/math/base/parallelize_threshold.h
new file mode 100644
index 0000000000..f4e5006097
--- /dev/null
+++ b/tachyon/math/base/parallelize_threshold.h
@@ -0,0 +1,15 @@
+#ifndef TACHYON_MATH_BASE_PARALLELIZE_THRESHOLD_H_
+#define TACHYON_MATH_BASE_PARALLELIZE_THRESHOLD_H_
+
+namespace tachyon::math {
+
+struct ParallelizeThreshold {
+  // The threshold for parallelizing a loop. If the size of the loop is less
+  // than this threshold, the loop will be executed sequentially.
+  static constexpr int kFieldInit = 1e6;
+  static constexpr int kFieldSimpleOp = 1e5;
+};
+
+}  // namespace tachyon::math
+
+#endif  // TACHYON_MATH_BASE_PARALLELIZE_THRESHOLD_H_
diff --git a/tachyon/math/base/rational_field.h b/tachyon/math/base/rational_field.h
index 8e92812789..eb1bc8f41e 100644
--- a/tachyon/math/base/rational_field.h
+++ b/tachyon/math/base/rational_field.h
@@ -30,7 +30,7 @@ class RationalField : public Field<RationalField<F>> {
       : numerator_(std::move(numerator)),
         denominator_(std::move(denominator)) {}
 
-  constexpr static RationalField Zero() { return RationalField(); }
+  constexpr static RationalField Zero() { return RationalField(F::Zero()); }
 
   constexpr static RationalField One() { return RationalField(F::One()); }
 
diff --git a/tachyon/math/elliptic_curves/bn/generator/generator.cc b/tachyon/math/elliptic_curves/bn/generator/generator.cc
index f7306afe68..d25745dbea 100644
--- a/tachyon/math/elliptic_curves/bn/generator/generator.cc
+++ b/tachyon/math/elliptic_curves/bn/generator/generator.cc
@@ -13,7 +13,7 @@ namespace tachyon {
 
 template <size_t N>
 std::vector<int8_t> ComputeAteLoopCount(const mpz_class& six_x_plus_2) {
-  math::BigInt<N> x;
+  math::BigInt<N> x(0);
   math::gmp::CopyLimbs(six_x_plus_2, x.limbs);
   return x.ToNAF();
 }
diff --git a/tachyon/math/finite_fields/fp12.h b/tachyon/math/finite_fields/fp12.h
index f50557d43e..1d5d1a54ed 100644
--- a/tachyon/math/finite_fields/fp12.h
+++ b/tachyon/math/finite_fields/fp12.h
@@ -341,7 +341,7 @@ class Fp12 final : public QuadraticExtensionField<Fp12<Config>> {
     // |kFrobeniusCoeffs[0]| = q^((P⁰ - 1) / 6)
     Config::kFrobeniusCoeffs[0] = FrobeniusCoefficient::One();
 #define SET_FROBENIUS_COEFF(d)                \
-  BigInt<d * N> exp##d;                       \
+  BigInt<d * N> exp##d(0);                    \
   gmp::CopyLimbs(exp##d##_gmp, exp##d.limbs); \
   Config::kFrobeniusCoeffs[d] = BaseFieldConfig::kNonResidue.Pow(exp##d)
     // |kFrobeniusCoeffs[1]| = q^(exp₁) = q^((P¹ - 1) / 6)
diff --git a/tachyon/math/finite_fields/fp3.h b/tachyon/math/finite_fields/fp3.h
index 9f6d533210..3fd7790789 100644
--- a/tachyon/math/finite_fields/fp3.h
+++ b/tachyon/math/finite_fields/fp3.h
@@ -81,7 +81,7 @@ class Fp3 final : public CubicExtensionField<Fp3<Config>> {
     // |kFrobeniusCoeffs[0]| = q^((P⁰ - 1) / 3) = 1
     Config::kFrobeniusCoeffs[0] = FrobeniusCoefficient::One();
 #define SET_FROBENIUS_COEFF(d)                \
-  BigInt<d * N> exp##d;                       \
+  BigInt<d * N> exp##d(0);                    \
   gmp::CopyLimbs(exp##d##_gmp, exp##d.limbs); \
   Config::kFrobeniusCoeffs[d] = Config::kNonResidue.Pow(exp##d)
 
diff --git a/tachyon/math/finite_fields/fp4.h b/tachyon/math/finite_fields/fp4.h
index 45e52b2ec3..82e3508f55 100644
--- a/tachyon/math/finite_fields/fp4.h
+++ b/tachyon/math/finite_fields/fp4.h
@@ -89,7 +89,7 @@ class Fp4<Config, std::enable_if_t<Config::kDegreeOverBaseField == 2>> final
     // |kFrobeniusCoeffs[0]| = q^((P⁰ - 1) / 4) = 1
     Config::kFrobeniusCoeffs[0] = FrobeniusCoefficient::One();
 #define SET_FROBENIUS_COEFF(d)                \
-  BigInt<d * N> exp##d;                       \
+  BigInt<d * N> exp##d(0);                    \
   gmp::CopyLimbs(exp##d##_gmp, exp##d.limbs); \
   Config::kFrobeniusCoeffs[d] = BaseFieldConfig::kNonResidue.Pow(exp##d)
 
@@ -176,7 +176,7 @@ class Fp4<Config, std::enable_if_t<Config::kDegreeOverBaseField == 4>> final
     // |kFrobeniusCoeffs[0]| = q^((P⁰ - 1) / 4) = 1
     Config::kFrobeniusCoeffs[0] = FrobeniusCoefficient::One();
 #define SET_FROBENIUS_COEFF(d)                \
-  BigInt<d * N> exp##d;                       \
+  BigInt<d * N> exp##d(0);                    \
   gmp::CopyLimbs(exp##d##_gmp, exp##d.limbs); \
   Config::kFrobeniusCoeffs[d] = Config::kNonResidue.Pow(exp##d)
 
diff --git a/tachyon/math/finite_fields/fp6.h b/tachyon/math/finite_fields/fp6.h
index 3dce0e49a9..3d640ec393 100644
--- a/tachyon/math/finite_fields/fp6.h
+++ b/tachyon/math/finite_fields/fp6.h
@@ -203,7 +203,7 @@ class Fp6<Config, std::enable_if_t<Config::kDegreeOverBaseField == 2>> final
     // |kFrobeniusCoeffs[0]| = q^((P⁰ - 1) / 6) = 1
     Config::kFrobeniusCoeffs[0] = FrobeniusCoefficient::One();
 #define SET_FROBENIUS_COEFF(d)                \
-  BigInt<d * N> exp##d;                       \
+  BigInt<d * N> exp##d(0);                    \
   gmp::CopyLimbs(exp##d##_gmp, exp##d.limbs); \
   Config::kFrobeniusCoeffs[d] = BaseFieldConfig::kNonResidue.Pow(exp##d)
 
@@ -365,7 +365,7 @@ class Fp6<Config, std::enable_if_t<Config::kDegreeOverBaseField == 3>> final
     // |kFrobeniusCoeffs[0]| = q^((P⁰ - 1) / 3)
     Config::kFrobeniusCoeffs[0] = FrobeniusCoefficient::One();
 #define SET_FROBENIUS_COEFF(d)                \
-  BigInt<d * N> exp##d;                       \
+  BigInt<d * N> exp##d(0);                    \
   gmp::CopyLimbs(exp##d##_gmp, exp##d.limbs); \
   Config::kFrobeniusCoeffs[d] = Config::kNonResidue.Pow(exp##d)
 
diff --git a/tachyon/math/finite_fields/generator/generator_util.h b/tachyon/math/finite_fields/generator/generator_util.h
index ee1c791743..6a0ee4ed8e 100644
--- a/tachyon/math/finite_fields/generator/generator_util.h
+++ b/tachyon/math/finite_fields/generator/generator_util.h
@@ -18,7 +18,7 @@ std::string MpzClassToString(const mpz_class& m);
 
 template <size_t N>
 std::string MpzClassToMontString(const mpz_class& v_in, const mpz_class& m_in) {
-  BigInt<N> m;
+  BigInt<N> m(0);
 
   if constexpr (N == 1) {
     // NOLINTNEXTLINE(runtime/int)
@@ -31,7 +31,7 @@ std::string MpzClassToMontString(const mpz_class& v_in, const mpz_class& m_in) {
   BigInt<N> r2 = Modulus<N>::MontgomeryR2(m);
   uint64_t inv = Modulus<N>::template Inverse<uint64_t>(m);
 
-  BigInt<N> v;
+  BigInt<N> v(0);
   if (gmp::IsNegative(v_in)) {
     gmp::CopyLimbs(m_in + v_in, v.limbs);
   } else {
diff --git a/tachyon/math/finite_fields/generator/prime_field_generator/BUILD.bazel b/tachyon/math/finite_fields/generator/prime_field_generator/BUILD.bazel
index d5f7b028be..e1068db663 100644
--- a/tachyon/math/finite_fields/generator/prime_field_generator/BUILD.bazel
+++ b/tachyon/math/finite_fields/generator/prime_field_generator/BUILD.bazel
@@ -26,6 +26,7 @@ tachyon_cc_binary(
         "//tachyon/base/flag:flag_parser",
         "//tachyon/build:cc_writer",
         "//tachyon/math/base:bit_iterator",
+        "//tachyon/math/base:const_init",
         "//tachyon/math/base/gmp:bit_traits",
         "//tachyon/math/finite_fields:prime_field_util",
         "//tachyon/math/finite_fields/generator:generator_util",
diff --git a/tachyon/math/finite_fields/generator/prime_field_generator/build_defs.bzl b/tachyon/math/finite_fields/generator/prime_field_generator/build_defs.bzl
index 358d5afb4f..6f363bb4f9 100644
--- a/tachyon/math/finite_fields/generator/prime_field_generator/build_defs.bzl
+++ b/tachyon/math/finite_fields/generator/prime_field_generator/build_defs.bzl
@@ -200,6 +200,7 @@ def _do_generate_prime_fields(
         "//tachyon:export",
         "//tachyon/build:build_config",
         "//tachyon/math/base:big_int",
+        "//tachyon/math/base:const_init",
         "@com_google_absl//absl/base",
     ]
 
diff --git a/tachyon/math/finite_fields/generator/prime_field_generator/prime_field_generator.cc b/tachyon/math/finite_fields/generator/prime_field_generator/prime_field_generator.cc
index 1184db899e..7400ef27d9 100644
--- a/tachyon/math/finite_fields/generator/prime_field_generator/prime_field_generator.cc
+++ b/tachyon/math/finite_fields/generator/prime_field_generator/prime_field_generator.cc
@@ -50,7 +50,7 @@ struct ModulusInfo {
 
   template <size_t N>
   static ModulusInfo From(const mpz_class& m_in) {
-    math::BigInt<N> m;
+    math::BigInt<N> m(0);
     math::gmp::CopyLimbs(m_in, m.limbs);
 
     ModulusInfo ret;
diff --git a/tachyon/math/finite_fields/generator/prime_field_generator/prime_field_x86.h.tpl b/tachyon/math/finite_fields/generator/prime_field_generator/prime_field_x86.h.tpl
index d9fe90394e..acafe6a243 100644
--- a/tachyon/math/finite_fields/generator/prime_field_generator/prime_field_x86.h.tpl
+++ b/tachyon/math/finite_fields/generator/prime_field_generator/prime_field_x86.h.tpl
@@ -9,6 +9,7 @@
 #include "absl/base/call_once.h"
 
 #include "tachyon/math/base/big_int.h"
+#include "tachyon/math/base/const_init.h"
 #include "tachyon/math/base/gmp/gmp_util.h"
 #include "tachyon/math/finite_fields/prime_field_base.h"
 #include "tachyon/math/base/byinverter.h"
@@ -42,11 +43,11 @@ class PrimeField<_Config, std::enable_if_t<_Config::%{asm_flag}>> final
 
   using CpuField = PrimeField<Config>;
   using GpuField = PrimeFieldGpu<Config>;
-  
+
   constexpr static BYInverter<N> inverter =
       BYInverter<N>(Config::kModulus, Config::kMontgomeryR2);
 
-  constexpr PrimeField() = default;
+  PrimeField() {}
   template <typename T,
             std::enable_if_t<std::is_constructible_v<BigInt<N>, T>>* = nullptr>
   constexpr explicit PrimeField(T value) : PrimeField(BigInt<N>(value)) {}
@@ -54,30 +55,23 @@ class PrimeField<_Config, std::enable_if_t<_Config::%{asm_flag}>> final
     DCHECK_LT(value, Config::kModulus);
     %{prefix}_rawToMontgomery(value_.limbs, value.limbs);
   }
+  constexpr explicit PrimeField(ZeroInitType zero_init) : value_(0) {}
+  constexpr explicit PrimeField(OneInitType one_init) : value_(Config::kOne) {}
+  constexpr explicit PrimeField(MinusOneInitType minus_one_init) : value_(Config::kMinusOne) {}
+  constexpr explicit PrimeField(TwoInvInitType two_inv_init) : value_(Config::kTwoInv) {}
+
   constexpr PrimeField(const PrimeField& other) = default;
   constexpr PrimeField& operator=(const PrimeField& other) = default;
   constexpr PrimeField(PrimeField&& other) = default;
   constexpr PrimeField& operator=(PrimeField&& other) = default;
 
-  constexpr static PrimeField Zero() { return PrimeField(); }
+  constexpr static PrimeField Zero() { return PrimeField(kZeroInit); }
 
-  constexpr static PrimeField One() {
-    PrimeField ret{};
-    ret.value_ = Config::kOne;
-    return ret;
-  }
+  constexpr static PrimeField One() { return PrimeField(kOneInit); }
 
-  constexpr static PrimeField MinusOne() {
-    PrimeField ret{};
-    ret.value_ = Config::kMinusOne;
-    return ret;
-  }
+  constexpr static PrimeField MinusOne() { return PrimeField(kMinusOneInit); }
 
-  constexpr static PrimeField TwoInv() {
-    PrimeField ret{};
-    ret.value_ = Config::kTwoInv;
-    return ret;
-  }
+  constexpr static PrimeField TwoInv() { return PrimeField(kTwoInvInit); }
 
   static PrimeField Random() {
     return PrimeField(BigInt<N>::Random(Config::kModulus));
@@ -113,7 +107,7 @@ class PrimeField<_Config, std::enable_if_t<_Config::%{asm_flag}>> final
   }
 
   static PrimeField FromMpzClass(const mpz_class& value) {
-    BigInt<N> big_int;
+    BigInt<N> big_int(0);
     gmp::CopyLimbs(value, big_int.limbs);
     return FromBigInt(big_int);
   }
diff --git a/tachyon/math/finite_fields/modulus.h b/tachyon/math/finite_fields/modulus.h
index 1bc1758078..c1c75e421e 100644
--- a/tachyon/math/finite_fields/modulus.h
+++ b/tachyon/math/finite_fields/modulus.h
@@ -43,21 +43,21 @@ class Modulus {
   }
 
   constexpr static BigInt<N> MontgomeryR(const BigInt<N>& modulus) {
-    BigInt<N + 1> two_pow_n_times_64;
+    BigInt<N + 1> two_pow_n_times_64(0);
     two_pow_n_times_64.biggest_limb() = uint64_t{1};
     return (two_pow_n_times_64 % modulus.template Extend<N + 1>())
         .template Shrink<N>();
   }
 
   constexpr static BigInt<N> MontgomeryR2(const BigInt<N>& modulus) {
-    BigInt<2 * N + 1> two_pow_n_times_64_square;
+    BigInt<2 * N + 1> two_pow_n_times_64_square(0);
     two_pow_n_times_64_square.biggest_limb() = uint64_t{1};
     return (two_pow_n_times_64_square % modulus.template Extend<2 * N + 1>())
         .template Shrink<N>();
   }
 
   constexpr static BigInt<N> MontgomeryR3(const BigInt<N>& modulus) {
-    BigInt<3 * N + 1> two_pow_n_times_64_cubic;
+    BigInt<3 * N + 1> two_pow_n_times_64_cubic(0);
     two_pow_n_times_64_cubic.biggest_limb() = uint64_t{1};
     return (two_pow_n_times_64_cubic % modulus.template Extend<3 * N + 1>())
         .template Shrink<N>();
diff --git a/tachyon/math/finite_fields/prime_field_fallback.h b/tachyon/math/finite_fields/prime_field_fallback.h
index febb7e0362..786957fcbb 100644
--- a/tachyon/math/finite_fields/prime_field_fallback.h
+++ b/tachyon/math/finite_fields/prime_field_fallback.h
@@ -119,7 +119,7 @@ class PrimeField<_Config, std::enable_if_t<!_Config::kUseAsm &&
   }
 
   static PrimeField FromMpzClass(const mpz_class& value) {
-    BigInt<N> big_int;
+    BigInt<N> big_int(0);
     gmp::CopyLimbs(value, big_int.limbs);
     return FromBigInt(big_int);
   }
@@ -330,7 +330,7 @@ class PrimeField<_Config, std::enable_if_t<!_Config::kUseAsm &&
 
   constexpr static void DoFastMul(const PrimeField& a, const PrimeField& b,
                                   PrimeField& c) {
-    BigInt<N> r;
+    BigInt<N> r(0);
     for (size_t i = 0; i < N; ++i) {
       MulResult<uint64_t> result;
       result = internal::u64::MulAddWithCarry(r[0], a[0], b[i]);
@@ -362,7 +362,7 @@ class PrimeField<_Config, std::enable_if_t<!_Config::kUseAsm &&
   }
 
   constexpr static void DoSquareImpl(const PrimeField& a, PrimeField& b) {
-    BigInt<2 * N> r;
+    BigInt<2 * N> r(0);
     MulResult<uint64_t> mul_result;
     for (size_t i = 0; i < N - 1; ++i) {
       for (size_t j = i + 1; j < N; ++j) {
diff --git a/tachyon/math/finite_fields/prime_field_gpu.h b/tachyon/math/finite_fields/prime_field_gpu.h
index 73d2824787..a5f2dcfb9d 100644
--- a/tachyon/math/finite_fields/prime_field_gpu.h
+++ b/tachyon/math/finite_fields/prime_field_gpu.h
@@ -103,7 +103,7 @@ class PrimeFieldGpu final : public PrimeFieldBase<PrimeFieldGpu<_Config>> {
   }
 
   static PrimeFieldGpu FromMpzClass(const mpz_class& value) {
-    BigInt<N> big_int;
+    BigInt<N> big_int(0);
     gmp::CopyLimbs(value, big_int.limbs);
     return FromBigInt(big_int);
   }
diff --git a/tachyon/math/finite_fields/prime_field_gpu_debug.h b/tachyon/math/finite_fields/prime_field_gpu_debug.h
index 9d003f5c12..bdc3ef8778 100644
--- a/tachyon/math/finite_fields/prime_field_gpu_debug.h
+++ b/tachyon/math/finite_fields/prime_field_gpu_debug.h
@@ -112,7 +112,7 @@ class PrimeFieldGpuDebug final
   }
 
   static PrimeFieldGpuDebug FromMpzClass(const mpz_class& value) {
-    BigInt<N> big_int;
+    BigInt<N> big_int(0);
     gmp::CopyLimbs(value, big_int.limbs);
     return FromBigInt(big_int);
   }
diff --git a/tachyon/math/polynomials/multivariate/BUILD.bazel b/tachyon/math/polynomials/multivariate/BUILD.bazel
index 94cc6f2f2a..ed0e48f4e9 100644
--- a/tachyon/math/polynomials/multivariate/BUILD.bazel
+++ b/tachyon/math/polynomials/multivariate/BUILD.bazel
@@ -35,9 +35,11 @@ tachyon_cc_library(
         "//tachyon/base:bits",
         "//tachyon/base:logging",
         "//tachyon/base:openmp_util",
+        "//tachyon/base:parallelize",
         "//tachyon/base/buffer:copyable",
         "//tachyon/base/containers:container_util",
         "//tachyon/base/strings:string_util",
+        "//tachyon/math/base:parallelize_threshold",
         "//tachyon/math/polynomials:polynomial",
     ],
 )
diff --git a/tachyon/math/polynomials/multivariate/multilinear_dense_evaluations.h b/tachyon/math/polynomials/multivariate/multilinear_dense_evaluations.h
index eb712fb0da..fda7c47942 100644
--- a/tachyon/math/polynomials/multivariate/multilinear_dense_evaluations.h
+++ b/tachyon/math/polynomials/multivariate/multilinear_dense_evaluations.h
@@ -18,7 +18,9 @@
 #include "tachyon/base/buffer/copyable.h"
 #include "tachyon/base/containers/container_util.h"
 #include "tachyon/base/logging.h"
+#include "tachyon/base/parallelize.h"
 #include "tachyon/base/strings/string_util.h"
+#include "tachyon/math/base/parallelize_threshold.h"
 #include "tachyon/math/polynomials/multivariate/support_poly_operators.h"
 
 namespace tachyon {
@@ -149,8 +151,10 @@ class MultilinearDenseEvaluations {
   // NOTE(chokobole): This creates a polynomial that contains |F::Zero()| up to
   // |degree| + 1.
   constexpr static MultilinearDenseEvaluations Zero(size_t degree) {
-    MultilinearDenseEvaluations ret{};
-    ret.evaluations_ = std::vector<F>(size_t{1} << degree);
+    MultilinearDenseEvaluations ret;
+    ret.evaluations_.resize(size_t{1} << degree);
+    base::ParallelizeFill(ret.evaluations_, F::Zero(),
+                          /*threshold=*/ParallelizeThreshold::kFieldInit);
     return ret;
   }
 
diff --git a/tachyon/math/polynomials/univariate/BUILD.bazel b/tachyon/math/polynomials/univariate/BUILD.bazel
index eb01ae5643..9e00d127d4 100644
--- a/tachyon/math/polynomials/univariate/BUILD.bazel
+++ b/tachyon/math/polynomials/univariate/BUILD.bazel
@@ -158,9 +158,11 @@ tachyon_cc_library(
         "//tachyon/base:logging",
         "//tachyon/base:openmp_util",
         "//tachyon/base:optional",
+        "//tachyon/base:parallelize",
         "//tachyon/base/buffer:copyable",
         "//tachyon/base/containers:container_util",
         "//tachyon/base/json",
+        "//tachyon/math/base:parallelize_threshold",
         "//tachyon/math/polynomials:polynomial",
     ],
 )
@@ -188,6 +190,7 @@ tachyon_cc_library(
         "//tachyon/base/ranges:algorithm",
         "//tachyon/base/strings:string_util",
         "//tachyon/math/base:arithmetics_results",
+        "//tachyon/math/base:parallelize_threshold",
         "//tachyon/math/polynomials:polynomial",
         "@com_google_absl//absl/hash",
         "@com_google_absl//absl/numeric:bits",
diff --git a/tachyon/math/polynomials/univariate/lagrange_interpolation.h b/tachyon/math/polynomials/univariate/lagrange_interpolation.h
index 9fe6805c52..098b9ca401 100644
--- a/tachyon/math/polynomials/univariate/lagrange_interpolation.h
+++ b/tachyon/math/polynomials/univariate/lagrange_interpolation.h
@@ -63,7 +63,10 @@ bool LagrangeInterpolate(
         // clang-format on
 
         std::vector<F> coeffs(points.size());
-        std::vector<F> coeffs_sum(points.size());
+
+        // NOTE(batzor): We initialize |coeffs_sum| serially since it is already
+        // inside a parallel loop.
+        std::vector<F> coeffs_sum(points.size(), F::Zero());
 
         size_t start = chunk_offset * chunk_size;
         for (size_t chunk_idx = 0; chunk_idx < chunk.size(); ++chunk_idx) {
diff --git a/tachyon/math/polynomials/univariate/univariate_dense_coefficients.h b/tachyon/math/polynomials/univariate/univariate_dense_coefficients.h
index 06a5c5d0aa..0f606846e4 100644
--- a/tachyon/math/polynomials/univariate/univariate_dense_coefficients.h
+++ b/tachyon/math/polynomials/univariate/univariate_dense_coefficients.h
@@ -22,6 +22,7 @@
 #include "tachyon/base/containers/container_util.h"
 #include "tachyon/base/parallelize.h"
 #include "tachyon/base/strings/string_util.h"
+#include "tachyon/math/base/parallelize_threshold.h"
 #include "tachyon/math/polynomials/univariate/support_poly_operators.h"
 #include "tachyon/math/polynomials/univariate/univariate_evaluation_domain_forwards.h"
 
@@ -108,6 +109,8 @@ class UnivariateDenseCoefficients {
     // clang-format on
 
     std::vector<F> coefficients(std::size(roots) + 1);
+    base::ParallelizeFill(coefficients, F::Zero(),
+                          /*threshold=*/ParallelizeThreshold::kFieldInit);
     coefficients[0] = F::One();
     for (size_t i = 0; i < std::size(roots); ++i) {
       for (size_t j = i + 1; j > 0; --j) {
@@ -179,6 +182,8 @@ class UnivariateDenseCoefficients {
   Fold(const Field& r) const {
     size_t size = coefficients_.size();
     std::vector<F> coefficients((size + 1) >> 1);
+    base::ParallelizeFill(coefficients, F::Zero(),
+                          /*threshold=*/ParallelizeThreshold::kFieldInit);
     OMP_PARALLEL_FOR(size_t i = 0; i < size; i += 2) {
       coefficients[i >> 1] = coefficients_[i + 1] * r;
       coefficients[i >> 1] += coefficients_[i];
@@ -242,7 +247,9 @@ class UnivariateDenseCoefficients {
   // with |IsZero()|, it returns false. So please use it carefully!
   constexpr static UnivariateDenseCoefficients Zero(size_t degree) {
     UnivariateDenseCoefficients ret;
-    ret.coefficients_ = std::vector<F>(degree + 1);
+    ret.coefficients_.resize(degree + 1);
+    base::ParallelizeFill(ret.coefficients_, F::Zero(),
+                          /*threshold=*/ParallelizeThreshold::kFieldInit);
     return ret;
   }
 
diff --git a/tachyon/math/polynomials/univariate/univariate_evaluations.h b/tachyon/math/polynomials/univariate/univariate_evaluations.h
index 4fd51d1969..6955448253 100644
--- a/tachyon/math/polynomials/univariate/univariate_evaluations.h
+++ b/tachyon/math/polynomials/univariate/univariate_evaluations.h
@@ -18,7 +18,9 @@
 #include "tachyon/base/containers/container_util.h"
 #include "tachyon/base/json/json.h"
 #include "tachyon/base/logging.h"
+#include "tachyon/base/parallelize.h"
 #include "tachyon/base/strings/string_util.h"
+#include "tachyon/math/base/parallelize_threshold.h"
 #include "tachyon/math/polynomials/polynomial.h"
 #include "tachyon/math/polynomials/univariate/univariate_evaluation_domain_forwards.h"
 
@@ -250,7 +252,9 @@ class UnivariateEvaluations final
   // |degree| + 1.
   constexpr static UnivariateEvaluations Zero(size_t degree) {
     UnivariateEvaluations ret;
-    ret.evaluations_ = std::vector<F>(degree + 1);
+    ret.evaluations_.resize(degree + 1);
+    base::ParallelizeFill(ret.evaluations_, F::Zero(),
+                          /*threshold=*/ParallelizeThreshold::kFieldInit);
     return ret;
   }
 
diff --git a/tachyon/math/polynomials/univariate/univariate_polynomial_ops.h b/tachyon/math/polynomials/univariate/univariate_polynomial_ops.h
index 5369811161..ed8e669891 100644
--- a/tachyon/math/polynomials/univariate/univariate_polynomial_ops.h
+++ b/tachyon/math/polynomials/univariate/univariate_polynomial_ops.h
@@ -15,8 +15,10 @@
 #include "tachyon/base/containers/container_util.h"
 #include "tachyon/base/openmp_util.h"
 #include "tachyon/base/optional.h"
+#include "tachyon/base/parallelize.h"
 #include "tachyon/base/sort.h"
 #include "tachyon/math/base/arithmetics_results.h"
+#include "tachyon/math/base/parallelize_threshold.h"
 #include "tachyon/math/polynomials/univariate/univariate_polynomial.h"
 
 namespace tachyon::math {
@@ -61,8 +63,11 @@ class UnivariatePolynomialOp<UnivariateDenseCoefficients<F, MaxDegree>> {
 
     std::vector<F>& l_coefficients = self.coefficients_.coefficients_;
     const std::vector<F>& r_coefficients = other.coefficients_.coefficients_;
-    l_coefficients.resize(
-        std::max(l_coefficients.size(), r_coefficients.size()));
+    base::ParallelizeResize(
+        l_coefficients, std::max(l_coefficients.size(), r_coefficients.size()),
+        F::Zero(),
+        /*threshold=*/ParallelizeThreshold::kFieldInit);
+
     OMP_PARALLEL_FOR(size_t i = 0; i < r_coefficients.size(); ++i) {
       l_coefficients[i] += r_coefficients[i];
     }
@@ -104,7 +109,8 @@ class UnivariatePolynomialOp<UnivariateDenseCoefficients<F, MaxDegree>> {
     size_t other_degree = other.Degree();
     std::vector<F> upper_coeffs;
     if (degree < other_degree) {
-      upper_coeffs = std::vector<F>(other_degree - degree);
+      base::ParallelizeResize(upper_coeffs, other_degree - degree, F::Zero(),
+                              /*threshold=*/ParallelizeThreshold::kFieldInit);
     }
 
     std::vector<F>& l_coefficients = self.coefficients_.coefficients_;
@@ -154,8 +160,10 @@ class UnivariatePolynomialOp<UnivariateDenseCoefficients<F, MaxDegree>> {
 
     std::vector<F>& l_coefficients = self.coefficients_.coefficients_;
     const std::vector<F>& r_coefficients = other.coefficients_.coefficients_;
-    l_coefficients.resize(
-        std::max(l_coefficients.size(), r_coefficients.size()));
+    base::ParallelizeResize(
+        l_coefficients, std::max(l_coefficients.size(), r_coefficients.size()),
+        F::Zero(),
+        /*threshold=*/ParallelizeThreshold::kFieldInit);
     OMP_PARALLEL_FOR(size_t i = 0; i < r_coefficients.size(); ++i) {
       l_coefficients[i] -= r_coefficients[i];
     }
@@ -439,7 +447,9 @@ class UnivariatePolynomialOp<UnivariateDenseCoefficients<F, MaxDegree>> {
   static UnivariatePolynomial<D>& Copy(UnivariatePolynomial<D>& self,
                                        const UnivariatePolynomial<S>& other) {
     std::vector<F>& l_coefficients = self.coefficients_.coefficients_;
-    l_coefficients = std::vector<F>(other.Degree() + 1);
+    l_coefficients.resize(other.Degree() + 1);
+    base::ParallelizeFill(l_coefficients, F::Zero(),
+                          /*threshold=*/ParallelizeThreshold::kFieldInit);
 
     const std::vector<Term>& r_terms = other.coefficients().terms_;
     OMP_PARALLEL_FOR(const Term& r_term : r_terms) {
@@ -462,6 +472,8 @@ class UnivariatePolynomialOp<UnivariateDenseCoefficients<F, MaxDegree>> {
     const std::vector<F>& a_coefficients = a.coefficients_.coefficients_;
     const std::vector<F>& b_coefficients = b.coefficients_.coefficients_;
     std::vector<F> c_coefficients(a_degree + b_degree + 1);
+    base::ParallelizeFill(c_coefficients, F::Zero(),
+                          /*threshold=*/ParallelizeThreshold::kFieldInit);
     for (size_t i = 0; i < b_coefficients.size(); ++i) {
       const F& b = b_coefficients[i];
       if (b.IsZero()) {
@@ -486,6 +498,8 @@ class UnivariatePolynomialOp<UnivariateDenseCoefficients<F, MaxDegree>> {
     const std::vector<F>& a_coefficients = a.coefficients_.coefficients_;
     const std::vector<Term>& b_terms = b.coefficients().terms_;
     std::vector<F> c_coefficients(a_degree + b_degree + 1);
+    base::ParallelizeFill(c_coefficients, F::Zero(),
+                          /*threshold=*/ParallelizeThreshold::kFieldInit);
     for (size_t i = 0; i < b_terms.size(); ++i) {
       const F& b = b_terms[i].coefficient;
       if (b.IsZero()) {
@@ -520,6 +534,8 @@ class UnivariatePolynomialOp<UnivariateDenseCoefficients<F, MaxDegree>> {
       return true;
     }
     std::vector<F> quotient(self.Degree() - other.Degree() + 1);
+    base::ParallelizeFill(quotient, F::Zero(),
+                          /*threshold=*/ParallelizeThreshold::kFieldInit);
     UnivariatePolynomial<D> remainder = self.ToDense();
     std::vector<F>& r_coefficients = remainder.coefficients_.coefficients_;
     F divisor_leading_inv = *other.GetLeadingCoefficient().Inverse();
diff --git a/tachyon/zk/lookup/halo2/BUILD.bazel b/tachyon/zk/lookup/halo2/BUILD.bazel
index 39f29e169f..f7aff30ee9 100644
--- a/tachyon/zk/lookup/halo2/BUILD.bazel
+++ b/tachyon/zk/lookup/halo2/BUILD.bazel
@@ -22,7 +22,9 @@ tachyon_cc_library(
     name = "permute_expression_pair",
     hdrs = ["permute_expression_pair.h"],
     deps = [
+        "//tachyon/base:parallelize",
         "//tachyon/base:sort",
+        "//tachyon/math/base:parallelize_threshold",
         "//tachyon/zk/base/entities:prover_base",
         "//tachyon/zk/lookup:pair",
         "@com_google_absl//absl/container:btree",
diff --git a/tachyon/zk/lookup/halo2/permute_expression_pair.h b/tachyon/zk/lookup/halo2/permute_expression_pair.h
index 5e75937d9a..79e87cd737 100644
--- a/tachyon/zk/lookup/halo2/permute_expression_pair.h
+++ b/tachyon/zk/lookup/halo2/permute_expression_pair.h
@@ -13,7 +13,9 @@
 
 #include "absl/container/btree_map.h"
 
+#include "tachyon/base/parallelize.h"
 #include "tachyon/base/sort.h"
+#include "tachyon/math/base/parallelize_threshold.h"
 #include "tachyon/zk/base/entities/prover_base.h"
 #include "tachyon/zk/lookup/pair.h"
 
@@ -53,6 +55,11 @@ template <typename PCS, typename Evals, typename F = typename Evals::Field>
   }
 
   std::vector<F> permuted_table_expressions(domain_size);
+  // NOTE(batzor): First |usable_rows| elements will be initialized afterwards.
+  auto init_span = absl::MakeSpan(permuted_table_expressions)
+                       .last(domain_size - usable_rows);
+  base::ParallelizeFill(init_span, F::Zero(),
+                        /*threshold=*/math::ParallelizeThreshold::kFieldInit);
 
   std::vector<RowIndex> repeated_input_rows;
   repeated_input_rows.reserve(usable_rows - 1);
diff --git a/tachyon/zk/lookup/log_derivative_halo2/prover.h b/tachyon/zk/lookup/log_derivative_halo2/prover.h
index e5aa4a85f3..69fb2eefb0 100644
--- a/tachyon/zk/lookup/log_derivative_halo2/prover.h
+++ b/tachyon/zk/lookup/log_derivative_halo2/prover.h
@@ -58,6 +58,8 @@ struct GrandSumPolysTempStorage {
   std::vector<F> inputs_log_derivatives;
   std::vector<F> table_log_derivatives;
 
+  // NOTE(batzor): Constructing these vectors without initialization is safe
+  // since |CreateGrandSumPoly| will overwrite them.
   explicit GrandSumPolysTempStorage(size_t usable_rows)
       : inputs_log_derivatives(usable_rows),
         table_log_derivatives(usable_rows) {}
diff --git a/tachyon/zk/lookup/log_derivative_halo2/prover_impl.h b/tachyon/zk/lookup/log_derivative_halo2/prover_impl.h
index 3c6e3c13ab..f981736cc8 100644
--- a/tachyon/zk/lookup/log_derivative_halo2/prover_impl.h
+++ b/tachyon/zk/lookup/log_derivative_halo2/prover_impl.h
@@ -112,6 +112,8 @@ BlindedPolynomial<Poly, Evals> Prover<Poly, Evals>::ComputeMPoly(
   base::StableSort(storage.sorted_table_with_indices.begin(),
                    storage.sorted_table_with_indices.end());
 
+  // NOTE(batzor): This vector is initialized below in the parallel loop so
+  // it is safe to keep it uninitialized here.
   std::vector<F> m_values(prover->pcs().N());
   OMP_PARALLEL {
     OMP_NESTED_FOR(size_t i = 0; i < compressed_inputs.size(); ++i) {
@@ -128,9 +130,10 @@ BlindedPolynomial<Poly, Evals> Prover<Poly, Evals>::ComputeMPoly(
     }
 
     // Convert atomic |m_values| to |Evals|.
-    OMP_FOR(RowIndex i = 0; i < usable_rows; ++i) {
-      m_values[i] =
-          F(storage.m_values_atomic[i].exchange(0, std::memory_order_relaxed));
+    OMP_FOR(RowIndex i = 0; i < m_values.size(); ++i) {
+      m_values[i] = i < usable_rows ? F(storage.m_values_atomic[i].exchange(
+                                          0, std::memory_order_relaxed))
+                                    : F::Zero();
     }
   }
 
@@ -232,6 +235,8 @@ BlindedPolynomial<Poly, Evals> Prover<Poly, Evals>::CreateGrandSumPoly(
   // 1 / τ(X)
   ComputeLogDerivatives(compressed_table, beta, storage.table_log_derivatives);
 
+  // NOTE(batzor): This vector is initialized below in the parallel loop so
+  // it is safe to keep it uninitialized here.
   std::vector<F> grand_sum(n);
   grand_sum[0] = F::Zero();
 
@@ -240,10 +245,14 @@ BlindedPolynomial<Poly, Evals> Prover<Poly, Evals>::CreateGrandSumPoly(
   // |storage.inputs_log_derivatives| since the current values of
   // |storage.inputs_log_derivatives| are not needed anymore.
   std::vector<F>& log_derivatives_diff = storage.inputs_log_derivatives;
-  OMP_PARALLEL_FOR(size_t i = 0; i < usable_rows; ++i) {
-    log_derivatives_diff[i] -= m_values[i] * storage.table_log_derivatives[i];
-    if (i != usable_rows - 1) {
-      grand_sum[i + 1] = log_derivatives_diff[i];
+  OMP_PARALLEL_FOR(size_t i = 0; i < n; ++i) {
+    if (i < usable_rows) {
+      log_derivatives_diff[i] -= m_values[i] * storage.table_log_derivatives[i];
+      if (i != usable_rows - 1) {
+        grand_sum[i + 1] = log_derivatives_diff[i];
+      }
+    } else {
+      grand_sum[i] = F::Zero();
     }
   }
 
diff --git a/tachyon/zk/plonk/constraint_system/BUILD.bazel b/tachyon/zk/plonk/constraint_system/BUILD.bazel
index 7cb876b2af..75dfe74f26 100644
--- a/tachyon/zk/plonk/constraint_system/BUILD.bazel
+++ b/tachyon/zk/plonk/constraint_system/BUILD.bazel
@@ -122,9 +122,11 @@ tachyon_cc_library(
         ":exclusion_matrix",
         ":selector_assignment",
         ":selector_description",
+        "//tachyon/base:parallelize",
         "//tachyon/base/containers:container_util",
         "//tachyon/base/containers:cxx20_erase",
         "//tachyon/base/functional:callback",
+        "//tachyon/math/base:parallelize_threshold",
         "//tachyon/zk/expressions:expression_factory",
         "@com_google_googletest//:gtest_prod",
     ],
diff --git a/tachyon/zk/plonk/constraint_system/selector_compressor.h b/tachyon/zk/plonk/constraint_system/selector_compressor.h
index 033f09775a..923e744fe5 100644
--- a/tachyon/zk/plonk/constraint_system/selector_compressor.h
+++ b/tachyon/zk/plonk/constraint_system/selector_compressor.h
@@ -20,6 +20,8 @@
 #include "tachyon/base/containers/container_util.h"
 #include "tachyon/base/containers/cxx20_erase_vector.h"
 #include "tachyon/base/functional/callback.h"
+#include "tachyon/base/parallelize.h"
+#include "tachyon/math/base/parallelize_threshold.h"
 #include "tachyon/zk/expressions/expression_factory.h"
 #include "tachyon/zk/plonk/constraint_system/exclusion_matrix.h"
 #include "tachyon/zk/plonk/constraint_system/selector_assignment.h"
@@ -238,6 +240,8 @@ class SelectorCompressor {
       size_t n, const std::vector<SelectorDescription>& combination) {
     // Now, compute the selector and combination assignments.
     std::vector<F> combination_assignment(n);
+    base::ParallelizeFill(combination_assignment, F::Zero(),
+                          /*threshold=*/math::ParallelizeThreshold::kFieldInit);
     size_t combination_len = combination.size();
     size_t combination_index = combination_assignments_.size();
     std::unique_ptr<Expression<F>> query = callback_.Run();
diff --git a/tachyon/zk/plonk/examples/shuffle_circuit.h b/tachyon/zk/plonk/examples/shuffle_circuit.h
index 8b94fee042..43195d3aeb 100644
--- a/tachyon/zk/plonk/examples/shuffle_circuit.h
+++ b/tachyon/zk/plonk/examples/shuffle_circuit.h
@@ -205,8 +205,10 @@ class ShuffleCircuit : public Circuit<ShuffleCircuitConfig<F, W>> {
 
   std::unique_ptr<Circuit<ShuffleCircuitConfig<F, W>>> WithoutWitness()
       const override {
-    std::vector<std::vector<F>> dummy_original_table(W, std::vector<F>(H));
-    std::vector<std::vector<F>> dummy_shuffled_table(W, std::vector<F>(H));
+    std::vector<std::vector<F>> dummy_original_table(
+        W, std::vector<F>(H, F::Zero()));
+    std::vector<std::vector<F>> dummy_shuffled_table(
+        W, std::vector<F>(H, F::Zero()));
     ShuffleCircuit dummy_circuit(std::move(dummy_original_table),
                                  std::move(dummy_shuffled_table));
     return std::make_unique<ShuffleCircuit>(std::move(dummy_circuit));
diff --git a/tachyon/zk/plonk/expressions/compress_expression_unittest.cc b/tachyon/zk/plonk/expressions/compress_expression_unittest.cc
index e764a090c9..2db7d2ae23 100644
--- a/tachyon/zk/plonk/expressions/compress_expression_unittest.cc
+++ b/tachyon/zk/plonk/expressions/compress_expression_unittest.cc
@@ -26,7 +26,7 @@ TEST_F(CompressExpressionTest, CompressExpressions) {
 
   size_t n = prover_->pcs().N();
   ProvingEvaluator<Evals> evaluator = evaluator_;
-  std::vector<F> expected(n);
+  std::vector<F> expected(n, F::Zero());
   for (size_t i = 0; i < expressions.size(); ++i) {
     F value = evaluator.Evaluate(expressions[i].get());
     for (size_t j = 0; j < n; ++j) {
diff --git a/tachyon/zk/plonk/halo2/BUILD.bazel b/tachyon/zk/plonk/halo2/BUILD.bazel
index 4f2ccc6273..7bb9ce8290 100644
--- a/tachyon/zk/plonk/halo2/BUILD.bazel
+++ b/tachyon/zk/plonk/halo2/BUILD.bazel
@@ -8,8 +8,10 @@ tachyon_cc_library(
     deps = [
         ":synthesizer",
         "//tachyon/base:logging",
+        "//tachyon/base:parallelize",
         "//tachyon/base/buffer:copyable",
         "//tachyon/base/containers:container_util",
+        "//tachyon/math/base:parallelize_threshold",
         "//tachyon/zk/plonk/base:multi_phase_ref_table",
         "@com_google_absl//absl/types:span",
     ],
diff --git a/tachyon/zk/plonk/halo2/argument_data.h b/tachyon/zk/plonk/halo2/argument_data.h
index 7fda2c4928..f5adc603b5 100644
--- a/tachyon/zk/plonk/halo2/argument_data.h
+++ b/tachyon/zk/plonk/halo2/argument_data.h
@@ -11,6 +11,8 @@
 #include "tachyon/base/buffer/copyable.h"
 #include "tachyon/base/containers/container_util.h"
 #include "tachyon/base/logging.h"
+#include "tachyon/base/parallelize.h"
+#include "tachyon/math/base/parallelize_threshold.h"
 #include "tachyon/zk/plonk/base/multi_phase_ref_table.h"
 #include "tachyon/zk/plonk/halo2/argument_data.h"
 #include "tachyon/zk/plonk/halo2/synthesizer.h"
@@ -63,7 +65,8 @@ class ArgumentData {
     size_t n = prover->pcs().N();
     for (size_t i = 0; i < num_circuits; ++i) {
       for (Evals& instance_column : instance_columns_vec[i]) {
-        instance_column.evaluations().resize(n);
+        base::ParallelizeResize(instance_column.evaluations(), n, F::Zero(),
+                                math::ParallelizeThreshold::kFieldInit);
       }
     }
 
diff --git a/tachyon/zk/plonk/halo2/pinned_evaluation_domain.h b/tachyon/zk/plonk/halo2/pinned_evaluation_domain.h
index 02715aa12a..82fb6d736a 100644
--- a/tachyon/zk/plonk/halo2/pinned_evaluation_domain.h
+++ b/tachyon/zk/plonk/halo2/pinned_evaluation_domain.h
@@ -22,7 +22,6 @@ namespace zk::plonk::halo2 {
 template <typename F>
 class PinnedEvaluationDomain {
  public:
-  PinnedEvaluationDomain() = default;
   PinnedEvaluationDomain(uint32_t k, uint32_t extended_k, const F& omega)
       : k_(k), extended_k_(extended_k), omega_(omega) {}
   PinnedEvaluationDomain(uint32_t k, uint32_t extended_k, F&& omega)
@@ -38,8 +37,8 @@ class PinnedEvaluationDomain {
   const F& omega() const { return omega_; }
 
  private:
-  uint32_t k_ = 0;
-  uint32_t extended_k_ = 0;
+  uint32_t k_;
+  uint32_t extended_k_;
   F omega_;
 };
 
diff --git a/tachyon/zk/plonk/halo2/proof.h b/tachyon/zk/plonk/halo2/proof.h
index 0ed4b6db25..de08e82ef8 100644
--- a/tachyon/zk/plonk/halo2/proof.h
+++ b/tachyon/zk/plonk/halo2/proof.h
@@ -23,18 +23,18 @@ template <typename F, typename C>
 struct ProofBase {
   std::vector<std::vector<C>> advices_commitments_vec;
   std::vector<F> challenges;
-  F theta;
-  F beta;
-  F gamma;
+  F theta = F::Zero();
+  F beta = F::Zero();
+  F gamma = F::Zero();
   std::vector<std::vector<C>> permutation_product_commitments_vec;
   C vanishing_random_poly_commitment;
-  F y;
+  F y = F::Zero();
   std::vector<C> vanishing_h_poly_commitments;
-  F x;
+  F x = F::Zero();
   std::vector<std::vector<F>> instance_evals_vec;
   std::vector<std::vector<F>> advice_evals_vec;
   std::vector<F> fixed_evals;
-  F vanishing_random_eval;
+  F vanishing_random_eval = F::Zero();
   std::vector<F> common_permutation_evals;
   std::vector<std::vector<F>> permutation_product_evals_vec;
   std::vector<std::vector<F>> permutation_product_next_evals_vec;
@@ -44,13 +44,13 @@ struct ProofBase {
   std::vector<std::vector<F>> shuffle_product_next_evals_vec;
 
   // auxiliary values
-  F l_first;
-  F l_blind;
-  F l_last;
-  F x_next;
-  F x_prev;
-  F x_last;
-  F x_n;
+  F l_first = F::Zero();
+  F l_blind = F::Zero();
+  F l_last = F::Zero();
+  F x_next = F::Zero();
+  F x_prev = F::Zero();
+  F x_last = F::Zero();
+  F x_n = F::Zero();
 
   bool operator==(const ProofBase& other) const {
     return advices_commitments_vec == other.advices_commitments_vec &&
diff --git a/tachyon/zk/plonk/halo2/proof_reader.h b/tachyon/zk/plonk/halo2/proof_reader.h
index d4ea924043..32db20da25 100644
--- a/tachyon/zk/plonk/halo2/proof_reader.h
+++ b/tachyon/zk/plonk/halo2/proof_reader.h
@@ -56,7 +56,9 @@ class ProofReader {
     proof_.advices_commitments_vec.resize(num_circuits_);
     size_t num_advice_columns = constraint_system.num_advice_columns();
     for (size_t i = 0; i < num_circuits_; ++i) {
-      proof_.advices_commitments_vec[i].resize(num_advice_columns);
+      base::ParallelizeResize(proof_.advices_commitments_vec[i],
+                              num_advice_columns, C::Zero(),
+                              math::ParallelizeThreshold::kFieldInit);
     }
     proof_.challenges.reserve(constraint_system.challenge_phases().size());
     for (Phase current_phase : constraint_system.GetPhases()) {
diff --git a/tachyon/zk/plonk/halo2/synthesizer.h b/tachyon/zk/plonk/halo2/synthesizer.h
index 924cffc7b6..3d36b9f20c 100644
--- a/tachyon/zk/plonk/halo2/synthesizer.h
+++ b/tachyon/zk/plonk/halo2/synthesizer.h
@@ -25,15 +25,13 @@ class Synthesizer {
   Synthesizer() = default;
   Synthesizer(size_t num_circuits, const ConstraintSystem<F>* constraint_system)
       : num_circuits_(num_circuits), constraint_system_(constraint_system) {
-    advice_columns_vec_.resize(num_circuits_);
-    advice_blinds_vec_.resize(num_circuits_);
-    for (size_t i = 0; i < num_circuits_; ++i) {
-      // And these may be assigned with random order.
-      advice_columns_vec_[i] =
-          std::vector<Evals>(constraint_system->num_advice_columns());
-      advice_blinds_vec_[i] =
-          std::vector<F>(constraint_system->num_advice_columns());
-    }
+    // NOTE(batzor): These vectors are initialized below in |SetAdviceColumn| so
+    // it is safe to keep it uninitialized here.
+    advice_columns_vec_.resize(
+        num_circuits_,
+        std::vector<Evals>(constraint_system->num_advice_columns()));
+    advice_blinds_vec_.resize(
+        num_circuits_, std::vector<F>(constraint_system->num_advice_columns()));
   }
 
   // Synthesize circuit and store advice columns.
diff --git a/tachyon/zk/plonk/halo2/verifier.h b/tachyon/zk/plonk/halo2/verifier.h
index 15e9a7727e..11715adb79 100644
--- a/tachyon/zk/plonk/halo2/verifier.h
+++ b/tachyon/zk/plonk/halo2/verifier.h
@@ -182,8 +182,19 @@ class Verifier : public VerifierBase<typename _PS::PCS> {
 
   std::vector<Commitment> CommitColumns(const std::vector<Evals>& columns) {
     return base::Map(columns, [this](const Evals& column) {
-      std::vector<F> expanded_evals = column.evaluations();
-      expanded_evals.resize(this->pcs_.N());
+      std::vector<F> expanded_evals(this->pcs_.N());
+      const std::vector<F>& evals = column.evaluations();
+      base::Parallelize(
+          expanded_evals,
+          [&evals](absl::Span<F> chunk, size_t chunk_offset,
+                   size_t chunk_size) {
+            size_t i = chunk_offset * chunk_size;
+            for (F& eval : chunk) {
+              eval = i < evals.size() ? evals[i] : F::Zero();
+              ++i;
+            }
+          },
+          /*threshold=*/math::ParallelizeThreshold::kFieldInit);
       return this->Commit(Evals(expanded_evals));
     });
   }
diff --git a/tachyon/zk/plonk/keys/key.h b/tachyon/zk/plonk/keys/key.h
index 8c333422c1..046d631624 100644
--- a/tachyon/zk/plonk/keys/key.h
+++ b/tachyon/zk/plonk/keys/key.h
@@ -93,6 +93,8 @@ class TACHYON_EXPORT Key {
 
     result->fixed_columns =
         base::Map(assembly.fixed_columns(), [](const RationalEvals& evals) {
+          // NOTE(batzor): This vector is initialized below in |BatchEvaluate|
+          // so it is safe to keep it uninitialized here.
           std::vector<F> result(evals.evaluations().size());
           CHECK(math::RationalField<F>::BatchEvaluate(evals.evaluations(),
                                                       &result));
diff --git a/tachyon/zk/plonk/permutation/grand_product_argument.h b/tachyon/zk/plonk/permutation/grand_product_argument.h
index fd83079d9f..3249f9d8ef 100644
--- a/tachyon/zk/plonk/permutation/grand_product_argument.h
+++ b/tachyon/zk/plonk/permutation/grand_product_argument.h
@@ -22,6 +22,10 @@ class GrandProductArgument {
 
     // NOTE(chokobole): It's safe to downcast because domain is already checked.
     RowIndex size = static_cast<RowIndex>(prover->pcs().N());
+
+    // NOTE(batzor): This vector is initialized below in |denominator_callback|
+    // so it is safe to keep it uninitialized here. First element is set to
+    // |last_z| in |DoCreatePoly|.
     std::vector<F> z(size + 1);
     absl::Span<F> grand_product = absl::MakeSpan(z).subspan(1);
 
@@ -44,6 +48,10 @@ class GrandProductArgument {
 
     // NOTE(chokobole): It's safe to downcast because domain is already checked.
     RowIndex size = static_cast<RowIndex>(prover->pcs().N());
+
+    // NOTE(batzor): This vector is initialized below in |denominator_callback|
+    // so it is safe to keep it uninitialized here. First element is set to
+    // |last_z| in |DoCreatePoly|.
     std::vector<F> z(size + 1);
     absl::Span<F> grand_product = absl::MakeSpan(z).subspan(1);
 
diff --git a/tachyon/zk/plonk/vanishing/BUILD.bazel b/tachyon/zk/plonk/vanishing/BUILD.bazel
index 55af77f93b..a6c038e636 100644
--- a/tachyon/zk/plonk/vanishing/BUILD.bazel
+++ b/tachyon/zk/plonk/vanishing/BUILD.bazel
@@ -109,8 +109,10 @@ tachyon_cc_library(
     deps = [
         ":vanishing_argument",
         ":vanishing_utils",
+        "//tachyon/base:parallelize",
         "//tachyon/base:profiler",
         "//tachyon/crypto/commitments:polynomial_openings",
+        "//tachyon/math/base:parallelize_threshold",
         "//tachyon/zk/lookup:prover",
         "//tachyon/zk/plonk/base:multi_phase_ref_table",
         "//tachyon/zk/plonk/keys:proving_key",
diff --git a/tachyon/zk/plonk/vanishing/circuit_polynomial_builder.h b/tachyon/zk/plonk/vanishing/circuit_polynomial_builder.h
index 78f96ab6bf..d26417bb3a 100644
--- a/tachyon/zk/plonk/vanishing/circuit_polynomial_builder.h
+++ b/tachyon/zk/plonk/vanishing/circuit_polynomial_builder.h
@@ -213,6 +213,8 @@ class CircuitPolynomialBuilder {
       UpdateLPolyCosets();
 
       std::vector<F> value_part(static_cast<size_t>(n_));
+      base::ParallelizeFill(value_part, F::Zero(),
+                            math::ParallelizeThreshold::kFieldInit);
       size_t circuit_num = poly_tables_.size();
       for (size_t j = 0; j < circuit_num; ++j) {
         VLOG(1) << "BuildExtendedCircuitColumn part: " << i << " circuit: ("
diff --git a/tachyon/zk/plonk/vanishing/vanishing_prover_impl.h b/tachyon/zk/plonk/vanishing/vanishing_prover_impl.h
index 9cbf58c5dc..dd88acf119 100644
--- a/tachyon/zk/plonk/vanishing/vanishing_prover_impl.h
+++ b/tachyon/zk/plonk/vanishing/vanishing_prover_impl.h
@@ -10,7 +10,9 @@
 #include <utility>
 #include <vector>
 
+#include "tachyon/base/parallelize.h"
 #include "tachyon/base/profiler.h"
+#include "tachyon/math/base/parallelize_threshold.h"
 #include "tachyon/zk/plonk/vanishing/vanishing_argument.h"
 #include "tachyon/zk/plonk/vanishing/vanishing_prover.h"
 
@@ -179,6 +181,8 @@ void VanishingProver<Poly, Evals, ExtendedPoly, ExtendedEvals>::BatchEvaluate(
       base::Map(h_chunks.begin(), h_chunks.end(),
                 [](absl::Span<F> h_piece) { return h_piece; });
   std::vector<F> coeffs(n);
+  base::ParallelizeFill(coeffs, F::Zero(),
+                        /*threshold=*/math::ParallelizeThreshold::kFieldInit);
   for (size_t i = h_pieces.size() - 1; i != SIZE_MAX; --i) {
     OMP_PARALLEL_FOR(size_t j = 0; j < n; ++j) {
       coeffs[j] *= x_n;
diff --git a/tachyon/zk/r1cs/constraint_system/quadratic_arithmetic_program.h b/tachyon/zk/r1cs/constraint_system/quadratic_arithmetic_program.h
index 95bb323448..a35758f1de 100644
--- a/tachyon/zk/r1cs/constraint_system/quadratic_arithmetic_program.h
+++ b/tachyon/zk/r1cs/constraint_system/quadratic_arithmetic_program.h
@@ -28,7 +28,7 @@ template <typename F>
 F EvaluateConstraint(const std::vector<Cell<F>>& cells,
                      absl::Span<const F> assignments) {
   TRACE_EVENT("ProofGeneration", "QAP::EvaluateConstraint");
-  F sum;
+  F sum = F::Zero();
   for (const Cell<F>& cell : cells) {
     if (cell.coefficient.IsOne()) {
       sum += assignments[cell.index];
@@ -71,6 +71,11 @@ class QuadraticArithmeticProgram {
     std::vector<F> a(num_qap_variables + 1);
     std::vector<F> b(num_qap_variables + 1);
     std::vector<F> c(num_qap_variables + 1);
+    OMP_PARALLEL_FOR(size_t i = 0; i <= num_qap_variables; ++i) {
+      a[i] = F::Zero();
+      b[i] = F::Zero();
+      c[i] = F::Zero();
+    }
 
     // clang-format off
     // |a[i]| = lₙ₊ᵢ(x) +  Σⱼ₌₀..ₙ₋₁ (lⱼ(x) * Aⱼ,ᵢ) (if i < |num_instance_variables|)
@@ -133,6 +138,12 @@ class QuadraticArithmeticProgram {
     std::vector<F> a(domain->size());
     std::vector<F> b(domain->size());
     std::vector<F> c(domain->size());
+    OMP_PARALLEL_FOR(size_t i = matrices.num_constraints; i < domain->size();
+                     ++i) {
+      a[i] = F::Zero();
+      b[i] = F::Zero();
+      c[i] = F::Zero();
+    }
 
     // clang-format off
     // |a[i]| = Σⱼ₌₀..ₘ (xⱼ * Aᵢ,ⱼ)    (if i < |num_constraints|)
diff --git a/vendors/circom/circomlib/circuit/quadratic_arithmetic_program.h b/vendors/circom/circomlib/circuit/quadratic_arithmetic_program.h
index c4540a9f14..1fb6e53acb 100644
--- a/vendors/circom/circomlib/circuit/quadratic_arithmetic_program.h
+++ b/vendors/circom/circomlib/circuit/quadratic_arithmetic_program.h
@@ -32,7 +32,12 @@ class QuadraticArithmeticProgram {
     TRACE_EVENT_BEGIN("Subtask", "ABC Initialization");
     std::vector<F> a(domain->size());
     std::vector<F> b(domain->size());
+    // NOTE: |c| is initialized later with |c[i]| = |a[i]| * |b[i]|.
     std::vector<F> c(domain->size());
+    OMP_PARALLEL_FOR(size_t i = 0; i < domain->size(); ++i) {
+      a[i] = F::Zero();
+      b[i] = F::Zero();
+    }
     TRACE_EVENT_END("Subtask");
 
     TRACE_EVENT_BEGIN("Subtask", "AB Calculation");