int n = 0;
double val = 0;
for (int i = 0; i < p_count; ++i) {
std::vector<double> cpoint;
for (unsigned int j = 0; j < dim; ++j) {
cpoint.push_back(p[dim * i + j]);
}
if (parea(cpoint)) {
n++;
val += pfunction(cpoint);
}
}
val = val / n;
double area = std::pow(sd, dim) * n / p_count;