#include <algorithm>
#include <future>
#include <iostream>
#include <random>
#include <string>

// A function to compute the scalar product of vectors a and b using the
// specified number of threads.
double dot_product(std::vector<double> const &a, std::vector<double> const &b,
                   int dot_threads = 4);

int main(int argc, char *argv[]) {
  if (argc != 3) {
    std::cerr << "Please, give the number of elements and number of threads in "
                 "the command line.\n";
    return 1;
  }

  int const N = std::stoi(argv[1]);
  if (N <= 0) {
    std::cerr << "Number of elements must be positive.\n";
    return 2;
  }

  int const num_threads = std::stoi(argv[2]);
  if (num_threads <= 0) {
    std::cerr << "Number of threads must be positive.\n";
    return 2;
  }

  // Random number generation
  std::random_device rd;
  std::default_random_engine rand_gen(rd());
  std::uniform_real_distribution<double> dist(0, 1);
  auto gen_uniform_rand = [&]() { return dist(rand_gen); };

  // Initialize the vectors with random values from 0 to 1.
  std::vector<double> a(N), b(N);
  std::generate(begin(a), end(a), gen_uniform_rand);
  std::generate(begin(b), end(b), gen_uniform_rand);

  // Compute the dot product in parallel.
  double dotprod = dot_product(a, b, num_threads);

  std::cout << "The dot product is " << dotprod << std::endl;

  return 0;
}

// A function to compute the scalar product of vectors a and b using the
// specified number of threads.
double dot_product(std::vector<double> const &a, std::vector<double> const &b,
                   int dot_threads) {
  std::vector<std::future<double>> results;

  // We assign a block of consecutive elements to each thread.
  auto task_size = a.size() / dot_threads; // Size of the block.
  auto part_dot = [&a = std::as_const(a), &b = std::as_const(b), dot_threads,
                   task_size](int tid) {
    double s = 0.0;
    size_t istart = tid * task_size;
    // Last thread get all remaining elements.
    size_t ifinish = tid != dot_threads - 1 ? (tid + 1) * task_size : a.size();
    for (size_t j = istart; j < ifinish; ++j) {
      s += a[j] * b[j];
    }
    return s;
  };

  // Start the threads.
  for (int i = 0; i < dot_threads; ++i) {
    results.emplace_back(std::async(std::launch::async, part_dot, i));
  }

  // Wait for the results and sum.
  double s = 0.0;
  for (auto &res : results) {
    s += res.get();
  }

  return s;
}