#pragma once
#include "mayer_primes.h"
#include <cassert>
#include <vector>


/**
  This function returns the number of possible decompositions of an integer into a sum of two prime numbers.
  @param[in]    k            first integer
  @param[out]   count        number of decompositions
*/
size_t single_goldbach(size_t k);


/**
  This function returns the number of possible decompositions into a sum of two prime numbers of all even integers in the interval [4,n].
  @param[in]    n                   upper integer bound
  @param[out]   count_vector        vector containing the number of decompositions for a natural number the corresponding index
*/
std::vector<size_t> count_goldbach(size_t n);


/** 	 Vector @p b adds its elements to vector @p a .
	@param[in] a	vector
	@param[in] b	vector
	@return 	     a+=b componentwise
*/
template<class T>
std::vector<T> &operator+=(std::vector<T> &a, std::vector<T> const &b)
{
    assert(a.size()==b.size());
    for (size_t k = 0; k < a.size(); ++k) {
        a[k] += b[k];
    }
    return a;
}

// Declare the reduction operation in OpenMP for an STL-vector
//   omp_out += omp_in  requires operator+=(vector<int> &, vector<int> const &) from above
// ------------------------------------------------------------
// https://scc.ustc.edu.cn/zlsc/tc4600/intel/2016.0.109/compiler_c/common/core/GUID-7312910C-D175-4544-99C5-29C12D980744.htm
// https://gist.github.com/eruffaldi/7180bdec4c8c9a11f019dd0ba9a2d68c
// https://stackoverflow.com/questions/29633531/user-defined-reduction-on-vector-of-varying-size
//  see also p.74ff in  https://www.fz-juelich.de/ias/jsc/EN/AboutUs/Staff/Hagemeier_A/docs-parallel-programming/OpenMP-Slides.pdf
#pragma omp declare reduction(VecAdd : std::vector<size_t>  : omp_out += omp_in) initializer (omp_priv=omp_orig)