 
 }
 
-Point3D reduce(PointCloud<Point3D> &cloud, func& f) {
+Point3D reduce(PointCloud<Point3D> &cloud, func& f, double identity_element) {
 
     Point3D *d_points;
     Point3D *d_result;
-    Point3D *d_reduce;
 
     auto num_points = cloud.points.size();
+
+    auto size = next_power_of_two(num_points);
 
     func* h_f;
     func* d_f;


    cudaFree(d_points);
    cudaFree(d_result);
    cudaFree(d_reduce);

    cudaFree(d_f);
    free(h_f);

#include <float.h>
#include "kernels.h"
#include "../utils/utils.h"

typedef double (*func)(double n, double m);

__global__ void reduce(const Point3D * const d_points, const uint number,
                       Point3D *d_result, func* f, double id_element) {

    const int idx = threadIdx.x;
    const int size = blockDim.x;

    extern __shared__ Point3D d_reduce[];

    if (idx < number) {
        d_reduce[idx] = d_points[idx];
    } else {
        d_reduce[idx] = Point3D(id_element,id_element,id_element);
    }

    __syncthreads();

    for (unsigned int s = size / 2; s > 0; s >>= 1)
    {
        if (idx < s)
        {
