feat: parallelization update for polynomials

barretenberg/cpp/src/barretenberg/common/thread_utils.cpp

@@ -0,0 +1,20 @@
+#include "thread_utils.hpp"
+
+/**
+ * @brief calculates number of threads to create based on minimum iterations per thread
+ * @details Finds the number of cpus with get_num_cpus(), and calculates `desired_num_threads`
+ * Returns the min of `desired_num_threads` and `max_num_theads`.
+ * Note that it will not calculate a power of 2 necessarily, use `calc_num_threads_pow2` instead
+ *
+ * @param num_iterations
+ * @param min_iterations_per_thread
+ * @return size_t
+ */
+size_t calc_num_threads(size_t num_iterations, size_t min_iterations_per_thread)
+{
+    size_t max_num_threads = get_num_cpus(); // number of available threads
+    size_t desired_num_threads = num_iterations / min_iterations_per_thread;
+    size_t num_threads = std::min(desired_num_threads, max_num_threads); // fewer than max if justified
+    num_threads = num_threads > 0 ? num_threads : 1;
+    return num_threads;
+}

barretenberg/cpp/src/barretenberg/common/thread_utils.hpp

@@ -0,0 +1,15 @@
+#include "thread.hpp"
+
+const size_t MIN_ITERS_PER_THREAD = 1 << 4;
+
+/**
+ * @brief calculates number of threads to create based on minimum iterations per thread
+ * @details Finds the number of cpus with get_num_cpus(), and calculates `desired_num_threads`
+ * Returns the min of `desired_num_threads` and `max_num_theads`.
+ * Note that it will not calculate a power of 2 necessarily, use `calc_num_threads_pow2` instead
+ *
+ * @param num_iterations
+ * @param min_iterations_per_thread
+ * @return size_t
+ */
+size_t calc_num_threads(size_t num_iterations, size_t min_iterations_per_thread);

barretenberg/cpp/src/barretenberg/polynomials/polynomial.cpp

@@ -1,7 +1,3 @@
-#include "polynomial.hpp"
-#include "barretenberg/common/assert.hpp"
-#include "barretenberg/common/slab_allocator.hpp"
-#include "polynomial_arithmetic.hpp"
 #include <cstddef>
 #include <fcntl.h>
 #include <list>
@@ -11,6 +7,13 @@
 #include <unordered_map>
 #include <utility>
 
+#include "barretenberg/common/assert.hpp"
+#include "barretenberg/common/slab_allocator.hpp"
+#include "barretenberg/common/thread.hpp"
+#include "barretenberg/common/thread_utils.hpp"
+#include "polynomial.hpp"
+#include "polynomial_arithmetic.hpp"
+
 namespace barretenberg {
 /**
  * Constructors / Destructors
@@ -272,25 +275,34 @@ template <typename Fr> void Polynomial<Fr>::add_scaled(std::span<const Fr> other
     const size_t other_size = other.size();
     ASSERT(in_place_operation_viable(other_size));
 
-    /** TODO parallelize using some kind of generic evaluation domain
-     *  we really only need to know the thread size, but we don't need all the FFT roots
+    /** Calculates number of threads with calc_num_threads
+     *  Possible improvements: standardize this parallelization code
      */
-    for (size_t i = 0; i < other_size; ++i) {
-        coefficients_.get()[i] += scaling_factor * other[i];
-    }
+    size_t num_threads = calc_num_threads(other_size, MIN_ITERS_PER_THREAD);
+    size_t range_per_thread = other_size / num_threads;
+    size_t leftovers = other_size - (range_per_thread * num_threads);
+    parallel_for(num_threads, [&](size_t j) {
+        size_t offset = j * range_per_thread;
+        size_t end = (j == num_threads - 1) ? offset + range_per_thread + leftovers : offset + range_per_thread;
+        for (size_t i = offset; i < end; ++i)
+            coefficients_.get()[i] += scaling_factor * other[i];
+    });
 }
 
 template <typename Fr> Polynomial<Fr>& Polynomial<Fr>::operator+=(std::span<const Fr> other)
 {
     const size_t other_size = other.size();
     ASSERT(in_place_operation_viable(other_size));
 
-    /** TODO parallelize using some kind of generic evaluation domain
-     *  we really only need to know the thread size, but we don't need all the FFT roots
-     */
-    for (size_t i = 0; i < other_size; ++i) {
-        coefficients_.get()[i] += other[i];
-    }
+    size_t num_threads = calc_num_threads(other_size, MIN_ITERS_PER_THREAD);
+    size_t range_per_thread = other_size / num_threads;
+    size_t leftovers = other_size - (range_per_thread * num_threads);
+    parallel_for(num_threads, [&](size_t j) {
+        size_t offset = j * range_per_thread;
+        size_t end = (j == num_threads - 1) ? offset + range_per_thread + leftovers : offset + range_per_thread;
+        for (size_t i = offset; i < end; ++i)
+            coefficients_.get()[i] += other[i];
+    });
 
     return *this;
 }
@@ -300,23 +312,33 @@ template <typename Fr> Polynomial<Fr>& Polynomial<Fr>::operator-=(std::span<cons
     const size_t other_size = other.size();
     ASSERT(in_place_operation_viable(other_size));
 
-    /** TODO parallelize using some kind of generic evaluation domain
-     *  we really only need to know the thread size, but we don't need all the FFT roots
-     */
-    for (size_t i = 0; i < other_size; ++i) {
-        coefficients_.get()[i] -= other[i];
-    }
+    size_t num_threads = calc_num_threads(other_size, MIN_ITERS_PER_THREAD);
+    size_t range_per_thread = other_size / num_threads;
+    size_t leftovers = other_size - (range_per_thread * num_threads);
+    parallel_for(num_threads, [&](size_t j) {
+        size_t offset = j * range_per_thread;
+        size_t end = (j == num_threads - 1) ? offset + range_per_thread + leftovers : offset + range_per_thread;
+        for (size_t i = offset; i < end; ++i)
+            coefficients_.get()[i] -= other[i];
+    });
 
     return *this;
 }
 
-template <typename Fr> Polynomial<Fr>& Polynomial<Fr>::operator*=(const Fr scaling_facor)
+template <typename Fr> Polynomial<Fr>& Polynomial<Fr>::operator*=(const Fr scaling_factor)
 {
     ASSERT(in_place_operation_viable());
 
-    for (size_t i = 0; i < size_; ++i) {
-        coefficients_.get()[i] *= scaling_facor;
-    }
+    size_t num_threads = calc_num_threads(size_, MIN_ITERS_PER_THREAD);
+    size_t range_per_thread = size_ / num_threads;
+    size_t leftovers = size_ - (range_per_thread * num_threads);
+    parallel_for(num_threads, [&](size_t j) {
+        size_t offset = j * range_per_thread;
+        size_t end = (j == num_threads - 1) ? offset + range_per_thread + leftovers : offset + range_per_thread;
+        for (size_t i = offset; i < end; ++i)
+            coefficients_.get()[i] *= scaling_factor;
+    });
+
     return *this;
 }
 

barretenberg/cpp/src/barretenberg/polynomials/polynomial.hpp

@@ -176,7 +176,7 @@ template <typename Fr> class Polynomial {
      *
      * @param scaling_factor s
      */
-    Polynomial& operator*=(const Fr scaling_facor);
+    Polynomial& operator*=(const Fr scaling_factor);
 
     /**
      * @brief evaluates p(X) = ∑ᵢ aᵢ⋅Xⁱ considered as multi-linear extension p(X₀,…,Xₘ₋₁) = ∑ᵢ aᵢ⋅Lᵢ(X₀,…,Xₘ₋₁)