  * are identical, then the longest prefix of i and j is used as a tie breaker.
  */
 __device__ inline int delta(const uint32_t *d_morton_codes,
-                            unsigned int i,
-                            unsigned int j,
-                            int min_idx,
-                            int max_idx)
-    {
+                            const int i,
+                            const int j,
+                            const int min_idx,
+                            const int max_idx)
+    {
+    assert(i >= min_idx && i <= max_idx);
+
     if (j > max_idx || j < min_idx)
         {
         return -1;

 * \param d_morton_codes Array of Morton codes
 * \param min_idx The smallest Morton code index considered "in range" (inclusive)
 * \param max_idx The last Morton code index considered "in range" (inclusive)
 * \param idx Current node (Morton code) index
 *
 * \returns the minimum and maximum leafs covered by this node
 * \note This is a literal implementation of the Karras pseudocode, with no optimizations or refinement.
 *       Tero Karras, "Maximizing parallelism in the construction of BVHs, octrees, and k-d trees",
 *       High Performance Graphics (2012).
 */
__device__ inline uint2 determineRange(const uint32_t *d_morton_codes,
                                       const int min_idx,
                                       const int max_idx,
                                       const int idx)
    {
    int forward_prefix = delta(d_morton_codes, idx, idx+1, min_idx, max_idx);
    int backward_prefix = delta(d_morton_codes, idx, idx-1, min_idx, max_idx);

    // get direction of the range based on sign
    int d = ((forward_prefix - backward_prefix) > 0) ? 1 : -1;

    // get minimum prefix
    int min_prefix = delta(d_morton_codes, idx, idx-d, min_idx, max_idx);

    // get maximum prefix by binary search
    int lmax = 2;
    while( delta(d_morton_codes, idx, idx + d*lmax, min_idx, max_idx) > min_prefix)
        {
        lmax = lmax << 1;
        }

    unsigned int len = 0;
    unsigned int step = lmax;
    do
        {
        step = step >> 1;
        unsigned int new_len = len + step;
        if (delta(d_morton_codes, idx, idx + d*new_len, min_idx, max_idx) > min_prefix)
            len = new_len;
        }
    while (step > 1);

   // order range based on direction
    uint2 range;
    if (d > 0)
        {
        range.x = idx;
        range.y = idx + len;
        }
    else
        {
        range.x = idx - len;
        range.y = idx;
        }
    return range;
    }

//! Finds the split position in Morton codes covered by a range
/*!
 * \param d_morton_codes Array of Morton codes
 * \param first First leaf node in the range
 * \param last Last leaf node in the range
 *
 * \returns the leaf index corresponding to the split in Morton codes
 * See determineRange for original source of algorithm.
 */
__device__ inline unsigned int findSplit(const uint32_t *d_morton_codes,
                                         const unsigned int first,
                                         const unsigned int last)
    {
    uint32_t first_code = d_morton_codes[first];
    uint32_t last_code = d_morton_codes[last];

    // if codes match, then just split evenly
    if (first_code == last_code)
        return (first + last) >> 1;

    // get the length of the common prefix
    int common_prefix = __clz(first_code ^ last_code);

    // assume split starts at first, and begin binary search
    unsigned int split = first;
    unsigned int step = last - first;
    do
        {
        // exponential decrease (is factor of 2 best?)
        step = (step + 1) >> 1;
        unsigned int new_split = split + step;

        // if proposed split lies within range
        if (new_split < last)
            {
            unsigned int split_code = d_morton_codes[new_split];
            int split_prefix = __clz(first_code ^ split_code);

            // if new split shares a longer number of bits, accept it
            if (split_prefix > common_prefix)
                {
                split = new_split;
                }
            }
        }
    while (step > 1);

    return split;
    }

//! Kernel to generate the parent-child-sibling relationships between nodes

        }

    // enact the magical split determining
    uint3 range_split = determineRangeSplit(d_morton_codes, origin, end, node_idx);
    unsigned int first = range_split.x;
    unsigned int last = range_split.y;
    unsigned int split = range_split.z;

    uint2 children;
    // set the children, shifting ahead by nleafs - cur_type to account for leaf shifting
