diff --git a/pykdtree/_kdtree_core.c b/pykdtree/_kdtree_core.c
index 2d6862b..4d130bd 100644
--- a/pykdtree/_kdtree_core.c
+++ b/pykdtree/_kdtree_core.c
@@ -29,7 +29,7 @@ Anne M. Archibald and libANN by David M. Mount and Sunil Arya.
 #include <float.h>
 
 #define PA(i,d)			(pa[no_dims * pidx[i] + d])
-#define PASWAP(a,b) { uint32_t tmp = pidx[a]; pidx[a] = pidx[b]; pidx[b] = tmp; }
+#define PASWAP(a,b) { uint64_t tmp = pidx[a]; pidx[a] = pidx[b]; pidx[b] = tmp; }
 
 #ifdef _MSC_VER
 #define restrict __restrict
@@ -40,8 +40,8 @@ typedef struct
 {
     float cut_val;
     int8_t cut_dim;
-    uint32_t start_idx;
-    uint32_t n;
+    uint64_t start_idx;
+    uint64_t n;
     float cut_bounds_lv;
     float cut_bounds_hv;
     struct Node_float *left_child;
@@ -52,7 +52,7 @@ typedef struct
 {
     float *bbox;
     int8_t no_dims;
-    uint32_t *pidx;
+    uint64_t *pidx;
     struct Node_float *root;
 } Tree_float;
 
@@ -61,8 +61,8 @@ typedef struct
 {
     double cut_val;
     int8_t cut_dim;
-    uint32_t start_idx;
-    uint32_t n;
+    uint64_t start_idx;
+    uint64_t n;
     double cut_bounds_lv;
     double cut_bounds_hv;
     struct Node_double *left_child;
@@ -73,58 +73,58 @@ typedef struct
 {
     double *bbox;
     int8_t no_dims;
-    uint32_t *pidx;
+    uint64_t *pidx;
     struct Node_double *root;
 } Tree_double;
 
 
 
-void insert_point_float(uint32_t *closest_idx, float *closest_dist, uint32_t pidx, float cur_dist, uint32_t k);
-void get_bounding_box_float(float *pa, uint32_t *pidx, int8_t no_dims, uint32_t n, float *bbox);
-int partition_float(float *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, float *bbox, int8_t *cut_dim,
-              float *cut_val, uint32_t *n_lo);
-Tree_float* construct_tree_float(float *pa, int8_t no_dims, uint32_t n, uint32_t bsp);
-Node_float* construct_subtree_float(float *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, uint32_t bsp, float *bbox);
-Node_float * create_node_float(uint32_t start_idx, uint32_t n, int is_leaf);
+void insert_point_float(uint64_t *closest_idx, float *closest_dist, uint64_t pidx, float cur_dist, uint64_t k);
+void get_bounding_box_float(float *pa, uint64_t *pidx, int8_t no_dims, uint64_t n, float *bbox);
+int partition_float(float *pa, uint64_t *pidx, int8_t no_dims, uint64_t start_idx, uint64_t n, float *bbox, int8_t *cut_dim,
+              float *cut_val, uint64_t *n_lo);
+Tree_float* construct_tree_float(float *pa, int8_t no_dims, uint64_t n, uint64_t bsp);
+Node_float* construct_subtree_float(float *pa, uint64_t *pidx, int8_t no_dims, uint64_t start_idx, uint64_t n, uint64_t bsp, float *bbox);
+Node_float * create_node_float(uint64_t start_idx, uint64_t n, int is_leaf);
 void delete_subtree_float(Node_float *root);
 void delete_tree_float(Tree_float *tree);
 void print_tree_float(Node_float *root, int level);
 float calc_dist_float(float *point1_coord, float *point2_coord, int8_t no_dims);
 float get_cube_offset_float(int8_t dim, float *point_coord, float *bbox);
 float get_min_dist_float(float *point_coord, int8_t no_dims, float *bbox);
-void search_leaf_float(float *restrict pa, uint32_t *restrict pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, float *restrict point_coord,
-                 uint32_t k, uint32_t *restrict closest_idx, float *restrict closest_dist);
-void search_leaf_float_mask(float *restrict pa, uint32_t *restrict pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, float *restrict point_coord,
-                 uint32_t k, uint8_t *restrict mask, uint32_t *restrict closest_idx, float *restrict closest_dist);
-void search_splitnode_float(Node_float *root, float *pa, uint32_t *pidx, int8_t no_dims, float *point_coord,
-                      float min_dist, uint32_t k, float distance_upper_bound, float eps_fac, uint8_t *mask, uint32_t *  closest_idx, float *closest_dist);
+void search_leaf_float(float *restrict pa, uint64_t *restrict pidx, int8_t no_dims, uint64_t start_idx, uint64_t n, float *restrict point_coord,
+                 uint64_t k, uint64_t *restrict closest_idx, float *restrict closest_dist);
+void search_leaf_float_mask(float *restrict pa, uint64_t *restrict pidx, int8_t no_dims, uint64_t start_idx, uint64_t n, float *restrict point_coord,
+                 uint64_t k, uint8_t *restrict mask, uint64_t *restrict closest_idx, float *restrict closest_dist);
+void search_splitnode_float(Node_float *root, float *pa, uint64_t *pidx, int8_t no_dims, float *point_coord,
+                      float min_dist, uint64_t k, float distance_upper_bound, float eps_fac, uint8_t *mask, uint64_t *  closest_idx, float *closest_dist);
 void search_tree_float(Tree_float *tree, float *pa, float *point_coords,
-                 uint32_t num_points, uint32_t k,  float distance_upper_bound,
-                 float eps, uint8_t *mask, uint32_t *closest_idxs, float *closest_dists);
+                 uint64_t num_points, uint64_t k,  float distance_upper_bound,
+                 float eps, uint8_t *mask, uint64_t *closest_idxs, float *closest_dists);
 
 
-void insert_point_double(uint32_t *closest_idx, double *closest_dist, uint32_t pidx, double cur_dist, uint32_t k);
-void get_bounding_box_double(double *pa, uint32_t *pidx, int8_t no_dims, uint32_t n, double *bbox);
-int partition_double(double *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, double *bbox, int8_t *cut_dim,
-              double *cut_val, uint32_t *n_lo);
-Tree_double* construct_tree_double(double *pa, int8_t no_dims, uint32_t n, uint32_t bsp);
-Node_double* construct_subtree_double(double *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, uint32_t bsp, double *bbox);
-Node_double * create_node_double(uint32_t start_idx, uint32_t n, int is_leaf);
+void insert_point_double(uint64_t *closest_idx, double *closest_dist, uint64_t pidx, double cur_dist, uint64_t k);
+void get_bounding_box_double(double *pa, uint64_t *pidx, int8_t no_dims, uint64_t n, double *bbox);
+int partition_double(double *pa, uint64_t *pidx, int8_t no_dims, uint64_t start_idx, uint64_t n, double *bbox, int8_t *cut_dim,
+              double *cut_val, uint64_t *n_lo);
+Tree_double* construct_tree_double(double *pa, int8_t no_dims, uint64_t n, uint64_t bsp);
+Node_double* construct_subtree_double(double *pa, uint64_t *pidx, int8_t no_dims, uint64_t start_idx, uint64_t n, uint64_t bsp, double *bbox);
+Node_double * create_node_double(uint64_t start_idx, uint64_t n, int is_leaf);
 void delete_subtree_double(Node_double *root);
 void delete_tree_double(Tree_double *tree);
 void print_tree_double(Node_double *root, int level);
 double calc_dist_double(double *point1_coord, double *point2_coord, int8_t no_dims);
 double get_cube_offset_double(int8_t dim, double *point_coord, double *bbox);
 double get_min_dist_double(double *point_coord, int8_t no_dims, double *bbox);
-void search_leaf_double(double *restrict pa, uint32_t *restrict pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, double *restrict point_coord,
-                 uint32_t k, uint32_t *restrict closest_idx, double *restrict closest_dist);
-void search_leaf_double_mask(double *restrict pa, uint32_t *restrict pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, double *restrict point_coord,
-                 uint32_t k, uint8_t *restrict mask, uint32_t *restrict closest_idx, double *restrict closest_dist);
-void search_splitnode_double(Node_double *root, double *pa, uint32_t *pidx, int8_t no_dims, double *point_coord,
-                      double min_dist, uint32_t k, double distance_upper_bound, double eps_fac, uint8_t *mask, uint32_t *  closest_idx, double *closest_dist);
+void search_leaf_double(double *restrict pa, uint64_t *restrict pidx, int8_t no_dims, uint64_t start_idx, uint64_t n, double *restrict point_coord,
+                 uint64_t k, uint64_t *restrict closest_idx, double *restrict closest_dist);
+void search_leaf_double_mask(double *restrict pa, uint64_t *restrict pidx, int8_t no_dims, uint64_t start_idx, uint64_t n, double *restrict point_coord,
+                 uint64_t k, uint8_t *restrict mask, uint64_t *restrict closest_idx, double *restrict closest_dist);
+void search_splitnode_double(Node_double *root, double *pa, uint64_t *pidx, int8_t no_dims, double *point_coord,
+                      double min_dist, uint64_t k, double distance_upper_bound, double eps_fac, uint8_t *mask, uint64_t *  closest_idx, double *closest_dist);
 void search_tree_double(Tree_double *tree, double *pa, double *point_coords,
-                 uint32_t num_points, uint32_t k,  double distance_upper_bound,
-                 double eps, uint8_t *mask, uint32_t *closest_idxs, double *closest_dists);
+                 uint64_t num_points, uint64_t k,  double distance_upper_bound,
+                 double eps, uint8_t *mask, uint64_t *closest_idxs, double *closest_dists);
 
 
 
@@ -137,7 +137,7 @@ Insert point into priority queue
     cur_dist : distance to point inserted
     k : number of neighbours
 ************************************************/
-void insert_point_float(uint32_t *closest_idx, float *closest_dist, uint32_t pidx, float cur_dist, uint32_t k)
+void insert_point_float(uint64_t *closest_idx, float *closest_dist, uint64_t pidx, float cur_dist, uint64_t k)
 {
     int i;
     for (i = k - 1; i > 0; i--)
@@ -165,11 +165,11 @@ Get the bounding box of a set of points
     n : number of points
     bbox : bounding box (return)
 ************************************************/
-void get_bounding_box_float(float *pa, uint32_t *pidx, int8_t no_dims, uint32_t n, float *bbox)
+void get_bounding_box_float(float *pa, uint64_t *pidx, int8_t no_dims, uint64_t n, float *bbox)
 {
     float cur;
     int8_t i, j;
-    uint32_t bbox_idx, i2;
+    uint64_t bbox_idx, i2;
 
     /* Use first data point to initialize */
     for (i = 0; i < no_dims; i++)
@@ -210,12 +210,12 @@ The sliding midpoint rule is used for the partitioning.
     cut_val : value of cutting point (return)
     n_lo : number of point below cutting plane (return)
 ************************************************/
-int partition_float(float *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, float *bbox, int8_t *cut_dim, float *cut_val, uint32_t *n_lo)
+int partition_float(float *pa, uint64_t *pidx, int8_t no_dims, uint64_t start_idx, uint64_t n, float *bbox, int8_t *cut_dim, float *cut_val, uint64_t *n_lo)
 {
     int8_t dim = 0, i;
-    uint32_t p, q, i2;
+    uint64_t p, q, i2;
     float size = 0, min_val, max_val, split, side_len, cur_val;
-    uint32_t end_idx = start_idx + n - 1;
+    uint64_t end_idx = start_idx + n - 1;
 
     /* Find largest bounding box side */
     for (i = 0; i < no_dims; i++)
@@ -275,7 +275,7 @@ int partition_float(float *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_id
            Minimum 1 point will be in lower box.
         */
 
-        uint32_t j = start_idx;
+        uint64_t j = start_idx;
         split = PA(j, dim);
         for (i2 = start_idx + 1; i2 <= end_idx; i2++)
         {
@@ -297,7 +297,7 @@ int partition_float(float *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_id
            Minimum 1 point will be in higher box.
         */
 
-        uint32_t j = end_idx;
+        uint64_t j = end_idx;
         split = PA(j, dim);
         for (i2 = start_idx; i2 < end_idx; i2++)
         {
@@ -331,14 +331,14 @@ Construct a sub tree over a range of data points.
     bsp : number of points per leaf
     bbox : bounding box of set of data points
 ************************************************/
-Node_float* construct_subtree_float(float *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, uint32_t bsp, float *bbox)
+Node_float* construct_subtree_float(float *pa, uint64_t *pidx, int8_t no_dims, uint64_t start_idx, uint64_t n, uint64_t bsp, float *bbox)
 {
     /* Create new node */
     int is_leaf = (n <= bsp);
     Node_float *root = create_node_float(start_idx, n, is_leaf);
     int rval;
     int8_t cut_dim;
-    uint32_t n_lo;
+    uint64_t n_lo;
     float cut_val, lv, hv;
     if (is_leaf)
     {
@@ -387,17 +387,17 @@ Construct a tree over data points.
     n :  number of data points
     bsp : number of points per leaf
 ************************************************/
-Tree_float* construct_tree_float(float *pa, int8_t no_dims, uint32_t n, uint32_t bsp)
+Tree_float* construct_tree_float(float *pa, int8_t no_dims, uint64_t n, uint64_t bsp)
 {
     Tree_float *tree = (Tree_float *)malloc(sizeof(Tree_float));
-    uint32_t i;
-    uint32_t *pidx;
+    uint64_t i;
+    uint64_t *pidx;
     float *bbox;
 
     tree->no_dims = no_dims;
 
     /* Initialize permutation array */
-    pidx = (uint32_t *)malloc(sizeof(uint32_t) * n);
+    pidx = (uint64_t *)malloc(sizeof(uint64_t) * n);
     for (i = 0; i < n; i++)
     {
         pidx[i] = i;
@@ -420,7 +420,7 @@ Create a tree node.
     start_idx : index of first data point to use
     n :  number of data points
 ************************************************/
-Node_float* create_node_float(uint32_t start_idx, uint32_t n, int is_leaf)
+Node_float* create_node_float(uint64_t start_idx, uint64_t n, int is_leaf)
 {
     Node_float *new_node;
     if (is_leaf)
@@ -566,11 +566,11 @@ Search a leaf node for closest point
     closest_idx : index of closest data point found (return)
     closest_dist : distance to closest point (return)
 ************************************************/
-void search_leaf_float(float *restrict pa, uint32_t *restrict pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, float *restrict point_coord,
-                 uint32_t k, uint32_t *restrict closest_idx, float *restrict closest_dist)
+void search_leaf_float(float *restrict pa, uint64_t *restrict pidx, int8_t no_dims, uint64_t start_idx, uint64_t n, float *restrict point_coord,
+                 uint64_t k, uint64_t *restrict closest_idx, float *restrict closest_dist)
 {
     float cur_dist;
-    uint32_t i;
+    uint64_t i;
     /* Loop through all points in leaf */
     for (i = 0; i < n; i++)
     {
@@ -598,11 +598,11 @@ Search a leaf node for closest point with data point mask
     closest_idx : index of closest data point found (return)
     closest_dist : distance to closest point (return)
 ************************************************/
-void search_leaf_float_mask(float *restrict pa, uint32_t *restrict pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, float *restrict point_coord,
-                               uint32_t k, uint8_t *mask, uint32_t *restrict closest_idx, float *restrict closest_dist)
+void search_leaf_float_mask(float *restrict pa, uint64_t *restrict pidx, int8_t no_dims, uint64_t start_idx, uint64_t n, float *restrict point_coord,
+                               uint64_t k, uint8_t *mask, uint64_t *restrict closest_idx, float *restrict closest_dist)
 {
     float cur_dist;
-    uint32_t i;
+    uint64_t i;
     /* Loop through all points in leaf */
     for (i = 0; i < n; i++)
     {
@@ -634,9 +634,9 @@ Search subtree for nearest to query point
     closest_idx : index of closest data point found (return)
     closest_dist : distance to closest point (return)
 ************************************************/
-void search_splitnode_float(Node_float *root, float *pa, uint32_t *pidx, int8_t no_dims, float *point_coord, 
-                      float min_dist, uint32_t k, float distance_upper_bound, float eps_fac, uint8_t *mask,
-                      uint32_t *closest_idx, float *closest_dist)
+void search_splitnode_float(Node_float *root, float *pa, uint64_t *pidx, int8_t no_dims, float *point_coord, 
+                      float min_dist, uint64_t k, float distance_upper_bound, float eps_fac, uint8_t *mask,
+                      uint64_t *closest_idx, float *closest_dist)
 {
     int8_t dim;
     float dist_left, dist_right;
@@ -733,21 +733,21 @@ Search for nearest neighbour for a set of query points
     closest_dist : distance to closest point (return)
 ************************************************/
 void search_tree_float(Tree_float *tree, float *pa, float *point_coords,
-                 uint32_t num_points, uint32_t k, float distance_upper_bound,
-                 float eps, uint8_t *mask, uint32_t *closest_idxs, float *closest_dists)
+                 uint64_t num_points, uint64_t k, float distance_upper_bound,
+                 float eps, uint8_t *mask, uint64_t *closest_idxs, float *closest_dists)
 {
     float min_dist;
     float eps_fac = 1 / ((1 + eps) * (1 + eps));
     int8_t no_dims = tree->no_dims;
     float *bbox = tree->bbox;
-    uint32_t *pidx = tree->pidx;
-    uint32_t j = 0;
+    uint64_t *pidx = tree->pidx;
+    uint64_t j = 0;
 #if defined(_MSC_VER) && defined(_OPENMP)
-    int32_t i = 0;
-    int32_t local_num_points = (int32_t) num_points;
+    int64_t i = 0;
+    int64_t local_num_points = (int64_t) num_points;
 #else
-    uint32_t i;
-    uint32_t local_num_points = num_points;
+    uint64_t i;
+    uint64_t local_num_points = num_points;
 #endif
     Node_float *root = (Node_float *)tree->root;
 
@@ -781,7 +781,7 @@ Insert point into priority queue
     cur_dist : distance to point inserted
     k : number of neighbours
 ************************************************/
-void insert_point_double(uint32_t *closest_idx, double *closest_dist, uint32_t pidx, double cur_dist, uint32_t k)
+void insert_point_double(uint64_t *closest_idx, double *closest_dist, uint64_t pidx, double cur_dist, uint64_t k)
 {
     int i;
     for (i = k - 1; i > 0; i--)
@@ -809,11 +809,11 @@ Get the bounding box of a set of points
     n : number of points
     bbox : bounding box (return)
 ************************************************/
-void get_bounding_box_double(double *pa, uint32_t *pidx, int8_t no_dims, uint32_t n, double *bbox)
+void get_bounding_box_double(double *pa, uint64_t *pidx, int8_t no_dims, uint64_t n, double *bbox)
 {
     double cur;
     int8_t i, j;
-    uint32_t bbox_idx, i2;
+    uint64_t bbox_idx, i2;
 
     /* Use first data point to initialize */
     for (i = 0; i < no_dims; i++)
@@ -854,12 +854,12 @@ The sliding midpoint rule is used for the partitioning.
     cut_val : value of cutting point (return)
     n_lo : number of point below cutting plane (return)
 ************************************************/
-int partition_double(double *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, double *bbox, int8_t *cut_dim, double *cut_val, uint32_t *n_lo)
+int partition_double(double *pa, uint64_t *pidx, int8_t no_dims, uint64_t start_idx, uint64_t n, double *bbox, int8_t *cut_dim, double *cut_val, uint64_t *n_lo)
 {
     int8_t dim = 0, i;
-    uint32_t p, q, i2;
+    uint64_t p, q, i2;
     double size = 0, min_val, max_val, split, side_len, cur_val;
-    uint32_t end_idx = start_idx + n - 1;
+    uint64_t end_idx = start_idx + n - 1;
 
     /* Find largest bounding box side */
     for (i = 0; i < no_dims; i++)
@@ -919,7 +919,7 @@ int partition_double(double *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_
            Minimum 1 point will be in lower box.
         */
 
-        uint32_t j = start_idx;
+        uint64_t j = start_idx;
         split = PA(j, dim);
         for (i2 = start_idx + 1; i2 <= end_idx; i2++)
         {
@@ -941,7 +941,7 @@ int partition_double(double *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_
            Minimum 1 point will be in higher box.
         */
 
-        uint32_t j = end_idx;
+        uint64_t j = end_idx;
         split = PA(j, dim);
         for (i2 = start_idx; i2 < end_idx; i2++)
         {
@@ -975,14 +975,14 @@ Construct a sub tree over a range of data points.
     bsp : number of points per leaf
     bbox : bounding box of set of data points
 ************************************************/
-Node_double* construct_subtree_double(double *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, uint32_t bsp, double *bbox)
+Node_double* construct_subtree_double(double *pa, uint64_t *pidx, int8_t no_dims, uint64_t start_idx, uint64_t n, uint64_t bsp, double *bbox)
 {
     /* Create new node */
     int is_leaf = (n <= bsp);
     Node_double *root = create_node_double(start_idx, n, is_leaf);
     int rval;
     int8_t cut_dim;
-    uint32_t n_lo;
+    uint64_t n_lo;
     double cut_val, lv, hv;
     if (is_leaf)
     {
@@ -1031,17 +1031,17 @@ Construct a tree over data points.
     n :  number of data points
     bsp : number of points per leaf
 ************************************************/
-Tree_double* construct_tree_double(double *pa, int8_t no_dims, uint32_t n, uint32_t bsp)
+Tree_double* construct_tree_double(double *pa, int8_t no_dims, uint64_t n, uint64_t bsp)
 {
     Tree_double *tree = (Tree_double *)malloc(sizeof(Tree_double));
-    uint32_t i;
-    uint32_t *pidx;
+    uint64_t i;
+    uint64_t *pidx;
     double *bbox;
 
     tree->no_dims = no_dims;
 
     /* Initialize permutation array */
-    pidx = (uint32_t *)malloc(sizeof(uint32_t) * n);
+    pidx = (uint64_t *)malloc(sizeof(uint64_t) * n);
     for (i = 0; i < n; i++)
     {
         pidx[i] = i;
@@ -1064,7 +1064,7 @@ Create a tree node.
     start_idx : index of first data point to use
     n :  number of data points
 ************************************************/
-Node_double* create_node_double(uint32_t start_idx, uint32_t n, int is_leaf)
+Node_double* create_node_double(uint64_t start_idx, uint64_t n, int is_leaf)
 {
     Node_double *new_node;
     if (is_leaf)
@@ -1210,11 +1210,11 @@ Search a leaf node for closest point
     closest_idx : index of closest data point found (return)
     closest_dist : distance to closest point (return)
 ************************************************/
-void search_leaf_double(double *restrict pa, uint32_t *restrict pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, double *restrict point_coord,
-                 uint32_t k, uint32_t *restrict closest_idx, double *restrict closest_dist)
+void search_leaf_double(double *restrict pa, uint64_t *restrict pidx, int8_t no_dims, uint64_t start_idx, uint64_t n, double *restrict point_coord,
+                 uint64_t k, uint64_t *restrict closest_idx, double *restrict closest_dist)
 {
     double cur_dist;
-    uint32_t i;
+    uint64_t i;
     /* Loop through all points in leaf */
     for (i = 0; i < n; i++)
     {
@@ -1242,11 +1242,11 @@ Search a leaf node for closest point with data point mask
     closest_idx : index of closest data point found (return)
     closest_dist : distance to closest point (return)
 ************************************************/
-void search_leaf_double_mask(double *restrict pa, uint32_t *restrict pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, double *restrict point_coord,
-                               uint32_t k, uint8_t *mask, uint32_t *restrict closest_idx, double *restrict closest_dist)
+void search_leaf_double_mask(double *restrict pa, uint64_t *restrict pidx, int8_t no_dims, uint64_t start_idx, uint64_t n, double *restrict point_coord,
+                               uint64_t k, uint8_t *mask, uint64_t *restrict closest_idx, double *restrict closest_dist)
 {
     double cur_dist;
-    uint32_t i;
+    uint64_t i;
     /* Loop through all points in leaf */
     for (i = 0; i < n; i++)
     {
@@ -1278,9 +1278,9 @@ Search subtree for nearest to query point
     closest_idx : index of closest data point found (return)
     closest_dist : distance to closest point (return)
 ************************************************/
-void search_splitnode_double(Node_double *root, double *pa, uint32_t *pidx, int8_t no_dims, double *point_coord, 
-                      double min_dist, uint32_t k, double distance_upper_bound, double eps_fac, uint8_t *mask,
-                      uint32_t *closest_idx, double *closest_dist)
+void search_splitnode_double(Node_double *root, double *pa, uint64_t *pidx, int8_t no_dims, double *point_coord, 
+                      double min_dist, uint64_t k, double distance_upper_bound, double eps_fac, uint8_t *mask,
+                      uint64_t *closest_idx, double *closest_dist)
 {
     int8_t dim;
     double dist_left, dist_right;
@@ -1377,21 +1377,21 @@ Search for nearest neighbour for a set of query points
     closest_dist : distance to closest point (return)
 ************************************************/
 void search_tree_double(Tree_double *tree, double *pa, double *point_coords,
-                 uint32_t num_points, uint32_t k, double distance_upper_bound,
-                 double eps, uint8_t *mask, uint32_t *closest_idxs, double *closest_dists)
+                 uint64_t num_points, uint64_t k, double distance_upper_bound,
+                 double eps, uint8_t *mask, uint64_t *closest_idxs, double *closest_dists)
 {
     double min_dist;
     double eps_fac = 1 / ((1 + eps) * (1 + eps));
     int8_t no_dims = tree->no_dims;
     double *bbox = tree->bbox;
-    uint32_t *pidx = tree->pidx;
-    uint32_t j = 0;
+    uint64_t *pidx = tree->pidx;
+    uint64_t j = 0;
 #if defined(_MSC_VER) && defined(_OPENMP)
-    int32_t i = 0;
-    int32_t local_num_points = (int32_t) num_points;
+    int64_t i = 0;
+    int64_t local_num_points = (int64_t) num_points;
 #else
-    uint32_t i;
-    uint32_t local_num_points = num_points;
+    uint64_t i;
+    uint64_t local_num_points = num_points;
 #endif
     Node_double *root = (Node_double *)tree->root;
 
diff --git a/pykdtree/_kdtree_core.c.mako b/pykdtree/_kdtree_core.c.mako
index 6e3eb8b..0ff0ba4 100644
--- a/pykdtree/_kdtree_core.c.mako
+++ b/pykdtree/_kdtree_core.c.mako
@@ -29,7 +29,7 @@ Anne M. Archibald and libANN by David M. Mount and Sunil Arya.
 #include <float.h>
 
 #define PA(i,d)			(pa[no_dims * pidx[i] + d])
-#define PASWAP(a,b) { uint32_t tmp = pidx[a]; pidx[a] = pidx[b]; pidx[b] = tmp; }
+#define PASWAP(a,b) { uint64_t tmp = pidx[a]; pidx[a] = pidx[b]; pidx[b] = tmp; }
 
 #ifdef _MSC_VER
 #define restrict __restrict
@@ -41,8 +41,8 @@ typedef struct
 {
     ${DTYPE} cut_val;
     int8_t cut_dim;
-    uint32_t start_idx;
-    uint32_t n;
+    uint64_t start_idx;
+    uint64_t n;
     ${DTYPE} cut_bounds_lv;
     ${DTYPE} cut_bounds_hv;
     struct Node_${DTYPE} *left_child;
@@ -53,7 +53,7 @@ typedef struct
 {
     ${DTYPE} *bbox;
     int8_t no_dims;
-    uint32_t *pidx;
+    uint64_t *pidx;
     struct Node_${DTYPE} *root;
 } Tree_${DTYPE};
 
@@ -61,28 +61,28 @@ typedef struct
 
 % for DTYPE in ['float', 'double']:
 
-void insert_point_${DTYPE}(uint32_t *closest_idx, ${DTYPE} *closest_dist, uint32_t pidx, ${DTYPE} cur_dist, uint32_t k);
-void get_bounding_box_${DTYPE}(${DTYPE} *pa, uint32_t *pidx, int8_t no_dims, uint32_t n, ${DTYPE} *bbox);
-int partition_${DTYPE}(${DTYPE} *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, ${DTYPE} *bbox, int8_t *cut_dim,
-              ${DTYPE} *cut_val, uint32_t *n_lo);
-Tree_${DTYPE}* construct_tree_${DTYPE}(${DTYPE} *pa, int8_t no_dims, uint32_t n, uint32_t bsp);
-Node_${DTYPE}* construct_subtree_${DTYPE}(${DTYPE} *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, uint32_t bsp, ${DTYPE} *bbox);
-Node_${DTYPE} * create_node_${DTYPE}(uint32_t start_idx, uint32_t n, int is_leaf);
+void insert_point_${DTYPE}(uint64_t *closest_idx, ${DTYPE} *closest_dist, uint64_t pidx, ${DTYPE} cur_dist, uint64_t k);
+void get_bounding_box_${DTYPE}(${DTYPE} *pa, uint64_t *pidx, int8_t no_dims, uint64_t n, ${DTYPE} *bbox);
+int partition_${DTYPE}(${DTYPE} *pa, uint64_t *pidx, int8_t no_dims, uint64_t start_idx, uint64_t n, ${DTYPE} *bbox, int8_t *cut_dim,
+              ${DTYPE} *cut_val, uint64_t *n_lo);
+Tree_${DTYPE}* construct_tree_${DTYPE}(${DTYPE} *pa, int8_t no_dims, uint64_t n, uint64_t bsp);
+Node_${DTYPE}* construct_subtree_${DTYPE}(${DTYPE} *pa, uint64_t *pidx, int8_t no_dims, uint64_t start_idx, uint64_t n, uint64_t bsp, ${DTYPE} *bbox);
+Node_${DTYPE} * create_node_${DTYPE}(uint64_t start_idx, uint64_t n, int is_leaf);
 void delete_subtree_${DTYPE}(Node_${DTYPE} *root);
 void delete_tree_${DTYPE}(Tree_${DTYPE} *tree);
 void print_tree_${DTYPE}(Node_${DTYPE} *root, int level);
 ${DTYPE} calc_dist_${DTYPE}(${DTYPE} *point1_coord, ${DTYPE} *point2_coord, int8_t no_dims);
 ${DTYPE} get_cube_offset_${DTYPE}(int8_t dim, ${DTYPE} *point_coord, ${DTYPE} *bbox);
 ${DTYPE} get_min_dist_${DTYPE}(${DTYPE} *point_coord, int8_t no_dims, ${DTYPE} *bbox);
-void search_leaf_${DTYPE}(${DTYPE} *restrict pa, uint32_t *restrict pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, ${DTYPE} *restrict point_coord,
-                 uint32_t k, uint32_t *restrict closest_idx, ${DTYPE} *restrict closest_dist);
-void search_leaf_${DTYPE}_mask(${DTYPE} *restrict pa, uint32_t *restrict pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, ${DTYPE} *restrict point_coord,
-                 uint32_t k, uint8_t *restrict mask, uint32_t *restrict closest_idx, ${DTYPE} *restrict closest_dist);
-void search_splitnode_${DTYPE}(Node_${DTYPE} *root, ${DTYPE} *pa, uint32_t *pidx, int8_t no_dims, ${DTYPE} *point_coord,
-                      ${DTYPE} min_dist, uint32_t k, ${DTYPE} distance_upper_bound, ${DTYPE} eps_fac, uint8_t *mask, uint32_t *  closest_idx, ${DTYPE} *closest_dist);
+void search_leaf_${DTYPE}(${DTYPE} *restrict pa, uint64_t *restrict pidx, int8_t no_dims, uint64_t start_idx, uint64_t n, ${DTYPE} *restrict point_coord,
+                 uint64_t k, uint64_t *restrict closest_idx, ${DTYPE} *restrict closest_dist);
+void search_leaf_${DTYPE}_mask(${DTYPE} *restrict pa, uint64_t *restrict pidx, int8_t no_dims, uint64_t start_idx, uint64_t n, ${DTYPE} *restrict point_coord,
+                 uint64_t k, uint8_t *restrict mask, uint64_t *restrict closest_idx, ${DTYPE} *restrict closest_dist);
+void search_splitnode_${DTYPE}(Node_${DTYPE} *root, ${DTYPE} *pa, uint64_t *pidx, int8_t no_dims, ${DTYPE} *point_coord,
+                      ${DTYPE} min_dist, uint64_t k, ${DTYPE} distance_upper_bound, ${DTYPE} eps_fac, uint8_t *mask, uint64_t *  closest_idx, ${DTYPE} *closest_dist);
 void search_tree_${DTYPE}(Tree_${DTYPE} *tree, ${DTYPE} *pa, ${DTYPE} *point_coords,
-                 uint32_t num_points, uint32_t k,  ${DTYPE} distance_upper_bound,
-                 ${DTYPE} eps, uint8_t *mask, uint32_t *closest_idxs, ${DTYPE} *closest_dists);
+                 uint64_t num_points, uint64_t k,  ${DTYPE} distance_upper_bound,
+                 ${DTYPE} eps, uint8_t *mask, uint64_t *closest_idxs, ${DTYPE} *closest_dists);
 
 % endfor
 
@@ -97,7 +97,7 @@ Params:
     cur_dist : distance to point inserted
     k : number of neighbours
 ************************************************/
-void insert_point_${DTYPE}(uint32_t *closest_idx, ${DTYPE} *closest_dist, uint32_t pidx, ${DTYPE} cur_dist, uint32_t k)
+void insert_point_${DTYPE}(uint64_t *closest_idx, ${DTYPE} *closest_dist, uint64_t pidx, ${DTYPE} cur_dist, uint64_t k)
 {
     int i;
     for (i = k - 1; i > 0; i--)
@@ -125,11 +125,11 @@ Params:
     n : number of points
     bbox : bounding box (return)
 ************************************************/
-void get_bounding_box_${DTYPE}(${DTYPE} *pa, uint32_t *pidx, int8_t no_dims, uint32_t n, ${DTYPE} *bbox)
+void get_bounding_box_${DTYPE}(${DTYPE} *pa, uint64_t *pidx, int8_t no_dims, uint64_t n, ${DTYPE} *bbox)
 {
     ${DTYPE} cur;
     int8_t i, j;
-    uint32_t bbox_idx, i2;
+    uint64_t bbox_idx, i2;
 
     /* Use first data point to initialize */
     for (i = 0; i < no_dims; i++)
@@ -170,12 +170,12 @@ Params:
     cut_val : value of cutting point (return)
     n_lo : number of point below cutting plane (return)
 ************************************************/
-int partition_${DTYPE}(${DTYPE} *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, ${DTYPE} *bbox, int8_t *cut_dim, ${DTYPE} *cut_val, uint32_t *n_lo)
+int partition_${DTYPE}(${DTYPE} *pa, uint64_t *pidx, int8_t no_dims, uint64_t start_idx, uint64_t n, ${DTYPE} *bbox, int8_t *cut_dim, ${DTYPE} *cut_val, uint64_t *n_lo)
 {
     int8_t dim = 0, i;
-    uint32_t p, q, i2;
+    uint64_t p, q, i2;
     ${DTYPE} size = 0, min_val, max_val, split, side_len, cur_val;
-    uint32_t end_idx = start_idx + n - 1;
+    uint64_t end_idx = start_idx + n - 1;
 
     /* Find largest bounding box side */
     for (i = 0; i < no_dims; i++)
@@ -235,7 +235,7 @@ int partition_${DTYPE}(${DTYPE} *pa, uint32_t *pidx, int8_t no_dims, uint32_t st
            Minimum 1 point will be in lower box.
         */
 
-        uint32_t j = start_idx;
+        uint64_t j = start_idx;
         split = PA(j, dim);
         for (i2 = start_idx + 1; i2 <= end_idx; i2++)
         {
@@ -257,7 +257,7 @@ int partition_${DTYPE}(${DTYPE} *pa, uint32_t *pidx, int8_t no_dims, uint32_t st
            Minimum 1 point will be in higher box.
         */
 
-        uint32_t j = end_idx;
+        uint64_t j = end_idx;
         split = PA(j, dim);
         for (i2 = start_idx; i2 < end_idx; i2++)
         {
@@ -291,14 +291,14 @@ Params:
     bsp : number of points per leaf
     bbox : bounding box of set of data points
 ************************************************/
-Node_${DTYPE}* construct_subtree_${DTYPE}(${DTYPE} *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, uint32_t bsp, ${DTYPE} *bbox)
+Node_${DTYPE}* construct_subtree_${DTYPE}(${DTYPE} *pa, uint64_t *pidx, int8_t no_dims, uint64_t start_idx, uint64_t n, uint64_t bsp, ${DTYPE} *bbox)
 {
     /* Create new node */
     int is_leaf = (n <= bsp);
     Node_${DTYPE} *root = create_node_${DTYPE}(start_idx, n, is_leaf);
     int rval;
     int8_t cut_dim;
-    uint32_t n_lo;
+    uint64_t n_lo;
     ${DTYPE} cut_val, lv, hv;
     if (is_leaf)
     {
@@ -347,17 +347,17 @@ Params:
     n :  number of data points
     bsp : number of points per leaf
 ************************************************/
-Tree_${DTYPE}* construct_tree_${DTYPE}(${DTYPE} *pa, int8_t no_dims, uint32_t n, uint32_t bsp)
+Tree_${DTYPE}* construct_tree_${DTYPE}(${DTYPE} *pa, int8_t no_dims, uint64_t n, uint64_t bsp)
 {
     Tree_${DTYPE} *tree = (Tree_${DTYPE} *)malloc(sizeof(Tree_${DTYPE}));
-    uint32_t i;
-    uint32_t *pidx;
+    uint64_t i;
+    uint64_t *pidx;
     ${DTYPE} *bbox;
 
     tree->no_dims = no_dims;
 
     /* Initialize permutation array */
-    pidx = (uint32_t *)malloc(sizeof(uint32_t) * n);
+    pidx = (uint64_t *)malloc(sizeof(uint64_t) * n);
     for (i = 0; i < n; i++)
     {
         pidx[i] = i;
@@ -380,7 +380,7 @@ Params:
     start_idx : index of first data point to use
     n :  number of data points
 ************************************************/
-Node_${DTYPE}* create_node_${DTYPE}(uint32_t start_idx, uint32_t n, int is_leaf)
+Node_${DTYPE}* create_node_${DTYPE}(uint64_t start_idx, uint64_t n, int is_leaf)
 {
     Node_${DTYPE} *new_node;
     if (is_leaf)
@@ -526,11 +526,11 @@ Params:
     closest_idx : index of closest data point found (return)
     closest_dist : distance to closest point (return)
 ************************************************/
-void search_leaf_${DTYPE}(${DTYPE} *restrict pa, uint32_t *restrict pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, ${DTYPE} *restrict point_coord,
-                 uint32_t k, uint32_t *restrict closest_idx, ${DTYPE} *restrict closest_dist)
+void search_leaf_${DTYPE}(${DTYPE} *restrict pa, uint64_t *restrict pidx, int8_t no_dims, uint64_t start_idx, uint64_t n, ${DTYPE} *restrict point_coord,
+                 uint64_t k, uint64_t *restrict closest_idx, ${DTYPE} *restrict closest_dist)
 {
     ${DTYPE} cur_dist;
-    uint32_t i;
+    uint64_t i;
     /* Loop through all points in leaf */
     for (i = 0; i < n; i++)
     {
@@ -558,11 +558,11 @@ Params:
     closest_idx : index of closest data point found (return)
     closest_dist : distance to closest point (return)
 ************************************************/
-void search_leaf_${DTYPE}_mask(${DTYPE} *restrict pa, uint32_t *restrict pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, ${DTYPE} *restrict point_coord,
-                               uint32_t k, uint8_t *mask, uint32_t *restrict closest_idx, ${DTYPE} *restrict closest_dist)
+void search_leaf_${DTYPE}_mask(${DTYPE} *restrict pa, uint64_t *restrict pidx, int8_t no_dims, uint64_t start_idx, uint64_t n, ${DTYPE} *restrict point_coord,
+                               uint64_t k, uint8_t *mask, uint64_t *restrict closest_idx, ${DTYPE} *restrict closest_dist)
 {
     ${DTYPE} cur_dist;
-    uint32_t i;
+    uint64_t i;
     /* Loop through all points in leaf */
     for (i = 0; i < n; i++)
     {
@@ -594,9 +594,9 @@ Params:
     closest_idx : index of closest data point found (return)
     closest_dist : distance to closest point (return)
 ************************************************/
-void search_splitnode_${DTYPE}(Node_${DTYPE} *root, ${DTYPE} *pa, uint32_t *pidx, int8_t no_dims, ${DTYPE} *point_coord, 
-                      ${DTYPE} min_dist, uint32_t k, ${DTYPE} distance_upper_bound, ${DTYPE} eps_fac, uint8_t *mask,
-                      uint32_t *closest_idx, ${DTYPE} *closest_dist)
+void search_splitnode_${DTYPE}(Node_${DTYPE} *root, ${DTYPE} *pa, uint64_t *pidx, int8_t no_dims, ${DTYPE} *point_coord, 
+                      ${DTYPE} min_dist, uint64_t k, ${DTYPE} distance_upper_bound, ${DTYPE} eps_fac, uint8_t *mask,
+                      uint64_t *closest_idx, ${DTYPE} *closest_dist)
 {
     int8_t dim;
     ${DTYPE} dist_left, dist_right;
@@ -693,21 +693,21 @@ Params:
     closest_dist : distance to closest point (return)
 ************************************************/
 void search_tree_${DTYPE}(Tree_${DTYPE} *tree, ${DTYPE} *pa, ${DTYPE} *point_coords,
-                 uint32_t num_points, uint32_t k, ${DTYPE} distance_upper_bound,
-                 ${DTYPE} eps, uint8_t *mask, uint32_t *closest_idxs, ${DTYPE} *closest_dists)
+                 uint64_t num_points, uint64_t k, ${DTYPE} distance_upper_bound,
+                 ${DTYPE} eps, uint8_t *mask, uint64_t *closest_idxs, ${DTYPE} *closest_dists)
 {
     ${DTYPE} min_dist;
     ${DTYPE} eps_fac = 1 / ((1 + eps) * (1 + eps));
     int8_t no_dims = tree->no_dims;
     ${DTYPE} *bbox = tree->bbox;
-    uint32_t *pidx = tree->pidx;
-    uint32_t j = 0;
+    uint64_t *pidx = tree->pidx;
+    uint64_t j = 0;
 #if defined(_MSC_VER) && defined(_OPENMP)
-    int32_t i = 0;
-    int32_t local_num_points = (int32_t) num_points;
+    int64_t i = 0;
+    int64_t local_num_points = (int64_t) num_points;
 #else
-    uint32_t i;
-    uint32_t local_num_points = num_points;
+    uint64_t i;
+    uint64_t local_num_points = num_points;
 #endif
     Node_${DTYPE} *root = (Node_${DTYPE} *)tree->root;
 
diff --git a/pykdtree/kdtree.pyx b/pykdtree/kdtree.pyx
index 7104810..8c42485 100644
--- a/pykdtree/kdtree.pyx
+++ b/pykdtree/kdtree.pyx
@@ -17,7 +17,7 @@
 
 import numpy as np
 cimport numpy as np
-from libc.stdint cimport uint32_t, int8_t, uint8_t
+from libc.stdint cimport uint64_t, int8_t, uint8_t
 cimport cython
 
 np.import_array()
@@ -26,8 +26,8 @@ np.import_array()
 cdef struct node_float:
     float cut_val
     int8_t cut_dim
-    uint32_t start_idx
-    uint32_t n
+    uint64_t start_idx
+    uint64_t n
     float cut_bounds_lv
     float cut_bounds_hv
     node_float *left_child
@@ -36,14 +36,14 @@ cdef struct node_float:
 cdef struct tree_float:
     float *bbox
     int8_t no_dims
-    uint32_t *pidx
+    uint64_t *pidx
     node_float *root
 
 cdef struct node_double:
     double cut_val
     int8_t cut_dim
-    uint32_t start_idx
-    uint32_t n
+    uint64_t start_idx
+    uint64_t n
     double cut_bounds_lv
     double cut_bounds_hv
     node_double *left_child
@@ -52,15 +52,15 @@ cdef struct node_double:
 cdef struct tree_double:
     double *bbox
     int8_t no_dims
-    uint32_t *pidx
+    uint64_t *pidx
     node_double *root
 
-cdef extern tree_float* construct_tree_float(float *pa, int8_t no_dims, uint32_t n, uint32_t bsp) nogil
-cdef extern void search_tree_float(tree_float *kdtree, float *pa, float *point_coords, uint32_t num_points, uint32_t k, float distance_upper_bound, float eps_fac, uint8_t *mask, uint32_t *closest_idxs, float *closest_dists) nogil
+cdef extern tree_float* construct_tree_float(float *pa, int8_t no_dims, uint64_t n, uint64_t bsp) nogil
+cdef extern void search_tree_float(tree_float *kdtree, float *pa, float *point_coords, uint64_t num_points, uint64_t k, float distance_upper_bound, float eps_fac, uint8_t *mask, uint64_t *closest_idxs, float *closest_dists) nogil
 cdef extern void delete_tree_float(tree_float *kdtree)
 
-cdef extern tree_double* construct_tree_double(double *pa, int8_t no_dims, uint32_t n, uint32_t bsp) nogil
-cdef extern void search_tree_double(tree_double *kdtree, double *pa, double *point_coords, uint32_t num_points, uint32_t k, double distance_upper_bound, double eps_fac, uint8_t *mask, uint32_t *closest_idxs, double *closest_dists) nogil
+cdef extern tree_double* construct_tree_double(double *pa, int8_t no_dims, uint64_t n, uint64_t bsp) nogil
+cdef extern void search_tree_double(tree_double *kdtree, double *pa, double *point_coords, uint64_t num_points, uint64_t k, double distance_upper_bound, double eps_fac, uint8_t *mask, uint64_t *closest_idxs, double *closest_dists) nogil
 cdef extern void delete_tree_double(tree_double *kdtree)
 
 cdef class KDTree:
@@ -81,9 +81,9 @@ cdef class KDTree:
     cdef readonly np.ndarray data
     cdef float *_data_pts_data_float
     cdef double *_data_pts_data_double
-    cdef readonly uint32_t n
+    cdef readonly uint64_t n
     cdef readonly int8_t ndim
-    cdef readonly uint32_t leafsize
+    cdef readonly uint64_t leafsize
 
     def __cinit__(KDTree self):
         self._kdtree_float = NULL
@@ -116,8 +116,8 @@ cdef class KDTree:
         self.data = self.data_pts
 
         # Get tree info
-        self.n = <uint32_t>data_pts.shape[0]
-        self.leafsize = <uint32_t>leafsize
+        self.n = <uint64_t>data_pts.shape[0]
+        self.leafsize = <uint64_t>leafsize
         if data_pts.ndim == 1:
             self.ndim = 1
         elif data_pts.shape[1] > 127:
@@ -185,15 +185,15 @@ cdef class KDTree:
             raise TypeError('Type mismatch. query points must be of type float32 when data points are of type float32')
 
         # Get query info
-        cdef uint32_t num_qpoints = query_pts.shape[0]
-        cdef uint32_t num_n = k
-        cdef np.ndarray[uint32_t, ndim=1] closest_idxs = np.empty(num_qpoints * k, dtype=np.uint32)
+        cdef uint64_t num_qpoints = query_pts.shape[0]
+        cdef uint64_t num_n = k
+        cdef np.ndarray[uint64_t, ndim=1] closest_idxs = np.empty(num_qpoints * k, dtype=np.uint64)
         cdef np.ndarray[float, ndim=1] closest_dists_float
         cdef np.ndarray[double, ndim=1] closest_dists_double
 
 
         # Set up return arrays
-        cdef uint32_t *closest_idxs_data = <uint32_t *>closest_idxs.data
+        cdef uint64_t *closest_idxs_data = <uint64_t *>closest_idxs.data
         cdef float *closest_dists_data_float
         cdef double *closest_dists_data_double