                                                              scratch_d + offsets[l],
                                                              levels[l]);
       }
+      
+      // Now, level L-1 is processed, scan levels L-2...0 pulling start offsets from the level above.
+      for (int i = (int)L - 2; 0 <= i; i--) {
+        scan<false, false, false> << <levels[i], blockSize, 0, stream >> > (scratch_d + offsets[l],
+                                                                            nullptr,
+                                                                            nullptr,
+                                                                            scratch_d + offsets[l],
+                                                                            levels[l]);
+      }
+      // Now, level 0 is processed, scan input writing to output, pulling offsets from level 0.
+
+      scan<false, false, false> << <levels[0], blockSize, 0, stream >> > (offsets_d,
+                                                                          nullptr,
+                                                                          nullptr,
+                                                                          counts_d,
+                                                                          n);
 
 
     }

                                     size_t N,
                                     cudaStream_t stream)
{
  std::vector<size_t> levels;
  calcLevelSizes(levels, N);

  std::vector<size_t> offsets;
  offsets.push_back(0);
  size_t alignment = 128 / sizeof(uint32_t);
  for (size_t i = 0; i < levels.size(); i++) {
    size_t levelSize = (levels[i] + alignment - 1) & ~alignment;
    offsets.push_back(offsets[i] + levelSize);
  }

  if (N <= std::numeric_limits<uint32_t>::max()) {
    uint32_t blockSize = SCAN_WARPS * WARPSIZE;
    uint32_t n = static_cast<uint32_t>(N);

    if (levels.empty()) {


#include <limits>
#include <vector>
#include <algorithm>
#include "Scan.h"

#define WARPSIZE 32
