diff --git a/.envrc b/.envrc
index f59b4d9..32ae5f6 100644
--- a/.envrc
+++ b/.envrc
@@ -1,5 +1,5 @@
 if ! has nix_direnv_version || ! nix_direnv_version 3.0.4; then
     source_url "https://raw.githubusercontent.com/nix-community/nix-direnv/3.0.4/direnvrc" "sha256-DzlYZ33mWF/Gs8DDeyjr8mnVmQGx7ASYqA5WlxwvBG4="
 fi
-watch_file devShell.nix shell.nix flake.nix
+watch_file shell.nix flake.nix
 use flake || use nix
diff --git a/Makefile b/Makefile
index fb6f455..db623e9 100644
--- a/Makefile
+++ b/Makefile
@@ -5,14 +5,15 @@ SHARED_LIB = libsparsemap.so
 
 #CFLAGS = -Wall -Wextra -Wpedantic -Of -std=c11 -Iinclude/ -fPIC
 #CFLAGS = -Wall -Wextra -Wpedantic -Og -g -std=c11 -Iinclude/ -fPIC
-CFLAGS = -DSPARSEMAP_DIAGNOSTIC -DSPARSEMAP_ASSERT -Wall -Wextra -Wpedantic -Og -g -std=c11 -Iinclude/ -fPIC
-#CFLAGS = -Wall -Wextra -Wpedantic -Og -g -fsanitize=address,leak,object-size,pointer-compare,pointer-subtract,null,return,bounds,pointer-overflow,undefined -fsanitize-address-use-after-scope -std=c11 -Iinclude/ -fPIC
+CFLAGS = -DSPARSEMAP_DIAGNOSTIC -DDEBUG -Wall -Wextra -Wpedantic -Og -g -std=c11 -Iinclude/ -fPIC
+#CFLAGS = -DSPARSEMAP_DIAGNOSTIC -DDEBUG -Wall -Wextra -Wpedantic -Og -g -fsanitize=address,leak,object-size,pointer-compare,pointer-subtract,null,return,bounds,pointer-overflow,undefined -fsanitize-address-use-after-scope -std=c11 -Iinclude/ -fPIC
 #CFLAGS = -Wall -Wextra -Wpedantic -Og -g -fsanitize=all -fhardened -std=c11 -Iinclude/ -fPIC
 
-TEST_FLAGS = -Wall -Wextra -Wpedantic -Og -g -std=c11 -Iinclude/ -Itests/ -fPIC
+TEST_FLAGS = -DDEBUG -Wall -Wextra -Wpedantic -Og -g -std=c11 -Iinclude/ -Itests/ -fPIC
+#TEST_FLAGS = -DDEBUG -Wall -Wextra -Wpedantic -Og -g -fsanitize=address,leak,object-size,pointer-compare,pointer-subtract,null,return,bounds,pointer-overflow,undefined -fsanitize-address-use-after-scope -std=c11 -Iinclude/ -fPIC
 
 TESTS = tests/test
-TEST_OBJS = tests/test.o tests/munit.o tests/common.o
+TEST_OBJS = tests/test.o tests/munit.o tests/tdigest.o tests/common.o
 EXAMPLES = examples/ex_1 examples/ex_2 examples/ex_3 examples/ex_4
 
 .PHONY: all shared static clean test examples mls
@@ -39,7 +40,7 @@ check: test
 	env ASAN_OPTIONS=detect_leaks=1 LSAN_OPTIONS=verbosity=1:log_threads=1 ./tests/test
 
 tests/test: $(TEST_OBJS) $(STATIC_LIB)
-	$(CC) $^ -o $@ $(TEST_FLAGS)
+	$(CC) $^ -lm -o $@ $(TEST_FLAGS)
 
 clean:
 	rm -f $(OBJS)
@@ -76,4 +77,7 @@ examples/ex_3: examples/common.o examples/ex_3.o $(STATIC_LIB)
 examples/ex_4: examples/common.o examples/ex_4.o $(STATIC_LIB)
 	$(CC) $^ -o $@ $(CFLAGS) $(TEST_FLAGS)
 
+todo:
+	rg -i 'todo|gsb|abort'
+
 # cp src/sparsemap.c /tmp && clang-tidy src/sparsemap.c -fix -fix-errors -checks="readability-braces-around-statements" -- -DDEBUG -DSPARSEMAP_DIAGNOSTIC -DSPARSEMAP_ASSERT -Wall -Wextra -Wpedantic -Og -g -std=c11 -Iinclude/ -fPIC
diff --git a/README.md b/README.md
index 8bb7acd..2e91716 100644
--- a/README.md
+++ b/README.md
@@ -60,5 +60,6 @@ to an uncompressed bit vector (sometimes higher due to the bytes required for
 metadata). In such cases, other compression schemes are more efficient (i.e.
 http://lemire.me/blog/archives/2008/08/20/the-mythical-bitmap-index/).
 
-This library was originally created for hamsterdb [http://hamsterdb.com] in
-C++ and then translated to C99 code by Greg Burd <greg@burd.me>.
+This library was originally created for [hamsterdb](http://hamsterdb.com) in
+C++ and then translated to C and further improved by Greg Burd <greg@burd.me>
+for use in LMDB and OpenLDAP.
diff --git a/examples/ex_1.c b/examples/ex_1.c
index 5fa1ae4..dea68f1 100644
--- a/examples/ex_1.c
+++ b/examples/ex_1.c
@@ -1,4 +1,7 @@
 #include <assert.h>
+#include <stdbool.h>
+#include <stddef.h>
+#include <stdint.h>
 #include <stdio.h>
 
 #include "../include/sparsemap.h"
@@ -14,16 +17,18 @@
 
 /* !!! Duplicated here for testing purposes. Keep in sync, or suffer. !!! */
 struct sparsemap {
-  uint8_t *m_data;
   size_t m_capacity;
   size_t m_data_used;
+  uint8_t *m_data;
 };
 
 int
 main()
 {
   size_t size = 4;
-  setbuf(stderr, 0); // disable buffering
+  setvbuf(stdout, NULL, _IONBF, 0); // Disable buffering for stdout
+  setvbuf(stderr, NULL, _IONBF, 0); // Disable buffering for stdout
+
   __diag("Please wait a moment...");
   sparsemap_t mmap, *map = &mmap;
   uint8_t buffer[1024];
@@ -135,7 +140,7 @@ main()
     sparsemap_set(map, i, true);
   }
   for (int i = 0; i < 100000; i++) {
-    assert(sparsemap_select(map, i) == (unsigned)i);
+    assert(sparsemap_select(map, i, true) == (unsigned)i);
   }
 
   sparsemap_clear(map);
@@ -145,7 +150,7 @@ main()
     sparsemap_set(map, i, true);
   }
   for (int i = 1; i < 513; i++) {
-    assert(sparsemap_select(map, i - 1) == (unsigned)i);
+    assert(sparsemap_select(map, i - 1, true) == (unsigned)i);
   }
 
   sparsemap_clear(map);
@@ -155,7 +160,7 @@ main()
     sparsemap_set(map, i * 10, true);
   }
   for (size_t i = 0; i < 8; i++) {
-    assert(sparsemap_select(map, i) == i * 10);
+    assert(sparsemap_select(map, i, true) == (sparsemap_idx_t)i * 10);
   }
 
   // split and move, aligned to MiniMap capacity
diff --git a/examples/ex_2.c b/examples/ex_2.c
index 97f1c82..d2250b6 100644
--- a/examples/ex_2.c
+++ b/examples/ex_2.c
@@ -1,9 +1,7 @@
 #include <assert.h>
-#include <stdarg.h>
+#include <stdbool.h>
 #include <stdio.h>
 #include <stdlib.h>
-#include <time.h>
-#include <unistd.h>
 
 #include "../include/sparsemap.h"
 
@@ -16,29 +14,28 @@
   } while (0)
 #pragma GCC diagnostic pop
 
-#define SEED
-
 int
 main(void)
 {
-  int i = 0;
+  int i;
 
   // disable buffering
-  setbuf(stderr, 0);
+  setvbuf(stdout, NULL, _IONBF, 0); // Disable buffering for stdout
+  setvbuf(stderr, NULL, _IONBF, 0); // Disable buffering for stdout
 
   // start with a 1KiB buffer, 1024 bits
   uint8_t *buf = calloc(1024, sizeof(uint8_t));
 
   // create the sparse bitmap
-  sparsemap_t *map = sparsemap(buf, sizeof(uint8_t) * 1024);
+  sparsemap_t *map = sparsemap_wrap(buf, sizeof(uint8_t) * 1024);
 
   // Set every other bit (pathologically worst case) to see what happens
   // when the map is full.
   for (i = 0; i < 7744; i++) {
-    if (i % 2)
-      continue;
-    sparsemap_set(map, i, true);
-    assert(sparsemap_is_set(map, i) == true);
+    if (!i % 2) {
+      sparsemap_set(map, i, true);
+      assert(sparsemap_is_set(map, i) == true);
+    }
   }
   // On 1024 KiB of buffer with every other bit set the map holds 7744 bits
   // and then runs out of space.  This next _set() call will fail/abort.
diff --git a/examples/ex_3.c b/examples/ex_3.c
index be76803..095a8d2 100644
--- a/examples/ex_3.c
+++ b/examples/ex_3.c
@@ -1,9 +1,7 @@
 #include <assert.h>
-#include <stdarg.h>
+#include <stdbool.h>
 #include <stdio.h>
 #include <stdlib.h>
-#include <time.h>
-#include <unistd.h>
 
 #include "../include/sparsemap.h"
 #include "../tests/common.h"
@@ -11,7 +9,7 @@
 int
 main(void)
 {
-  int i = 0;
+  int i;
   int array[1024] = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37,
     38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76,
     77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112,
@@ -60,7 +58,7 @@ main(void)
   uint8_t *buf = calloc(1024, sizeof(uint8_t));
 
   // create the sparse bitmap
-  sparsemap_t *map = sparsemap(buf, sizeof(uint8_t) * 1024);
+  sparsemap_t *map = sparsemap_wrap(buf, sizeof(uint8_t) * 1024);
 
   // set all the bits on in a random order
   for (i = 0; i < 1024; i++) {
diff --git a/examples/ex_4.c b/examples/ex_4.c
index 41156e2..b9c0c5e 100644
--- a/examples/ex_4.c
+++ b/examples/ex_4.c
@@ -24,7 +24,7 @@ main(void)
   uint8_t *buf = calloc((size_t)3 * 1024, sizeof(uint8_t));
 
   // create the sparse bitmap
-  sparsemap_t *map = sparsemap(buf, sizeof(uint8_t) * 3 * 1024);
+  sparsemap_t *map = sparsemap_wrap(buf, sizeof(uint8_t) * 3 * 1024);
 
   // create an array of ints
   setup_test_array(array, TEST_ARRAY_SIZE, 1024 * 3);
@@ -60,7 +60,7 @@ main(void)
       assert(sparsemap_is_set(map, array[i]) == true);
     }
     has_span(map, array, TEST_ARRAY_SIZE, (int)len);
-    size_t l = sparsemap_span(map, 0, len);
+    size_t l = sparsemap_span(map, 0, len, true);
     if (l != (size_t)-1) {
       __diag("Found span in map starting at %lu of length %lu\n", l, len);
       __diag("is_span(%lu, %lu) == %s\n", l, len, is_span(array, TEST_ARRAY_SIZE, l, len) ? "yes" : "no");
diff --git a/flake.lock b/flake.lock
index c8aedb6..f0ad25e 100644
--- a/flake.lock
+++ b/flake.lock
@@ -1,43 +1,25 @@
 {
   "nodes": {
-    "flake-utils": {
-      "inputs": {
-        "systems": "systems"
-      },
-      "locked": {
-        "lastModified": 1710146030,
-        "narHash": "sha256-SZ5L6eA7HJ/nmkzGG7/ISclqe6oZdOZTNoesiInkXPQ=",
-        "owner": "numtide",
-        "repo": "flake-utils",
-        "rev": "b1d9ab70662946ef0850d488da1c9019f3a9752a",
-        "type": "github"
-      },
-      "original": {
-        "owner": "numtide",
-        "repo": "flake-utils",
-        "type": "github"
-      }
-    },
     "nixpkgs": {
       "locked": {
-        "lastModified": 1712192574,
-        "narHash": "sha256-LbbVOliJKTF4Zl2b9salumvdMXuQBr2kuKP5+ZwbYq4=",
+        "lastModified": 1701282334,
+        "narHash": "sha256-MxCVrXY6v4QmfTwIysjjaX0XUhqBbxTWWB4HXtDYsdk=",
         "owner": "NixOS",
         "repo": "nixpkgs",
-        "rev": "f480f9d09e4b4cf87ee6151eba068197125714de",
+        "rev": "057f9aecfb71c4437d2b27d3323df7f93c010b7e",
         "type": "github"
       },
       "original": {
         "owner": "NixOS",
-        "ref": "nixpkgs-unstable",
+        "ref": "23.11",
         "repo": "nixpkgs",
         "type": "github"
       }
     },
     "root": {
       "inputs": {
-        "flake-utils": "flake-utils",
-        "nixpkgs": "nixpkgs"
+        "nixpkgs": "nixpkgs",
+        "utils": "utils"
       }
     },
     "systems": {
@@ -54,6 +36,24 @@
         "repo": "default",
         "type": "github"
       }
+    },
+    "utils": {
+      "inputs": {
+        "systems": "systems"
+      },
+      "locked": {
+        "lastModified": 1710146030,
+        "narHash": "sha256-SZ5L6eA7HJ/nmkzGG7/ISclqe6oZdOZTNoesiInkXPQ=",
+        "owner": "numtide",
+        "repo": "flake-utils",
+        "rev": "b1d9ab70662946ef0850d488da1c9019f3a9752a",
+        "type": "github"
+      },
+      "original": {
+        "owner": "numtide",
+        "repo": "flake-utils",
+        "type": "github"
+      }
     }
   },
   "root": "root",
diff --git a/flake.nix b/flake.nix
index e270e61..0c39844 100644
--- a/flake.nix
+++ b/flake.nix
@@ -1,59 +1,55 @@
 {
-  description = "A Concurrent Skip List library for key/value pairs.";
+  description = "A sparse bitmapped index library in C.";
 
   inputs = {
-    nixpkgs.url = "github:NixOS/nixpkgs/nixpkgs-unstable";
-    flake-utils.url = "github:numtide/flake-utils";
+    # nixpkgs.url = "github:NixOS/nixpkgs/nixpkgs-unstable";
+    nixpkgs.url = "github:NixOS/nixpkgs/23.11";
+    utils.url = "github:numtide/flake-utils";
+    utils.inputs.nixpkgs.follows = "nixpkgs";
   };
 
-  outputs =
-    { self
-    , nixpkgs
-    , flake-utils
-    , ...
-    }:
-    flake-utils.lib.eachDefaultSystem (system:
-      let
-        pkgs = import nixpkgs {
-          inherit system;
-          config = { allowUnfree = true; };
-        };
-        supportedSystems = [ "x86_64-linux" ];
-        forAllSystems = nixpkgs.lib.genAttrs supportedSystems;
-        nixpkgsFor = forAllSystems (system: import nixpkgs {
-          inherit system;
-          overlays = [ self.overlay ];
-        });
-      in {
-        pkgs = import nixpkgs {
-          inherit system;
-          devShell = nixpkgs.legacyPackages.${system} {
-            pkgs.mkShell = {
-              nativeBuildInputs = with pkgs.buildPackages; [
-                act
-                autoconf
-                clang
-                ed
-                gcc
-                gdb
-                gettext
-                graphviz-nox
-                libtool
-                m4
-                perl
-                pkg-config
-                python3
-                ripgrep
-                valgrind
-              ];
-              buildInputs = with pkgs; [
-                libbacktrace
-                glibc.out
-                glibc.static
-              ];
-            };
-            DOCKER_BUILDKIT = 1;
+  outputs = { self, nixpkgs, ... }
+    @inputs: inputs.utils.lib.eachSystem [
+      "x86_64-linux" "i686-linux" "aarch64-linux" "x86_64-darwin"
+    ] (system:
+      let pkgs = import nixpkgs {
+            inherit system;
+            overlays = [];
+            config.allowUnfree = true;
           };
-        };
-      });
+      in {
+      devShell = pkgs.mkShell rec {
+        name = "sparsemap";
+        packages = with pkgs; [
+          act
+          autoconf
+          clang
+          ed
+          gcc
+          gdb
+          gettext
+          graphviz-nox
+          libtool
+          m4
+          perl
+          pkg-config
+          python3
+          ripgrep
+          valgrind
+        ];
+
+        buildInputs = with pkgs; [
+          libbacktrace
+          glibc.out
+          glibc.static
+        ];
+
+        shellHook = let
+          icon = "f121";
+        in ''
+        export PS1="$(echo -e '\u${icon}') {\[$(tput sgr0)\]\[\033[38;5;228m\]\w\[$(tput sgr0)\]\[\033[38;5;15m\]} (${name}) \\$ \[$(tput sgr0)\]"
+        '';
+      };
+      DOCKER_BUILDKIT = 1;
+    });
 }
diff --git a/include/sparsemap.h b/include/sparsemap.h
index 548aa3d..689980b 100644
--- a/include/sparsemap.h
+++ b/include/sparsemap.h
@@ -69,14 +69,14 @@
 #ifndef SPARSEMAP_H
 #define SPARSEMAP_H
 
-#include <sys/types.h>
-
-#include <assert.h>
 #include <limits.h>
 #include <stdbool.h>
+#include <stddef.h>
 #include <stdint.h>
-#include <stdio.h>
-#include <string.h>
+
+#if defined(__cplusplus)
+extern "C" {
+#endif
 
 /*
  * The public interface for a sparse bit-mapped index, a "sparse map".
@@ -88,55 +88,119 @@
  */
 
 typedef struct sparsemap sparsemap_t;
+typedef long int sparsemap_idx_t;
+#define SPARSEMAP_IDX_MAX ((1UL << (sizeof(long) * CHAR_BIT - 1)) - 1)
+#define SPARSEMAP_IDX_MIN (-(SPARSEMAP_IDX_MAX)-1)
+#define SPARSEMAP_NOT_FOUND(_x) ((_x) == SPARSEMAP_IDX_MAX || (_x) == SPARSEMAP_IDX_MIN)
 typedef uint32_t sm_idx_t;
 typedef uint64_t sm_bitvec_t;
 
-/* Allocate on a sparsemap_t on the heap and initialize it. */
-sparsemap_t *sparsemap(uint8_t *data, size_t size);
+/**
+ * Create a new, empty sparsemap_t with a buffer of |size|.
+ * Default when set to 0 is 1024.
+ */
+sparsemap_t *sparsemap(size_t size);
 
-/* Initialize sparsemap_t with data. */
+/**
+ * Allocate on a sparsemap_t on the heap to wrap the provided fixed-size
+ * buffer (heap or stack allocated).
+ */
+sparsemap_t *sparsemap_wrap(uint8_t *data, size_t size);
+
+/**
+ * Initialize a (possibly stack allocated) sparsemap_t with data (potentially
+ * also on the stack).
+ */
 void sparsemap_init(sparsemap_t *map, uint8_t *data, size_t size);
 
-/* Clears the whole buffer. */
-void sparsemap_clear(sparsemap_t *map);
-
-/* Opens an existing sparsemap at the specified buffer. */
+/**
+ * Opens an existing sparsemap contained within the specified buffer.
+ */
 void sparsemap_open(sparsemap_t *, uint8_t *data, size_t data_size);
 
-/* Resizes the data range. */
-void sparsemap_set_data_size(sparsemap_t *map, size_t data_size);
+/**
+ * Resets values and empties the buffer making it ready to accept new data.
+ */
+void sparsemap_clear(sparsemap_t *map);
 
-/* Calculate remaining capacity, full when 0. */
+/**
+ * Resizes the data range within the limits of the provided buffer, the map may
+ * move to a new address returned iff the map was created with the sparsemap() API.
+ * Take care to use the new reference (think: realloc()).  NOTE: If the returned
+ * value equals NULL then the map was not resized.
+ */
+sparsemap_t *sparsemap_set_data_size(sparsemap_t *map, size_t data_size);
+
+/**
+ * Calculate remaining capacity, approaches 0 when full.
+ */
 double sparsemap_capacity_remaining(sparsemap_t *map);
 
-/* Returns the size of the underlying byte array. */
+/**
+ * Returns the capacity of the underlying byte array.
+ */
 size_t sparsemap_get_capacity(sparsemap_t *map);
 
-/* Returns the value of a bit at index |idx|. */
-bool sparsemap_is_set(sparsemap_t *map, size_t idx);
+/**
+ * Returns the value of a bit at index |idx|, either on/true/1 or off/false/0.
+ * When |idx| is negative it is an error.
+ */
+bool sparsemap_is_set(sparsemap_t *map, sparsemap_idx_t idx);
 
-/* Sets the bit at index |idx| to true or false, depending on |value|. */
-void sparsemap_set(sparsemap_t *map, size_t idx, bool value);
+/**
+ * Sets the bit at index |idx| to true or false, depending on |value|.
+ * When |idx| is negative is it an error.  Returns the |idx| supplied or
+ * SPARSEMAP_IDX_MAX on error with |errno| set to ENOSP when the map is full.
+ */
+sparsemap_idx_t sparsemap_set(sparsemap_t *map, sparsemap_idx_t idx, bool value);
 
-/* Returns the offset of the very first bit. */
-sm_idx_t sparsemap_get_start_offset(sparsemap_t *map);
+/**
+ * Returns the offset of the very first/last bit in the map.
+ */
+sm_idx_t sparsemap_get_starting_offset(sparsemap_t *map);
 
-/* Returns the used size in the data buffer. */
+/**
+ * Returns the used size in the data buffer in bytes.
+ */
 size_t sparsemap_get_size(sparsemap_t *map);
 
-/* Decompresses the whole bitmap; calls scanner for all bits. */
-void sparsemap_scan(sparsemap_t *map, void (*scanner)(sm_idx_t[], size_t), size_t skip);
+/**
+ * Decompresses the whole bitmap; calls scanner for all bits with a set of
+ * |n| vectors |vec| each a sm_bitmap_t which can be masked and read using
+ * bit operators to read the values for each position in the bitmap index.
+ * Setting |skip| will start the scan after "skip" bits.
+ */
+void sparsemap_scan(sparsemap_t *map, void (*scanner)(sm_idx_t vec[], size_t n), size_t skip);
 
-/* Appends all chunk maps from |map| starting at |sstart| to |other|, then
-   reduces the chunk map-count appropriately. */
-void sparsemap_split(sparsemap_t *map, size_t sstart, sparsemap_t *other);
+/**
+ * Appends all chunk maps from |map| starting at |offset| to |other|, then
+ * reduces the chunk map-count appropriately.
+ */
+void sparsemap_split(sparsemap_t *map, sparsemap_idx_t offset, sparsemap_t *other);
 
-/* Returns the index of the n'th set bit; uses a 0-based index. */
-size_t sparsemap_select(sparsemap_t *map, size_t n);
+/**
+ * Finds the offset of the n'th bit either set (|value| is true) or unset
+ * (|value| is false) from the start (positive |n|), or end (negative |n|),
+ * of the bitmap and returns that (uses a 0-based index).  Returns -inf or +inf
+ * if not found (where "inf" is SPARSEMAP_IDX_MAX and "-inf" is SPARSEMAP_IDX_MIN).
+ */
+sparsemap_idx_t sparsemap_select(sparsemap_t *map, sparsemap_idx_t n, bool value);
 
-/* Counts the set bits in the range [offset, idx]. */
-size_t sparsemap_rank(sparsemap_t *map, size_t offset, size_t idx);
+/**
+ * Counts the set (|value| is true) or unset (|value| is false) bits starting
+ * at |x| bits (0-based) in the range [x, y] (inclusive on either end).
+ */
+size_t sparsemap_rank(sparsemap_t *map, size_t x, size_t y, bool value);
 
-size_t sparsemap_span(sparsemap_t *map, size_t loc, size_t len);
+/**
+ * Finds the first span (i.e. a contiguous set of bits), in the bitmap that
+ * are set (|value| is true) or unset (|value| is false) and returns the
+ * starting offset for the span (0-based).
+ */
+size_t sparsemap_span(sparsemap_t *map, sparsemap_idx_t idx, size_t len, bool value);
 
+#if defined(__cplusplus)
+}
 #endif
+
+#endif /* !defined(SPARSEMAP_H) */
diff --git a/src/sparsemap.c b/src/sparsemap.c
index eb7fd43..c5d88c5 100644
--- a/src/sparsemap.c
+++ b/src/sparsemap.c
@@ -20,12 +20,19 @@
  * SOFTWARE.
  */
 
+#include <sys/types.h>
+
 #include <assert.h>
+#include <errno.h>
 #include <popcount.h>
 #include <sparsemap.h>
+#include <stdarg.h>
 #include <stdbool.h>
+#include <stddef.h>
 #include <stdint.h>
+#include <stdio.h>
 #include <stdlib.h>
+#include <string.h>
 
 #ifdef SPARSEMAP_DIAGNOSTIC
 #pragma GCC diagnostic push
@@ -42,19 +49,24 @@ void __attribute__((format(printf, 4, 5))) __sm_diag_(const char *file, int line
   vfprintf(stderr, format, args);
   va_end(args);
 }
-#else
-#define __sm_diag(file, line, func, format, ...) ((void)0)
-#endif
 
-#ifndef SPARSEMAP_ASSERT
-#define SPARSEMAP_ASSERT
 #define __sm_assert(expr) \
   if (!(expr))            \
   fprintf(stderr, "%s:%d:%s(): assertion failed! %s", __FILE__, __LINE__, __func__, #expr)
+
+#define __sm_when_diag(expr) \
+  if (1)                     \
+  expr
 #else
+#define __sm_diag(file, line, func, format, ...) ((void)0)
 #define __sm_assert(expr) ((void)0)
+#define __sm_when_diag(expr) \
+  if (0)                     \
+  expr
 #endif
 
+#define IS_8_BYTE_ALIGNED(addr) (((uintptr_t)(addr)&0x7) == 0)
+
 enum __SM_CHUNK_INFO {
   /* metadata overhead: 4 bytes for __sm_chunk_t count */
   SM_SIZEOF_OVERHEAD = sizeof(uint32_t),
@@ -83,7 +95,7 @@ enum __SM_CHUNK_INFO {
   /* sm_bitvec_t is not used (2#01) */
   SM_PAYLOAD_NONE = 1,
 
-  /* a mask for checking flags (2 bits) */
+  /* a mask for checking flags (2 bits, 2#11) */
   SM_FLAG_MASK = 3,
 
   /* return code for set(): ok, no further action required */
@@ -96,16 +108,16 @@ enum __SM_CHUNK_INFO {
   SM_NEEDS_TO_SHRINK = 2
 };
 
-#define SM_CHUNK_GET_FLAGS(from, at) (((from)) & ((sm_bitvec_t)SM_FLAG_MASK << ((at) * 2))) >> ((at) * 2)
+#define SM_CHUNK_GET_FLAGS(from, at) (((from)) & ((sm_bitvec_t)SM_FLAG_MASK << ((at)*2))) >> ((at)*2)
 
 typedef struct {
   sm_bitvec_t *m_data;
 } __sm_chunk_t;
 
-struct sparsemap {
-  uint8_t *m_data;    /* The serialized bitmap data */
+struct __attribute__((aligned(8))) sparsemap {
   size_t m_capacity;  /* The total size of m_data */
   size_t m_data_used; /* The used size of m_data */
+  uint8_t *m_data;    /* The serialized bitmap data */
 };
 
 /**
@@ -135,7 +147,7 @@ __sm_chunk_map_calc_vector_size(uint8_t b)
     0,  0,  1,  0,  0,  0,  1,  0,  1,  1,  2,  1,  0,  0,  1,  0
   };
   // clang-format on
-  return ((size_t)lookup[b]);
+  return (size_t)lookup[b];
 }
 
 /**
@@ -161,7 +173,7 @@ __sm_chunk_map_get_position(__sm_chunk_t *map, size_t bv)
     }
   }
 
-  return (position);
+  return position;
 }
 
 /**
@@ -193,7 +205,7 @@ __sm_chunk_map_get_capacity(__sm_chunk_t *map)
       }
     }
   }
-  return (capacity);
+  return capacity;
 }
 
 /**
@@ -232,7 +244,7 @@ __sm_chunk_map_is_empty(__sm_chunk_t *map)
 {
   /* The __sm_chunk_t is empty if all flags (in m_data[0]) are zero. */
   if (map->m_data[0] == 0) {
-    return (true);
+    return true;
   }
 
   /* It's also empty if all flags are Zero or None. */
@@ -242,12 +254,12 @@ __sm_chunk_map_is_empty(__sm_chunk_t *map)
       for (int j = 0; j < SM_FLAGS_PER_INDEX_BYTE; j++) {
         size_t flags = SM_CHUNK_GET_FLAGS(*p, j);
         if (flags != SM_PAYLOAD_NONE && flags != SM_PAYLOAD_ZEROS) {
-          return (false);
+          return false;
         }
       }
     }
   }
-  return (true);
+  return true;
 }
 
 /**
@@ -264,7 +276,7 @@ __sm_chunk_map_get_size(__sm_chunk_t *map)
     size += sizeof(sm_bitvec_t) * __sm_chunk_map_calc_vector_size(*p);
   }
 
-  return (size);
+  return size;
 }
 
 /**
@@ -282,9 +294,9 @@ __sm_chunk_map_is_set(__sm_chunk_t *map, size_t idx)
   switch (flags) {
   case SM_PAYLOAD_ZEROS:
   case SM_PAYLOAD_NONE:
-    return (false);
+    return false;
   case SM_PAYLOAD_ONES:
-    return (true);
+    return true;
   default:
     __sm_assert(flags == SM_PAYLOAD_MIXED);
     /* FALLTHROUGH */
@@ -293,7 +305,7 @@ __sm_chunk_map_is_set(__sm_chunk_t *map, size_t idx)
   /* get the sm_bitvec_t at |bv| */
   sm_bitvec_t w = map->m_data[1 + __sm_chunk_map_get_position(map, bv)];
   /* and finally check the bit in that sm_bitvec_t */
-  return ((w & ((sm_bitvec_t)1 << (idx % SM_BITS_PER_VECTOR))) > 0);
+  return (w & ((sm_bitvec_t)1 << (idx % SM_BITS_PER_VECTOR))) > 0;
 }
 
 /**
@@ -388,7 +400,7 @@ __sm_chunk_map_set(__sm_chunk_t *map, size_t idx, bool value, size_t *pos, sm_bi
  * value of |n|.
  */
 static size_t
-__sm_chunk_map_select(__sm_chunk_t *map, size_t n, ssize_t *pnew_n)
+__sm_chunk_map_select(__sm_chunk_t *map, size_t n, ssize_t *pnew_n, bool value)
 {
   size_t ret = 0;
   register uint8_t *p;
@@ -406,50 +418,86 @@ __sm_chunk_map_select(__sm_chunk_t *map, size_t n, ssize_t *pnew_n)
         continue;
       }
       if (flags == SM_PAYLOAD_ZEROS) {
-        ret += SM_BITS_PER_VECTOR;
-        continue;
+        if (value) {
+          ret += SM_BITS_PER_VECTOR;
+          continue;
+        } else {
+          if (n > SM_BITS_PER_VECTOR) {
+            n -= SM_BITS_PER_VECTOR;
+            ret += SM_BITS_PER_VECTOR;
+            continue;
+          }
+          *pnew_n = -1;
+          return ret + n;
+        }
       }
       if (flags == SM_PAYLOAD_ONES) {
-        if (n > SM_BITS_PER_VECTOR) {
-          n -= SM_BITS_PER_VECTOR;
+        if (value) {
+          if (n > SM_BITS_PER_VECTOR) {
+            n -= SM_BITS_PER_VECTOR;
+            ret += SM_BITS_PER_VECTOR;
+            continue;
+          }
+          *pnew_n = -1;
+          return ret + n;
+        } else {
           ret += SM_BITS_PER_VECTOR;
           continue;
         }
-
-        *pnew_n = -1;
-        return (ret + n);
       }
       if (flags == SM_PAYLOAD_MIXED) {
         sm_bitvec_t w = map->m_data[1 + __sm_chunk_map_get_position(map, i * SM_FLAGS_PER_INDEX_BYTE + j)];
         for (int k = 0; k < SM_BITS_PER_VECTOR; k++) {
-          if (w & ((sm_bitvec_t)1 << k)) {
-            if (n == 0) {
-              *pnew_n = -1;
-              return (ret);
+          if (value) {
+            if (w & ((sm_bitvec_t)1 << k)) {
+              if (n == 0) {
+                *pnew_n = -1;
+                return ret;
+              }
+              n--;
             }
-            n--;
+            ret++;
+          } else {
+            if (!(w & ((sm_bitvec_t)1 << k))) {
+              if (n == 0) {
+                *pnew_n = -1;
+                return ret;
+              }
+              n--;
+            }
+            ret++;
           }
-          ret++;
         }
       }
     }
   }
-
   *pnew_n = (ssize_t)n;
-  return (ret);
+  return ret;
 }
 
+extern void print_bits(char *name, uint64_t value); // GSB
+
 /**
  * Counts the set bits in the range [0, 'idx'] inclusive ignoring the first
- * '*offset' bits.  Modifies '*offset' decreasing it by the number of bits
- * ignored during the search.  The ranking (counting) will start after the
- * '*offset' has been reached 0.
+ * '*offset' bits in this chunk.  Modifies '*offset' decreasing it by the number
+ * of bits ignored during the search.  The ranking (counting) will start after
+ * the '*offset' has been reached 0.
  */
 static size_t
-__sm_chunk_map_rank(__sm_chunk_t *map, size_t *offset, size_t idx)
+__sm_chunk_map_rank(__sm_chunk_t *map, size_t *offset, size_t idx, size_t *pos, sm_bitvec_t *vec, bool value)
 {
   size_t ret = 0;
 
+  *pos = 0;
+
+  /* A chunk can only hold at most SM_CHUNK_MAX_CAPACITY bits, so if the
+     offset is larger than that, we're basically done. */
+  if (*offset > SM_CHUNK_MAX_CAPACITY) {
+    *pos = SM_CHUNK_MAX_CAPACITY;
+    *offset -= SM_CHUNK_MAX_CAPACITY;
+    return 0;
+  }
+
   register uint8_t *p = (uint8_t *)map->m_data;
   for (size_t i = 0; i < sizeof(sm_bitvec_t); i++, p++) {
     for (int j = 0; j < SM_FLAGS_PER_INDEX_BYTE; j++) {
@@ -458,49 +506,100 @@ __sm_chunk_map_rank(__sm_chunk_t *map, size_t *offset, size_t idx)
         continue;
       }
       if (flags == SM_PAYLOAD_ZEROS) {
-        if (idx > SM_BITS_PER_VECTOR) {
+        *vec = 0;
+        if (idx >= SM_BITS_PER_VECTOR) {
+          *pos += SM_BITS_PER_VECTOR;
+          idx -= SM_BITS_PER_VECTOR;
           if (*offset > SM_BITS_PER_VECTOR) {
             *offset = *offset - SM_BITS_PER_VECTOR;
           } else {
-            idx -= SM_BITS_PER_VECTOR - *offset;
-            *offset = 0;
-          }
-        } else {
-          return (ret);
-        }
-      } else if (flags == SM_PAYLOAD_ONES) {
-        if (idx > SM_BITS_PER_VECTOR) {
-          if (*offset > SM_BITS_PER_VECTOR) {
-            *offset = *offset - SM_BITS_PER_VECTOR;
-          } else {
-            idx -= SM_BITS_PER_VECTOR - *offset;
-            if (*offset == 0) {
-              ret += SM_BITS_PER_VECTOR;
+            if (value == false) {
+              ret += SM_BITS_PER_VECTOR - *offset;
             }
             *offset = 0;
           }
         } else {
-          return (ret + idx);
+          *pos += idx + 1;
+          if (value == false) {
+            if (*offset > idx) {
+              *offset = *offset - idx;
+            } else {
+              ret += idx + 1 - *offset;
+              *offset = 0;
+              return ret;
+            }
+          } else {
+            return ret;
+          }
+        }
+      } else if (flags == SM_PAYLOAD_ONES) {
+        *vec = UINT64_MAX;
+        if (idx >= SM_BITS_PER_VECTOR) {
+          *pos += SM_BITS_PER_VECTOR;
+          idx -= SM_BITS_PER_VECTOR;
+          if (*offset > SM_BITS_PER_VECTOR) {
+            *offset = *offset - SM_BITS_PER_VECTOR;
+          } else {
+            if (value == true) {
+              ret += SM_BITS_PER_VECTOR - *offset;
+            }
+            *offset = 0;
+          }
+        } else {
+          *pos += idx + 1;
+          if (value == true) {
+            if (*offset > idx) {
+              *offset = *offset - idx;
+            } else {
+              ret += idx + 1 - *offset;
+              *offset = 0;
+              return ret;
+            }
+          } else {
+            return ret;
+          }
         }
       } else if (flags == SM_PAYLOAD_MIXED) {
         sm_bitvec_t w = map->m_data[1 + __sm_chunk_map_get_position(map, i * SM_FLAGS_PER_INDEX_BYTE + j)];
-        if (idx > SM_BITS_PER_VECTOR) {
-          uint64_t mask_offset = ~(UINT64_MAX >> (SM_BITS_PER_VECTOR - *offset));
+        if (idx >= SM_BITS_PER_VECTOR) {
+          *pos += SM_BITS_PER_VECTOR;
           idx -= SM_BITS_PER_VECTOR;
-          ret += popcountll(w & mask_offset);
+          uint64_t mask = *offset == 0 ? UINT64_MAX : ~(UINT64_MAX >> (SM_BITS_PER_VECTOR - (*offset >= 64 ? 64 : *offset)));
+          sm_bitvec_t mw;
+          if (value == true) {
+            mw = w & mask;
+          } else {
+            mw = ~w & mask;
+          }
+          size_t pc = popcountll(mw);
+          ret += pc;
           *offset = (*offset > SM_BITS_PER_VECTOR) ? *offset - SM_BITS_PER_VECTOR : 0;
         } else {
-          /* Create a mask for the range between offset and idx inclusive [*offset, idx]. */
-          uint64_t offset_mask = (((uint64_t)1 << *offset) - 1);
-          uint64_t idx_mask = idx >= 63 ? UINT64_MAX : ((uint64_t)1 << (idx + 1)) - 1;
-          ret += popcountll(w & (idx_mask - offset_mask));
-          *offset = *offset > idx ? *offset - idx : 0;
-          return (ret);
+          *pos += idx + 1;
+          sm_bitvec_t mw;
+          uint64_t mask;
+          uint64_t idx_mask = (idx == 63) ? UINT64_MAX : ((uint64_t)1 << (idx + 1)) - 1;
+          uint64_t offset_mask = *offset == 0 ? UINT64_MAX : ~(UINT64_MAX >> (SM_BITS_PER_VECTOR - (*offset >= 64 ? 64 : *offset)));
+          /* To count the set bits we need to mask off the portion of the vector that we need
+             to count then call popcount().  So, let's create a mask for the range between
+             offset and idx inclusive [*offset, idx]. */
+          mask = idx_mask & offset_mask;
+          if (value == true) {
+            mw = w & mask;
+          } else {
+            mw = ~w & mask;
+          }
+          int pc = popcountll(mw);
+          ret += pc;
+          *offset = *offset > idx ? *offset - idx + 1 : 0;
+          *vec = mw;
+          (*vec) <<= *offset;
+          return ret;
         }
       }
     }
   }
-  return (ret);
+  return ret;
 }
 
 /**
@@ -571,7 +670,7 @@ __sm_chunk_map_scan(__sm_chunk_t *map, sm_idx_t start, void (*scanner)(sm_idx_t[
       }
     }
   }
-  return (ret);
+  return ret;
 }
 
 /*
@@ -584,7 +683,7 @@ __sm_chunk_map_scan(__sm_chunk_t *map, sm_idx_t start, void (*scanner)(sm_idx_t[
 static size_t
 __sm_get_chunk_map_count(sparsemap_t *map)
 {
-  return (*(uint32_t *)&map->m_data[0]);
+  return *(uint32_t *)&map->m_data[0];
 }
 
 /**
@@ -593,7 +692,7 @@ __sm_get_chunk_map_count(sparsemap_t *map)
 static inline uint8_t *
 __sm_get_chunk_map_data(sparsemap_t *map, size_t offset)
 {
-  return (uint8_t *)(&map->m_data[SM_SIZEOF_OVERHEAD + offset]);
+  return &map->m_data[SM_SIZEOF_OVERHEAD + offset];
 }
 
 /**
@@ -611,7 +710,7 @@ __sm_get_chunk_map_end(sparsemap_t *map)
     __sm_chunk_map_init(&chunk, p);
     p += __sm_chunk_map_get_size(&chunk);
   }
-  return (p);
+  return p;
 }
 
 /**
@@ -630,7 +729,7 @@ __sm_get_size_impl(sparsemap_t *map)
     __sm_chunk_map_init(&chunk, p);
     p += __sm_chunk_map_get_size(&chunk);
   }
-  return (SM_SIZEOF_OVERHEAD + p - start);
+  return SM_SIZEOF_OVERHEAD + p - start;
 }
 
 /**
@@ -640,37 +739,55 @@ static sm_idx_t
 __sm_get_aligned_offset(size_t idx)
 {
   const size_t capacity = SM_BITS_PER_VECTOR;
-  return ((idx / capacity) * capacity);
+  return (idx / capacity) * capacity;
 }
 
 /**
- * Returns the byte offset of a __sm_chunk_t in m_data
+ * Returns the byte offset of a __sm_chunk_t in m_data.
  */
 static ssize_t
-__sm_get_chunk_map_offset(sparsemap_t *map, size_t idx)
+__sm_get_chunk_map_offset(sparsemap_t *map, sparsemap_idx_t idx)
 {
-  size_t count;
+  int count;
 
   count = __sm_get_chunk_map_count(map);
   if (count == 0) {
-    return (-1);
+    return -1;
   }
 
-  uint8_t *start = __sm_get_chunk_map_data(map, 0);
-  uint8_t *p = start;
+  if (idx > 0 || idx == 0) {
+    uint8_t *start = __sm_get_chunk_map_data(map, 0);
+    uint8_t *p = start;
 
-  for (size_t i = 0; i < count - 1; i++) {
-    sm_idx_t start = *(sm_idx_t *)p;
-    __sm_assert(start == __sm_get_aligned_offset(start));
-    __sm_chunk_t chunk;
-    __sm_chunk_map_init(&chunk, p + sizeof(sm_idx_t));
-    if (start >= idx || idx < start + __sm_chunk_map_get_capacity(&chunk)) {
-      break;
+    for (sparsemap_idx_t i = 0; i < count - 1; i++) {
+      sm_idx_t s = *(sm_idx_t *)p;
+      __sm_assert(s == __sm_get_aligned_offset(s));
+      __sm_chunk_t chunk;
+      __sm_chunk_map_init(&chunk, p + sizeof(sm_idx_t));
+      if (s >= idx || (unsigned long)idx < s + __sm_chunk_map_get_capacity(&chunk)) {
+        break;
+      }
+      p += sizeof(sm_idx_t) + __sm_chunk_map_get_size(&chunk);
     }
-    p += sizeof(sm_idx_t) + __sm_chunk_map_get_size(&chunk);
-  }
 
-  return ((ssize_t)(p - start));
+    return (ssize_t)(p - start);
+  } else {
+    uint8_t *end = __sm_get_chunk_map_data(map, count - 1);
+    uint8_t *p = end;
+
+    for (sparsemap_idx_t i = count - 1; i >= 0; i--) {
+      sm_idx_t e = *(sm_idx_t *)p;
+      __sm_assert(e == __sm_get_aligned_offset(e));
+      __sm_chunk_t chunk;
+      __sm_chunk_map_init(&chunk, p + sizeof(sm_idx_t));
+      if (e >= idx || (unsigned long)idx < e + __sm_chunk_map_get_capacity(&chunk)) {
+        break;
+      }
+      p += sizeof(sm_idx_t) + __sm_chunk_map_get_size(&chunk);
+    }
+
+    return (ssize_t)(p - end);
+  }
 }
 
 /**
@@ -680,7 +797,7 @@ static sm_idx_t
 __sm_get_fully_aligned_offset(size_t idx)
 {
   const size_t capacity = SM_CHUNK_MAX_CAPACITY;
-  return ((idx / capacity) * capacity);
+  return (idx / capacity) * capacity;
 }
 
 /**
@@ -705,19 +822,13 @@ __sm_append_data(sparsemap_t *map, uint8_t *buffer, size_t buffer_size)
 /**
  * Inserts data somewhere in the middle of m_data.
  */
-static int
+void
 __sm_insert_data(sparsemap_t *map, size_t offset, uint8_t *buffer, size_t buffer_size)
 {
-  if (map->m_data_used + buffer_size > map->m_capacity) {
-    __sm_assert(!"buffer overflow");
-    abort();
-  }
-
   uint8_t *p = __sm_get_chunk_map_data(map, offset);
   memmove(p + buffer_size, p, map->m_data_used - offset);
   memcpy(p, buffer, buffer_size);
   map->m_data_used += buffer_size;
-  return 0;
 }
 
 /**
@@ -738,16 +849,41 @@ __sm_remove_data(sparsemap_t *map, size_t offset, size_t gap_size)
 void
 sparsemap_clear(sparsemap_t *map)
 {
+  if (map == NULL) {
+    return;
+  }
   memset(map->m_data, 0, map->m_capacity);
   map->m_data_used = SM_SIZEOF_OVERHEAD;
   __sm_set_chunk_map_count(map, 0);
 }
 
-/**
- * Allocate on a sparsemap_t on the heap and initialize it.
- */
 sparsemap_t *
-sparsemap(uint8_t *data, size_t size)
+sparsemap(size_t size)
+{
+  if (size == 0) {
+    size = 1024;
+  }
+
+  size_t data_size = (size * sizeof(uint8_t));
+
+  /* Ensure that m_data is 8-byte aligned. */
+  size_t total_size = sizeof(sparsemap_t) + data_size;
+  size_t padding = total_size % 8 == 0 ? 0 : 8 - (total_size % 8);
+  total_size += padding;
+
+  sparsemap_t *map = (sparsemap_t *)calloc(1, total_size);
+  if (map) {
+    uint8_t *data = (uint8_t *)(((uintptr_t)map + sizeof(sparsemap_t)) & ~(uintptr_t)7);
+    sparsemap_init(map, data, size);
+    __sm_when_diag({ __sm_assert(IS_8_BYTE_ALIGNED(map->m_data)); });
+  }
+  sparsemap_clear(map);
+
+  return map;
+}
+
+sparsemap_t *
+sparsemap_wrap(uint8_t *data, size_t size)
 {
   sparsemap_t *map = (sparsemap_t *)calloc(1, sizeof(sparsemap_t));
   if (map) {
@@ -756,45 +892,53 @@ sparsemap(uint8_t *data, size_t size)
   return map;
 }
 
-/**
- * Initialize sparsemap_t with data.
- */
 void
 sparsemap_init(sparsemap_t *map, uint8_t *data, size_t size)
 {
   map->m_data = data;
   map->m_data_used = 0;
-  map->m_capacity = size == 0 ? UINT64_MAX : size;
+  map->m_capacity = size;
   sparsemap_clear(map);
 }
 
-/**
- * Opens an existing sparsemap at the specified buffer.
- */
 void
 sparsemap_open(sparsemap_t *map, uint8_t *data, size_t data_size)
 {
   map->m_data = data;
-  map->m_data_used = 0;
+  map->m_data_used = map->m_data_used > 0 ? map->m_data_used : 0;
   map->m_capacity = data_size;
 }
 
-/**
- * Resizes the data range.
- *
+/*
  * TODO/NOTE: This is a dangerous operation because we cannot verify that
  *       data_size is not exceeding the size of the underlying buffer.
  */
-void
-sparsemap_set_data_size(sparsemap_t *map, size_t data_size)
+sparsemap_t *
+sparsemap_set_data_size(sparsemap_t *map, size_t size)
 {
-  map->m_capacity = data_size;
+  if ((uintptr_t)map->m_data == (uintptr_t)map + sizeof(sparsemap_t) && size > map->m_capacity) {
+    /* This sparsemap was allocated by the sparsemap() API, we can resize it. */
+    size_t data_size = (size * sizeof(uint8_t));
+
+    /* Ensure that m_data is 8-byte aligned. */
+    size_t total_size = sizeof(sparsemap_t) + data_size;
+    size_t padding = total_size % 8 == 0 ? 0 : 8 - (total_size % 8);
+    total_size += padding;
+
+    sparsemap_t *m = (sparsemap_t *)realloc(map, total_size);
+    if (!m) {
+      return NULL;
+    }
+    memset(((uint8_t *)m) + sizeof(sparsemap_t) + (m->m_capacity * sizeof(uint8_t)), 0, size - m->m_capacity + padding);
+    m->m_capacity = data_size;
+    m->m_data = (uint8_t *)(((uintptr_t)m + sizeof(sparsemap_t)) & ~(uintptr_t)7);
+    __sm_when_diag({ __sm_assert(IS_8_BYTE_ALIGNED(m->m_data)); }) return m;
+  } else {
+    map->m_capacity = size;
+    return map;
+  }
 }
 
-/**
- * Calculates the remaining capacity as an integer that approaches 0 to
- * indicate full.
- */
 double
 sparsemap_capacity_remaining(sparsemap_t *map)
 {
@@ -807,29 +951,27 @@ sparsemap_capacity_remaining(sparsemap_t *map)
   return 100 - (((double)map->m_data_used / (double)map->m_capacity) * 100);
 }
 
-/**
- * Returns the size of the underlying byte array.
- */
 size_t
 sparsemap_get_capacity(sparsemap_t *map)
 {
-  return (map->m_capacity);
+  return map->m_capacity;
 }
 
-/**
- * Returns the value of a bit at index |idx|.
- */
 bool
-sparsemap_is_set(sparsemap_t *map, size_t idx)
+sparsemap_is_set(sparsemap_t *map, sparsemap_idx_t idx)
 {
   __sm_assert(sparsemap_get_size(map) >= SM_SIZEOF_OVERHEAD);
 
+  if (idx < 0) {
+    return false;
+  }
+
   /* Get the __sm_chunk_t which manages this index */
   ssize_t offset = __sm_get_chunk_map_offset(map, idx);
 
   /* No __sm_chunk_t's available -> the bit is not set */
   if (offset == -1) {
-    return (false);
+    return false;
   }
 
   /* Otherwise load the __sm_chunk_t */
@@ -839,20 +981,17 @@ sparsemap_is_set(sparsemap_t *map, size_t idx)
   __sm_chunk_map_init(&chunk, p + sizeof(sm_idx_t));
 
   /* Determine if the bit is out of bounds of the __sm_chunk_t; if yes then
-    the bit is not set. */
-  if (idx < start || idx - start >= __sm_chunk_map_get_capacity(&chunk)) {
-    return (false);
+  the bit is not set. */
+  if (idx < start || (unsigned long)idx - start >= __sm_chunk_map_get_capacity(&chunk)) {
+    return false;
   }
 
   /* Otherwise ask the __sm_chunk_t whether the bit is set. */
-  return (__sm_chunk_map_is_set(&chunk, idx - start));
+  return __sm_chunk_map_is_set(&chunk, idx - start);
 }
 
-/**
- * Sets the bit at index |idx| to true or false, depending on |value|.
- */
-void
-sparsemap_set(sparsemap_t *map, size_t idx, bool value)
+sparsemap_idx_t
+sparsemap_set(sparsemap_t *map, sparsemap_idx_t idx, bool value)
 {
   __sm_assert(sparsemap_get_size(map) >= SM_SIZEOF_OVERHEAD);
 
@@ -860,11 +999,16 @@ sparsemap_set(sparsemap_t *map, size_t idx, bool value)
   ssize_t offset = __sm_get_chunk_map_offset(map, idx);
   bool dont_grow = false;
 
+  if (map->m_data_used + sizeof(sm_idx_t) + sizeof(sm_bitvec_t) * 2 > map->m_capacity) {
+    errno = ENOSPC;
+    return SPARSEMAP_IDX_MAX;
+  }
+
   /* If there is no __sm_chunk_t and the bit is set to zero then return
      immediately; otherwise create an initial __sm_chunk_t. */
   if (offset == -1) {
     if (value == false) {
-      return;
+      return idx;
     }
 
     uint8_t buf[sizeof(sm_idx_t) + sizeof(sm_bitvec_t) * 2] = { 0 };
@@ -875,8 +1019,8 @@ sparsemap_set(sparsemap_t *map, size_t idx, bool value)
 
     __sm_set_chunk_map_count(map, 1);
 
-    /* We already inserted an additional sm_bitvec_t; later on there
-       is no need to grow the vector even further. */
+    /* We already inserted an additional sm_bitvec_t; given that has happened
+       there is no need to grow the vector even further. */
     dont_grow = true;
     offset = 0;
   }
@@ -890,7 +1034,7 @@ sparsemap_set(sparsemap_t *map, size_t idx, bool value)
   if (idx < start) {
     if (value == false) {
       /* nothing to do */
-      return;
+      return idx;
     }
 
     uint8_t buf[sizeof(sm_idx_t) + sizeof(sm_bitvec_t) * 2] = { 0 };
@@ -917,10 +1061,10 @@ sparsemap_set(sparsemap_t *map, size_t idx, bool value)
   else {
     __sm_chunk_t chunk;
     __sm_chunk_map_init(&chunk, p + sizeof(sm_idx_t));
-    if (idx - start >= __sm_chunk_map_get_capacity(&chunk)) {
+    if (idx - start >= (sparsemap_idx_t)__sm_chunk_map_get_capacity(&chunk)) {
       if (value == false) {
         /* nothing to do */
-        return;
+        return idx;
       }
 
       size_t size = __sm_chunk_map_get_size(&chunk);
@@ -931,7 +1075,7 @@ sparsemap_set(sparsemap_t *map, size_t idx, bool value)
       __sm_insert_data(map, offset, &buf[0], sizeof(buf));
 
       start += __sm_chunk_map_get_capacity(&chunk);
-      if ((size_t)start + SM_CHUNK_MAX_CAPACITY < idx) {
+      if ((sparsemap_idx_t)start + SM_CHUNK_MAX_CAPACITY < idx) {
         start = __sm_get_fully_aligned_offset(idx);
       }
       *(sm_idx_t *)p = start;
@@ -960,8 +1104,10 @@ sparsemap_set(sparsemap_t *map, size_t idx, bool value)
       offset += (ssize_t)(sizeof(sm_idx_t) + position * sizeof(sm_bitvec_t));
       __sm_insert_data(map, offset, (uint8_t *)&fill, sizeof(sm_bitvec_t));
     }
-    code = __sm_chunk_map_set(&chunk, idx - start, value, &position, &fill, true);
-    __sm_assert(code == SM_OK);
+    __sm_when_diag({
+      code = __sm_chunk_map_set(&chunk, idx - start, value, &position, &fill, true);
+      __sm_assert(code == SM_OK);
+    });
     break;
   case SM_NEEDS_TO_SHRINK:
     /* If the __sm_chunk_t is empty then remove it. */
@@ -982,18 +1128,18 @@ sparsemap_set(sparsemap_t *map, size_t idx, bool value)
     break;
   }
   __sm_assert(sparsemap_get_size(map) >= SM_SIZEOF_OVERHEAD);
+  return idx;
 }
 
-/**
- * Returns the offset of the very first bit.
- */
 sm_idx_t
-sparsemap_get_start_offset(sparsemap_t *map)
+sparsemap_get_starting_offset(sparsemap_t *map)
 {
-  if (__sm_get_chunk_map_count(map) == 0) {
-    return (0);
+  size_t count = __sm_get_chunk_map_count(map);
+  if (count == 0) {
+    return 0;
   }
-  return (*(sm_idx_t *)__sm_get_chunk_map_data(map, 0));
+  sm_idx_t *chunk = (sm_idx_t *)__sm_get_chunk_map_data(map, 0);
+  return *chunk;
 }
 
 /**
@@ -1003,10 +1149,13 @@ size_t
 sparsemap_get_size(sparsemap_t *map)
 {
   if (map->m_data_used) {
-    assert(map->m_data_used == __sm_get_size_impl(map));
-    return (map->m_data_used);
+    __sm_when_diag({
+      size_t used = __sm_get_size_impl(map);
+      __sm_assert(map->m_data_used == used);
+    });
+    return map->m_data_used;
   }
-  return (map->m_data_used = __sm_get_size_impl(map));
+  return map->m_data_used = __sm_get_size_impl(map);
 }
 
 /**
@@ -1032,14 +1181,13 @@ sparsemap_scan(sparsemap_t *map, void (*scanner)(sm_idx_t[], size_t), size_t ski
   }
 }
 
-/**
- * Appends all chunk maps from |sstart| to |other|, then reduces the chunk
- * map-count appropriately. |sstart| must be BitVector-aligned!
- */
 void
-sparsemap_split(sparsemap_t *map, size_t sstart, sparsemap_t *other)
+sparsemap_split(sparsemap_t *map, sparsemap_idx_t offset, sparsemap_t *other)
 {
-  assert(sstart % SM_BITS_PER_VECTOR == 0);
+  assert(offset % SM_BITS_PER_VECTOR == 0);
+
+  if (offset < 0)
+    return;
 
   /* |dst| points to the destination buffer */
   uint8_t *dst = __sm_get_chunk_map_end(other);
@@ -1047,9 +1195,9 @@ sparsemap_split(sparsemap_t *map, size_t sstart, sparsemap_t *other)
   /* |src| points to the source-chunk map */
   uint8_t *src = __sm_get_chunk_map_data(map, 0);
 
-  /* |sstart| is relative to the beginning of this sparsemap_t; best
+  /* |offset| is relative to the beginning of this sparsemap_t; best
      make it absolute. */
-  sstart += *(sm_idx_t *)src;
+  offset += *(sm_idx_t *)src;
 
   bool in_middle = false;
   uint8_t *prev = src;
@@ -1058,14 +1206,14 @@ sparsemap_split(sparsemap_t *map, size_t sstart, sparsemap_t *other)
     sm_idx_t start = *(sm_idx_t *)src;
     __sm_chunk_t chunk;
     __sm_chunk_map_init(&chunk, src + sizeof(sm_idx_t));
-    if (start == sstart) {
+    if (start == offset) {
       break;
     }
-    if (start + __sm_chunk_map_get_capacity(&chunk) > sstart) {
+    if (start + __sm_chunk_map_get_capacity(&chunk) > (unsigned long)offset) {
       in_middle = true;
       break;
     }
-    if (start > sstart) {
+    if (start > offset) {
       src = prev;
       i--;
       break;
@@ -1083,12 +1231,12 @@ sparsemap_split(sparsemap_t *map, size_t sstart, sparsemap_t *other)
   /* Now copy all the remaining chunks. */
   int moved = 0;
 
-  /* If |sstart| is in the middle of a chunk then this chunk has to be split */
+  /* If |offset| is in the middle of a chunk then this chunk has to be split */
   if (in_middle) {
     uint8_t buf[sizeof(sm_idx_t) + sizeof(sm_bitvec_t) * 2] = { 0 };
     memcpy(dst, &buf[0], sizeof(buf));
 
-    *(sm_idx_t *)dst = sstart;
+    *(sm_idx_t *)dst = offset;
     dst += sizeof(sm_idx_t);
 
     /* the |other| sparsemap_t now has one additional chunk */
@@ -1104,11 +1252,11 @@ sparsemap_split(sparsemap_t *map, size_t sstart, sparsemap_t *other)
 
     __sm_chunk_t d_chunk;
     __sm_chunk_map_init(&d_chunk, dst);
-    __sm_chunk_map_set_capacity(&d_chunk, capacity - (sstart % capacity));
+    __sm_chunk_map_set_capacity(&d_chunk, capacity - (offset % capacity));
 
     /* Now copy the bits. */
-    size_t d = sstart;
-    for (size_t j = sstart % capacity; j < capacity; j++, d++) {
+    sparsemap_idx_t d = offset;
+    for (size_t j = offset % capacity; j < capacity; j++, d++) {
       if (__sm_chunk_map_is_set(&s_chunk, j)) {
         sparsemap_set(other, d, true);
       }
@@ -1120,7 +1268,7 @@ sparsemap_split(sparsemap_t *map, size_t sstart, sparsemap_t *other)
     i++;
 
     /* Reduce the capacity of the source-chunk map. */
-    __sm_chunk_map_set_capacity(&s_chunk, sstart % capacity);
+    __sm_chunk_map_set_capacity(&s_chunk, offset % capacity);
   }
 
   /* Now continue with all remaining minimaps. */
@@ -1152,93 +1300,141 @@ sparsemap_split(sparsemap_t *map, size_t sstart, sparsemap_t *other)
   assert(sparsemap_get_size(other) > SM_SIZEOF_OVERHEAD);
 }
 
-/**
- * Returns the index of the n'th set bit; uses a 0-based index,
- * i.e. n == 0 for the first bit which is set, n == 1 for the second bit etc.
- */
-size_t
-sparsemap_select(sparsemap_t *map, size_t n)
+sparsemap_idx_t
+sparsemap_select(sparsemap_t *map, sparsemap_idx_t n, bool value)
 {
   assert(sparsemap_get_size(map) >= SM_SIZEOF_OVERHEAD);
   size_t result;
   size_t count = __sm_get_chunk_map_count(map);
-  uint8_t *p = __sm_get_chunk_map_data(map, 0);
+  if (n >= 0) {
+    uint8_t *p = __sm_get_chunk_map_data(map, 0);
 
-  for (size_t i = 0; i < count; i++) {
-    result = *(sm_idx_t *)p;
-    p += sizeof(sm_idx_t);
-    __sm_chunk_t chunk;
-    __sm_chunk_map_init(&chunk, p);
+    for (size_t i = 0; i < count; i++) {
+      result = *(sm_idx_t *)p;
+      p += sizeof(sm_idx_t);
+      __sm_chunk_t chunk;
+      __sm_chunk_map_init(&chunk, p);
 
-    ssize_t new_n = (ssize_t)n;
-    size_t index = __sm_chunk_map_select(&chunk, n, &new_n);
-    if (new_n == -1) {
-      return (result + index);
+      ssize_t new_n = (ssize_t)n;
+      size_t index = __sm_chunk_map_select(&chunk, n, &new_n, value);
+      if (new_n == -1) {
+        return result + index;
+      }
+      n = new_n;
+
+      p += __sm_chunk_map_get_size(&chunk);
     }
-    n = new_n;
-
-    p += __sm_chunk_map_get_size(&chunk);
+    return SPARSEMAP_IDX_MAX; // TODO... shouldn't be here?
+  } else {
+    return SPARSEMAP_IDX_MIN; // TODO... sparsemap_select(map, -n, value); seek from end, not start
   }
-#ifdef DEBUG
-  assert(!"shouldn't be here");
-#endif
-  return (size_t)-1;
 }
 
-/**
- * Counts the set bits starting at 'offset' until and including 'idx', meaning
- * [offset, idx] inclusive.
- */
 size_t
-sparsemap_rank(sparsemap_t *map, size_t offset, size_t idx)
+sparsemap_rank_vec(sparsemap_t *map, size_t x, size_t y, bool value, sm_bitvec_t *vec)
 {
   assert(sparsemap_get_size(map) >= SM_SIZEOF_OVERHEAD);
-  size_t result = 0, prev = 0, count = __sm_get_chunk_map_count(map);
-  uint8_t *p = __sm_get_chunk_map_data(map, 0);
+  size_t amt, gap, pos = 0, result = 0, prev = 0, count;
+  uint8_t *p;
+
+  if (x > y) {
+    return 0;
+  }
+
+  count = __sm_get_chunk_map_count(map);
+
+  if (count == 0) {
+    if (value == false) {
+      /* The count/rank of unset bits in an empty map is inf, so what you requested is the answer. */
+      return y - x + 1;
+    }
+  }
+
+  p = __sm_get_chunk_map_data(map, 0);
 
   for (size_t i = 0; i < count; i++) {
     sm_idx_t start = *(sm_idx_t *)p;
-    if (start > idx) {
-      return (result);
+    /* [prev, start + pos), prev is the last bit examined 0-based. */
+    gap = start - (prev + pos);
+    /* Start of this chunk is greater than the end of the desired range. */
+    if (start > y) {
+      /* This chunk starts after our range [x, y]. */
+      if (value == true) {
+        return result;
+      } else {
+        return result + (y - x) + 1;
+      }
+    } else {
+      /* The range and this chunk overlap. */
+      if (value == false) {
+        if (x > gap) {
+          x -= gap;
+        } else {
+          result += gap - x;
+          x = 0;
+        }
+      } else {
+        if (x > gap) {
+          x -= gap;
+        }
+      }
     }
-    offset -= start - prev;
     prev = start;
     p += sizeof(sm_idx_t);
     __sm_chunk_t chunk;
     __sm_chunk_map_init(&chunk, p);
 
-    result += __sm_chunk_map_rank(&chunk, &offset, idx - start);
+    /* Count all the set/unset inside this chunk. */
+    amt = __sm_chunk_map_rank(&chunk, &x, y - start, &pos, vec, value);
+    result += amt;
     p += __sm_chunk_map_get_size(&chunk);
   }
-  return (result);
+  /* Count any additional unset bits that fall outside the last chunk but
+     within the range. */
+  if (value == false) {
+    size_t last = prev - 1 + pos;
+    if (y > last) {
+      result += y - last - x;
+    }
+  }
+  return result;
 }
 
-/**
- * Finds a span of set bits of at least |len| after |loc|. Returns the index of
- * the n'th set bit that starts a span of at least |len| bits set to true.
- */
 size_t
-sparsemap_span(sparsemap_t *map, size_t loc, size_t len)
+sparsemap_rank(sparsemap_t *map, size_t x, size_t y, bool value)
 {
-  size_t offset, nth = 0, count;
-  (void)loc; // TODO
+  sm_bitvec_t vec;
+  return sparsemap_rank_vec(map, x, y, value, &vec);
+}
 
-  offset = sparsemap_select(map, 0);
+size_t
+sparsemap_span(sparsemap_t *map, sparsemap_idx_t idx, size_t len, bool value)
+{
+  size_t count, nth = 0;
+  sm_bitvec_t vec = 0;
+  sparsemap_idx_t offset;
+
+  offset = sparsemap_select(map, nth++, value);
   if (len == 1) {
     return offset;
   }
   do {
-    count = sparsemap_rank(map, offset, offset + len);
+    count = sparsemap_rank_vec(map, offset, offset + len, value, &vec);
     if (count == len) {
       return offset;
     } else {
-      count = len;
-      while (--count && sparsemap_is_set(map, offset)) {
+      // TODO: what is nth when len > SM_BITS_PER_VECTOR?
+      int c = len > SM_BITS_PER_VECTOR ? SM_BITS_PER_VECTOR : len;
+      for (int b = 0; b < c && (vec & 1 << b); b++) {
         nth++;
       }
     }
-    offset = sparsemap_select(map, nth);
-  } while (offset != ((size_t)-1));
+    if (count) {
+      nth++;
+    }
+    /* Use select to potentially jump very far forward in the map. */
+    offset = sparsemap_select(map, nth, value);
+  } while (offset != SPARSEMAP_IDX_MAX);
 
-  return offset;
+  return idx >= 0 ? SPARSEMAP_IDX_MAX : SPARSEMAP_IDX_MIN;
 }
diff --git a/tests/common.c b/tests/common.c
index 5a632b9..4450f4b 100644
--- a/tests/common.c
+++ b/tests/common.c
@@ -1,13 +1,20 @@
-#include <sys/types.h>
+#define _POSIX_C_SOURCE 199309L
+#define X86_INTRIN
 
 #include <assert.h>
-#include <pthread.h>
-#include <sparsemap.h>
+#include <pthread.h> // If using threads
 #include <stdbool.h>
 #include <stdio.h>
 #include <stdlib.h>
+#include <string.h>
+#include <time.h>
 #include <unistd.h>
+#ifdef X86_INTRIN
+#include <errno.h>
+#include <x86intrin.h>
+#endif
 
+#include "../include/sparsemap.h"
 #include "common.h"
 
 #pragma GCC diagnostic push
@@ -22,84 +29,25 @@
 uint64_t
 tsc(void)
 {
+#ifdef X86_INTRIN
+  return __rdtsc();
+#else
   uint32_t low, high;
   __asm__ volatile("rdtsc" : "=a"(low), "=d"(high));
   return ((uint64_t)high << 32) | low;
-}
-
-static uint64_t
-get_tsc_frequency()
-{
-  uint32_t high, low;
-  __asm__ volatile("rdtsc" : "=a"(low), "=d"(high));
-  __asm__ volatile("rdtsc");
-  return ((uint64_t)high << 32) | low;
+#endif
 }
 
 double
-tsc_ticks_to_ns(uint64_t tsc_ticks)
+nsts()
 {
-  static uint64_t tsc_freq = 0;
-  if (tsc_freq == 0) {
-    tsc_freq = get_tsc_frequency();
-  }
-  return (double)tsc_ticks / (double)tsc_freq * 1e9;
-}
+  struct timespec ts;
 
-void
-est_sift_up(uint64_t *heap, int child_index)
-{
-  while (child_index > 0) {
-    int parent_index = (child_index - 1) / 2;
-    if (heap[parent_index] > heap[child_index]) {
-      // Swap parent and child
-      uint64_t temp = heap[parent_index];
-      heap[parent_index] = heap[child_index];
-      heap[child_index] = temp;
-      child_index = parent_index;
-    } else {
-      break; // Heap property satisfied
-    }
-  }
-}
-
-void
-est_sift_down(uint64_t *heap, int heap_size, int parent_index)
-{
-  int child_index = 2 * parent_index + 1; // Left child
-  while (child_index < heap_size) {
-    // Right child exists and is smaller than left child
-    if (child_index + 1 < heap_size && heap[child_index + 1] < heap[child_index]) {
-      child_index++;
-    }
-    // If the smallest child is smaller than the parent, swap them
-    if (heap[child_index] < heap[parent_index]) {
-      uint64_t temp = heap[child_index];
-      heap[child_index] = heap[parent_index];
-      heap[parent_index] = temp;
-      parent_index = child_index;
-      child_index = 2 * parent_index + 1;
-    } else {
-      break; // Heap property satisfied
-    }
-  }
-}
-
-void
-est_insert_value(uint64_t *heap, int heap_max_size, int *heap_size, uint64_t value)
-{
-  if (*heap_size < heap_max_size) { // Heap not full, insert value
-    heap[*heap_size] = value;
-    est_sift_up(heap, *heap_size);
-    (*heap_size)++;
-  } else {
-    // Heap is full, replace root with new value with a certain probability
-    // This is a very naive approach to maintain a sample of the input
-    if (rand() % 2) {
-      heap[0] = value;
-      est_sift_down(heap, heap_max_size, 0);
-    }
+  if (clock_gettime(CLOCK_REALTIME, &ts) == -1) {
+    perror("clock_gettime");
+    return -1.0; // Return -1.0 on error
   }
+  return ts.tv_sec + ts.tv_nsec / 1e9;
 }
 
 int __xorshift32_state = 0;
@@ -170,7 +118,7 @@ has_sequential_set(int a[], int l, int r)
 int
 ensure_sequential_set(int a[], int l, int r)
 {
-  if (!a || l == 0 || r < 1 || r > l) {
+  if (!a || l == 0 || r < 1 || r > l - 1) {
     return 0;
   }
 
@@ -199,21 +147,6 @@ ensure_sequential_set(int a[], int l, int r)
   return value;
 }
 
-int
-create_sequential_set_in_empty_map(sparsemap_t *map, int s, int r)
-{
-  int placed_at;
-  if (s >= r + 1) {
-    placed_at = 0;
-  } else {
-    placed_at = random_uint32() % (s - r - 1);
-  }
-  for (int i = placed_at; i < placed_at + r; i++) {
-    sparsemap_set(map, i, true);
-  }
-  return placed_at;
-}
-
 void
 print_array(int *array, int l)
 {
@@ -340,6 +273,20 @@ is_unique(int a[], int l, int value)
   return 1; // Unique
 }
 
+int
+whats_set_uint64(uint64_t number, int pos[64])
+{
+  int length = 0;
+
+  for (int i = 0; i < 64; i++) {
+    if (number & ((uint64_t)1 << i)) {
+      pos[length++] = i;
+    }
+  }
+
+  return length;
+}
+
 void
 setup_test_array(int a[], int l, int max_value)
 {
@@ -364,15 +311,6 @@ bitmap_from_uint32(sparsemap_t *map, uint32_t number)
   }
 }
 
-void
-bitmap_from_uint64(sparsemap_t *map, uint64_t number)
-{
-  for (int i = 0; i < 64; i++) {
-    bool bit = number & (1 << i);
-    sparsemap_set(map, i, bit);
-  }
-}
-
 uint32_t
 rank_uint64(uint64_t number, int n, int p)
 {
@@ -400,22 +338,53 @@ rank_uint64(uint64_t number, int n, int p)
   return count;
 }
 
-int
-whats_set_uint64(uint64_t number, int pos[64])
+void
+print_bits(char *name, uint64_t value)
 {
-  int length = 0;
-
-  for (int i = 0; i < 64; i++) {
-    if (number & ((uint64_t)1 << i)) {
-      pos[length++] = i;
+  if (name) {
+    printf("%s\t", name);
+  }
+  for (int i = 63; i >= 0; i--) {
+    printf("%ld", (value >> i) & 1);
+    if (i % 8 == 0) {
+      printf(" "); // Add space for better readability
     }
   }
-
-  return length;
+  printf("\n");
 }
 
 void
-whats_set(sparsemap_t *map, int m)
+sm_bitmap_from_uint64(sparsemap_t *map, uint64_t number)
+{
+  for (int i = 0; i < 64; i++) {
+    bool bit = number & ((uint64_t)1 << i);
+    sparsemap_set(map, i, bit);
+  }
+}
+
+sparsemap_idx_t
+sm_add_span(sparsemap_t *map, int map_size, int span_length)
+{
+  int attempts = map_size / span_length;
+  sparsemap_idx_t placed_at;
+  do {
+    placed_at = random_uint32() % (map_size - span_length - 1);
+    if (sm_occupied(map, placed_at, span_length, true)) {
+      attempts--;
+    } else {
+      break;
+    }
+  } while (attempts);
+  for (int i = placed_at; i < placed_at + span_length; i++) {
+    if (sparsemap_set(map, i, true) != i) {
+      return placed_at; // TODO error?
+    }
+  }
+  return placed_at;
+}
+
+void
+sm_whats_set(sparsemap_t *map, int m)
 {
   logf("what's set in the range [0, %d): ", m);
   for (int i = 0; i < m; i++) {
@@ -425,3 +394,25 @@ whats_set(sparsemap_t *map, int m)
   }
   logf("\n");
 }
+
+bool
+sm_is_span(sparsemap_t *map, sparsemap_idx_t m, int len, bool value)
+{
+  for (sparsemap_idx_t i = m; i < m + len; i++) {
+    if (sparsemap_is_set(map, i) != value) {
+      return false;
+    }
+  }
+  return true;
+}
+
+bool
+sm_occupied(sparsemap_t *map, sparsemap_idx_t m, int len, bool value)
+{
+  for (sparsemap_idx_t i = m; i < (sparsemap_idx_t)len; i++) {
+    if (sparsemap_is_set(map, i) == value) {
+      return true;
+    }
+  }
+  return false;
+}
diff --git a/tests/common.h b/tests/common.h
index 17b27b1..87b75fc 100644
--- a/tests/common.h
+++ b/tests/common.h
@@ -23,20 +23,9 @@
 #define XORSHIFT_SEED_VALUE ((unsigned int)time(NULL) ^ getpid())
 #endif
 
-#define EST_MEDIAN_DECL(decl, size)   \
-  uint64_t heap_##decl[size] = { 0 }; \
-  int heap_##decl##_max_size = size;  \
-  int heap_##decl##_size = 0;
-
-#define EST_MEDIAN_ADD(decl, value) est_insert_value(heap_##decl, heap_##decl##_max_size, &heap_##decl##_size, (value));
-
-#define EST_MEDIAN_GET(decl) heap_##decl[0]
-
 uint64_t tsc(void);
 double tsc_ticks_to_ns(uint64_t tsc_ticks);
-void est_sift_up(uint64_t *heap, int child_index);
-void est_sift_down(uint64_t *heap, int heap_size, int parent_index);
-void est_insert_value(uint64_t *heap, int heap_max_size, int *heap_size, uint64_t value);
+double nsts();
 
 void xorshift32_seed();
 uint32_t xorshift32();
@@ -52,11 +41,16 @@ int is_unique(int a[], int l, int value);
 void setup_test_array(int a[], int l, int max_value);
 void shuffle(int *array, size_t n);
 int ensure_sequential_set(int a[], int l, int r);
-int create_sequential_set_in_empty_map(sparsemap_t *map, int s, int r);
+sparsemap_idx_t sm_add_span(sparsemap_t *map, int map_size, int span_length);
+
+void print_bits(char *name, uint64_t value);
 
 void bitmap_from_uint32(sparsemap_t *map, uint32_t number);
-void bitmap_from_uint64(sparsemap_t *map, uint64_t number);
+void sm_bitmap_from_uint64(sparsemap_t *map, uint64_t number);
 uint32_t rank_uint64(uint64_t number, int n, int p);
 int whats_set_uint64(uint64_t number, int bitPositions[64]);
 
-void whats_set(sparsemap_t *map, int m);
+void sm_whats_set(sparsemap_t *map, int m);
+
+bool sm_is_span(sparsemap_t *map, sparsemap_idx_t m, int len, bool value);
+bool sm_occupied(sparsemap_t *map, sparsemap_idx_t m, int len, bool value);
diff --git a/tests/tdigest.c b/tests/tdigest.c
new file mode 100644
index 0000000..71d3eed
--- /dev/null
+++ b/tests/tdigest.c
@@ -0,0 +1,680 @@
+#include <stdlib.h>
+#include <stdbool.h>
+#include <string.h>
+#include <math.h>
+#include "tdigest.h"
+#include <errno.h>
+#include <stdint.h>
+
+#ifndef TD_MALLOC_INCLUDE
+#define TD_MALLOC_INCLUDE "td_malloc.h"
+#endif
+
+#ifndef TD_ALLOC_H
+#define TD_ALLOC_H
+#define __td_malloc malloc
+#define __td_calloc calloc
+#define __td_realloc realloc
+#define __td_free free
+#endif
+
+#define __td_max(x, y) (((x) > (y)) ? (x) : (y))
+#define __td_min(x, y) (((x) < (y)) ? (x) : (y))
+
+static inline double weighted_average_sorted(double x1, double w1, double x2, double w2) {
+    const double x = (x1 * w1 + x2 * w2) / (w1 + w2);
+    return __td_max(x1, __td_min(x, x2));
+}
+
+static inline bool _tdigest_long_long_add_safe(long long a, long long b) {
+    if (b < 0) {
+        return (a >= __LONG_LONG_MAX__ - b);
+    } else {
+        return (a <= __LONG_LONG_MAX__ - b);
+    }
+}
+
+static inline double weighted_average(double x1, double w1, double x2, double w2) {
+    if (x1 <= x2) {
+        return weighted_average_sorted(x1, w1, x2, w2);
+    } else {
+        return weighted_average_sorted(x2, w2, x1, w1);
+    }
+}
+
+static inline void swap(double *arr, int i, int j) {
+    const double temp = arr[i];
+    arr[i] = arr[j];
+    arr[j] = temp;
+}
+
+static inline void swap_l(long long *arr, int i, int j) {
+    const long long temp = arr[i];
+    arr[i] = arr[j];
+    arr[j] = temp;
+}
+
+static unsigned int partition(double *means, long long *weights, unsigned int start,
+                              unsigned int end, unsigned int pivot_idx) {
+    const double pivotMean = means[pivot_idx];
+    swap(means, pivot_idx, end);
+    swap_l(weights, pivot_idx, end);
+
+    int i = start - 1;
+
+    for (unsigned int j = start; j < end; j++) {
+        // If current element is smaller than the pivot
+        if (means[j] < pivotMean) {
+            // increment index of smaller element
+            i++;
+            swap(means, i, j);
+            swap_l(weights, i, j);
+        }
+    }
+    swap(means, i + 1, end);
+    swap_l(weights, i + 1, end);
+    return i + 1;
+}
+
+/**
+ * Standard quick sort except that sorting rearranges parallel arrays
+ *
+ * @param means  Values to sort on
+ * @param weights The auxillary values to sort.
+ * @param start  The beginning of the values to sort
+ * @param end    The value after the last value to sort
+ */
+static void td_qsort(double *means, long long *weights, unsigned int start, unsigned int end) {
+    if (start < end) {
+        // two elements can be directly compared
+        if ((end - start) == 1) {
+            if (means[start] > means[end]) {
+                swap(means, start, end);
+                swap_l(weights, start, end);
+            }
+            return;
+        }
+        // generating a random number as a pivot was very expensive vs the array size
+        // const unsigned int pivot_idx = start + rand()%(end - start + 1);
+        const unsigned int pivot_idx = (end + start) / 2; // central pivot
+        const unsigned int new_pivot_idx = partition(means, weights, start, end, pivot_idx);
+        if (new_pivot_idx > start) {
+            td_qsort(means, weights, start, new_pivot_idx - 1);
+        }
+        td_qsort(means, weights, new_pivot_idx + 1, end);
+    }
+}
+
+static inline size_t cap_from_compression(double compression) {
+    if ((size_t)compression > ((SIZE_MAX / sizeof(double) / 6) - 10)) {
+        return 0;
+    }
+
+    return (6 * (size_t)(compression)) + 10;
+}
+
+static inline bool should_td_compress(td_histogram_t *h) {
+    return ((h->merged_nodes + h->unmerged_nodes) >= (h->cap - 1));
+}
+
+static inline int next_node(td_histogram_t *h) { return h->merged_nodes + h->unmerged_nodes; }
+
+int td_compress(td_histogram_t *h);
+
+static inline int _check_overflow(const double v) {
+    // double-precision overflow detected on h->unmerged_weight
+    if (v == INFINITY) {
+        return EDOM;
+    }
+    return 0;
+}
+
+static inline int _check_td_overflow(const double new_unmerged_weight,
+                                     const double new_total_weight) {
+    // double-precision overflow detected on h->unmerged_weight
+    if (new_unmerged_weight == INFINITY) {
+        return EDOM;
+    }
+    if (new_total_weight == INFINITY) {
+        return EDOM;
+    }
+    const double denom = 2 * MM_PI * new_total_weight * log(new_total_weight);
+    if (denom == INFINITY) {
+        return EDOM;
+    }
+
+    return 0;
+}
+
+int td_centroid_count(td_histogram_t *h) { return next_node(h); }
+
+void td_reset(td_histogram_t *h) {
+    if (!h) {
+        return;
+    }
+    h->min = __DBL_MAX__;
+    h->max = -h->min;
+    h->merged_nodes = 0;
+    h->merged_weight = 0;
+    h->unmerged_nodes = 0;
+    h->unmerged_weight = 0;
+    h->total_compressions = 0;
+}
+
+int td_init(double compression, td_histogram_t **result) {
+
+    const size_t capacity = cap_from_compression(compression);
+    if (capacity < 1) {
+        return 1;
+    }
+    td_histogram_t *histogram;
+    histogram = (td_histogram_t *)__td_malloc(sizeof(td_histogram_t));
+    if (!histogram) {
+        return 1;
+    }
+    histogram->cap = capacity;
+    histogram->compression = (double)compression;
+    td_reset(histogram);
+    histogram->nodes_mean = (double *)__td_calloc(capacity, sizeof(double));
+    if (!histogram->nodes_mean) {
+        td_free(histogram);
+        return 1;
+    }
+    histogram->nodes_weight = (long long *)__td_calloc(capacity, sizeof(long long));
+    if (!histogram->nodes_weight) {
+        td_free(histogram);
+        return 1;
+    }
+    *result = histogram;
+
+    return 0;
+}
+
+td_histogram_t *td_new(double compression) {
+    td_histogram_t *mdigest = NULL;
+    td_init(compression, &mdigest);
+    return mdigest;
+}
+
+void td_free(td_histogram_t *histogram) {
+    if (histogram->nodes_mean) {
+        __td_free((void *)(histogram->nodes_mean));
+    }
+    if (histogram->nodes_weight) {
+        __td_free((void *)(histogram->nodes_weight));
+    }
+    __td_free((void *)(histogram));
+}
+
+int td_merge(td_histogram_t *into, td_histogram_t *from) {
+    if (td_compress(into) != 0)
+        return EDOM;
+    if (td_compress(from) != 0)
+        return EDOM;
+    const int pos = from->merged_nodes + from->unmerged_nodes;
+    for (int i = 0; i < pos; i++) {
+        const double mean = from->nodes_mean[i];
+        const long long weight = from->nodes_weight[i];
+        if (td_add(into, mean, weight) != 0) {
+            return EDOM;
+        }
+    }
+    return 0;
+}
+
+long long td_size(td_histogram_t *h) { return h->merged_weight + h->unmerged_weight; }
+
+double td_cdf(td_histogram_t *h, double val) {
+    td_compress(h);
+    // no data to examine
+    if (h->merged_nodes == 0) {
+        return NAN;
+    }
+    // bellow lower bound
+    if (val < h->min) {
+        return 0;
+    }
+    // above upper bound
+    if (val > h->max) {
+        return 1;
+    }
+    if (h->merged_nodes == 1) {
+        // exactly one centroid, should have max==min
+        const double width = h->max - h->min;
+        if (val - h->min <= width) {
+            // min and max are too close together to do any viable interpolation
+            return 0.5;
+        } else {
+            // interpolate if somehow we have weight > 0 and max != min
+            return (val - h->min) / width;
+        }
+    }
+    const int n = h->merged_nodes;
+    // check for the left tail
+    const double left_centroid_mean = h->nodes_mean[0];
+    const double left_centroid_weight = (double)h->nodes_weight[0];
+    const double merged_weight_d = (double)h->merged_weight;
+    if (val < left_centroid_mean) {
+        // note that this is different than h->nodes_mean[0] > min
+        // ... this guarantees we divide by non-zero number and interpolation works
+        const double width = left_centroid_mean - h->min;
+        if (width > 0) {
+            // must be a sample exactly at min
+            if (val == h->min) {
+                return 0.5 / merged_weight_d;
+            } else {
+                return (1 + (val - h->min) / width * (left_centroid_weight / 2 - 1)) /
+                       merged_weight_d;
+            }
+        } else {
+            // this should be redundant with the check val < h->min
+            return 0;
+        }
+    }
+    // and the right tail
+    const double right_centroid_mean = h->nodes_mean[n - 1];
+    const double right_centroid_weight = (double)h->nodes_weight[n - 1];
+    if (val > right_centroid_mean) {
+        const double width = h->max - right_centroid_mean;
+        if (width > 0) {
+            if (val == h->max) {
+                return 1 - 0.5 / merged_weight_d;
+            } else {
+                // there has to be a single sample exactly at max
+                const double dq = (1 + (h->max - val) / width * (right_centroid_weight / 2 - 1)) /
+                                  merged_weight_d;
+                return 1 - dq;
+            }
+        } else {
+            return 1;
+        }
+    }
+    // we know that there are at least two centroids and mean[0] < x < mean[n-1]
+    // that means that there are either one or more consecutive centroids all at exactly x
+    // or there are consecutive centroids, c0 < x < c1
+    double weightSoFar = 0;
+    for (int it = 0; it < n - 1; it++) {
+        // weightSoFar does not include weight[it] yet
+        if (h->nodes_mean[it] == val) {
+            // we have one or more centroids == x, treat them as one
+            // dw will accumulate the weight of all of the centroids at x
+            double dw = 0;
+            while (it < n && h->nodes_mean[it] == val) {
+                dw += (double)h->nodes_weight[it];
+                it++;
+            }
+            return (weightSoFar + dw / 2) / (double)h->merged_weight;
+        } else if (h->nodes_mean[it] <= val && val < h->nodes_mean[it + 1]) {
+            const double node_weight = (double)h->nodes_weight[it];
+            const double node_weight_next = (double)h->nodes_weight[it + 1];
+            const double node_mean = h->nodes_mean[it];
+            const double node_mean_next = h->nodes_mean[it + 1];
+            // landed between centroids ... check for floating point madness
+            if (node_mean_next - node_mean > 0) {
+                // note how we handle singleton centroids here
+                // the point is that for singleton centroids, we know that their entire
+                // weight is exactly at the centroid and thus shouldn't be involved in
+                // interpolation
+                double leftExcludedW = 0;
+                double rightExcludedW = 0;
+                if (node_weight == 1) {
+                    if (node_weight_next == 1) {
+                        // two singletons means no interpolation
+                        // left singleton is in, right is out
+                        return (weightSoFar + 1) / merged_weight_d;
+                    } else {
+                        leftExcludedW = 0.5;
+                    }
+                } else if (node_weight_next == 1) {
+                    rightExcludedW = 0.5;
+                }
+                double dw = (node_weight + node_weight_next) / 2;
+
+                // adjust endpoints for any singleton
+                double dwNoSingleton = dw - leftExcludedW - rightExcludedW;
+
+                double base = weightSoFar + node_weight / 2 + leftExcludedW;
+                return (base + dwNoSingleton * (val - node_mean) / (node_mean_next - node_mean)) /
+                       merged_weight_d;
+            } else {
+                // this is simply caution against floating point madness
+                // it is conceivable that the centroids will be different
+                // but too near to allow safe interpolation
+                double dw = (node_weight + node_weight_next) / 2;
+                return (weightSoFar + dw) / merged_weight_d;
+            }
+        } else {
+            weightSoFar += (double)h->nodes_weight[it];
+        }
+    }
+    return 1 - 0.5 / merged_weight_d;
+}
+
+static double td_internal_iterate_centroids_to_index(const td_histogram_t *h, const double index,
+                                                     const double left_centroid_weight,
+                                                     const int total_centroids, double *weightSoFar,
+                                                     int *node_pos) {
+    if (left_centroid_weight > 1 && index < left_centroid_weight / 2) {
+        // there is a single sample at min so we interpolate with less weight
+        return h->min + (index - 1) / (left_centroid_weight / 2 - 1) * (h->nodes_mean[0] - h->min);
+    }
+
+    // usually the last centroid will have unit weight so this test will make it moot
+    if (index > h->merged_weight - 1) {
+        return h->max;
+    }
+
+    // if the right-most centroid has more than one sample, we still know
+    // that one sample occurred at max so we can do some interpolation
+    const double right_centroid_weight = (double)h->nodes_weight[total_centroids - 1];
+    const double right_centroid_mean = h->nodes_mean[total_centroids - 1];
+    if (right_centroid_weight > 1 &&
+        (double)h->merged_weight - index <= right_centroid_weight / 2) {
+        return h->max - ((double)h->merged_weight - index - 1) / (right_centroid_weight / 2 - 1) *
+                            (h->max - right_centroid_mean);
+    }
+
+    for (; *node_pos < total_centroids - 1; (*node_pos)++) {
+        const int i = *node_pos;
+        const double node_weight = (double)h->nodes_weight[i];
+        const double node_weight_next = (double)h->nodes_weight[i + 1];
+        const double node_mean = h->nodes_mean[i];
+        const double node_mean_next = h->nodes_mean[i + 1];
+        const double dw = (node_weight + node_weight_next) / 2;
+        if (*weightSoFar + dw > index) {
+            // centroids i and i+1 bracket our current point
+            // check for unit weight
+            double leftUnit = 0;
+            if (node_weight == 1) {
+                if (index - *weightSoFar < 0.5) {
+                    // within the singleton's sphere
+                    return node_mean;
+                } else {
+                    leftUnit = 0.5;
+                }
+            }
+            double rightUnit = 0;
+            if (node_weight_next == 1) {
+                if (*weightSoFar + dw - index <= 0.5) {
+                    // no interpolation needed near singleton
+                    return node_mean_next;
+                }
+                rightUnit = 0.5;
+            }
+            const double z1 = index - *weightSoFar - leftUnit;
+            const double z2 = *weightSoFar + dw - index - rightUnit;
+            return weighted_average(node_mean, z2, node_mean_next, z1);
+        }
+        *weightSoFar += dw;
+    }
+
+    // weightSoFar = totalWeight - weight[total_centroids-1]/2 (very nearly)
+    // so we interpolate out to max value ever seen
+    const double z1 = index - h->merged_weight - right_centroid_weight / 2.0;
+    const double z2 = right_centroid_weight / 2 - z1;
+    return weighted_average(right_centroid_mean, z1, h->max, z2);
+}
+
+double td_quantile(td_histogram_t *h, double q) {
+    td_compress(h);
+    // q should be in [0,1]
+    if (q < 0.0 || q > 1.0 || h->merged_nodes == 0) {
+        return NAN;
+    }
+    // with one data point, all quantiles lead to Rome
+    if (h->merged_nodes == 1) {
+        return h->nodes_mean[0];
+    }
+
+    // if values were stored in a sorted array, index would be the offset we are interested in
+    const double index = q * (double)h->merged_weight;
+
+    // beyond the boundaries, we return min or max
+    // usually, the first centroid will have unit weight so this will make it moot
+    if (index < 1) {
+        return h->min;
+    }
+
+    // we know that there are at least two centroids now
+    const int n = h->merged_nodes;
+
+    // if the left centroid has more than one sample, we still know
+    // that one sample occurred at min so we can do some interpolation
+    const double left_centroid_weight = (double)h->nodes_weight[0];
+
+    // in between extremes we interpolate between centroids
+    double weightSoFar = left_centroid_weight / 2;
+    int i = 0;
+    return td_internal_iterate_centroids_to_index(h, index, left_centroid_weight, n, &weightSoFar,
+                                                  &i);
+}
+
+int td_quantiles(td_histogram_t *h, const double *quantiles, double *values, size_t length) {
+    td_compress(h);
+
+    if (NULL == quantiles || NULL == values) {
+        return EINVAL;
+    }
+
+    const int n = h->merged_nodes;
+    if (n == 0) {
+        for (size_t i = 0; i < length; i++) {
+            values[i] = NAN;
+        }
+        return 0;
+    }
+    if (n == 1) {
+        for (size_t i = 0; i < length; i++) {
+            const double requested_quantile = quantiles[i];
+
+            // q should be in [0,1]
+            if (requested_quantile < 0.0 || requested_quantile > 1.0) {
+                values[i] = NAN;
+            } else {
+                // with one data point, all quantiles lead to Rome
+                values[i] = h->nodes_mean[0];
+            }
+        }
+        return 0;
+    }
+
+    // we know that there are at least two centroids now
+    // if the left centroid has more than one sample, we still know
+    // that one sample occurred at min so we can do some interpolation
+    const double left_centroid_weight = (double)h->nodes_weight[0];
+
+    // in between extremes we interpolate between centroids
+    double weightSoFar = left_centroid_weight / 2;
+    int node_pos = 0;
+
+    // to avoid allocations we use the values array for intermediate computation
+    // i.e. to store the expected cumulative count at each percentile
+    for (size_t qpos = 0; qpos < length; qpos++) {
+        const double index = quantiles[qpos] * (double)h->merged_weight;
+        values[qpos] = td_internal_iterate_centroids_to_index(h, index, left_centroid_weight, n,
+                                                              &weightSoFar, &node_pos);
+    }
+    return 0;
+}
+
+static double td_internal_trimmed_mean(const td_histogram_t *h, const double leftmost_weight,
+                                       const double rightmost_weight) {
+    double count_done = 0;
+    double trimmed_sum = 0;
+    double trimmed_count = 0;
+    for (int i = 0; i < h->merged_nodes; i++) {
+
+        const double n_weight = (double)h->nodes_weight[i];
+        // Assume the whole centroid falls into the range
+        double count_add = n_weight;
+
+        // If we haven't reached the low threshold yet, skip appropriate part of the centroid.
+        count_add -= __td_min(__td_max(0, leftmost_weight - count_done), count_add);
+
+        // If we have reached the upper threshold, ignore the overflowing part of the centroid.
+
+        count_add = __td_min(__td_max(0, rightmost_weight - count_done), count_add);
+
+        // consider the whole centroid processed
+        count_done += n_weight;
+
+        // increment the sum / count
+        trimmed_sum += h->nodes_mean[i] * count_add;
+        trimmed_count += count_add;
+
+        // break once we cross the high threshold
+        if (count_done >= rightmost_weight)
+            break;
+    }
+
+    return trimmed_sum / trimmed_count;
+}
+
+double td_trimmed_mean_symmetric(td_histogram_t *h, double proportion_to_cut) {
+    td_compress(h);
+    // proportion_to_cut should be in [0,1]
+    if (h->merged_nodes == 0 || proportion_to_cut < 0.0 || proportion_to_cut > 1.0) {
+        return NAN;
+    }
+    // with one data point, all values lead to Rome
+    if (h->merged_nodes == 1) {
+        return h->nodes_mean[0];
+    }
+
+    /* translate the percentiles to counts */
+    const double leftmost_weight = floor((double)h->merged_weight * proportion_to_cut);
+    const double rightmost_weight = ceil((double)h->merged_weight * (1.0 - proportion_to_cut));
+
+    return td_internal_trimmed_mean(h, leftmost_weight, rightmost_weight);
+}
+
+double td_trimmed_mean(td_histogram_t *h, double leftmost_cut, double rightmost_cut) {
+    td_compress(h);
+    // leftmost_cut and rightmost_cut should be in [0,1]
+    if (h->merged_nodes == 0 || leftmost_cut < 0.0 || leftmost_cut > 1.0 || rightmost_cut < 0.0 ||
+        rightmost_cut > 1.0) {
+        return NAN;
+    }
+    // with one data point, all values lead to Rome
+    if (h->merged_nodes == 1) {
+        return h->nodes_mean[0];
+    }
+
+    /* translate the percentiles to counts */
+    const double leftmost_weight = floor((double)h->merged_weight * leftmost_cut);
+    const double rightmost_weight = ceil((double)h->merged_weight * rightmost_cut);
+
+    return td_internal_trimmed_mean(h, leftmost_weight, rightmost_weight);
+}
+
+int td_add(td_histogram_t *h, double mean, long long weight) {
+    if (should_td_compress(h)) {
+        const int overflow_res = td_compress(h);
+        if (overflow_res != 0)
+            return overflow_res;
+    }
+    const int pos = next_node(h);
+    if (pos >= h->cap)
+        return EDOM;
+    if (_tdigest_long_long_add_safe(h->unmerged_weight, weight) == false)
+        return EDOM;
+    const long long new_unmerged_weight = h->unmerged_weight + weight;
+    if (_tdigest_long_long_add_safe(new_unmerged_weight, h->merged_weight) == false)
+        return EDOM;
+    const long long new_total_weight = new_unmerged_weight + h->merged_weight;
+    // double-precision overflow detected
+    const int overflow_res =
+        _check_td_overflow((double)new_unmerged_weight, (double)new_total_weight);
+    if (overflow_res != 0)
+        return overflow_res;
+
+    if (mean < h->min) {
+        h->min = mean;
+    }
+    if (mean > h->max) {
+        h->max = mean;
+    }
+    h->nodes_mean[pos] = mean;
+    h->nodes_weight[pos] = weight;
+    h->unmerged_nodes++;
+    h->unmerged_weight = new_unmerged_weight;
+    return 0;
+}
+
+int td_compress(td_histogram_t *h) {
+    if (h->unmerged_nodes == 0) {
+        return 0;
+    }
+    int N = h->merged_nodes + h->unmerged_nodes;
+    td_qsort(h->nodes_mean, h->nodes_weight, 0, N - 1);
+    const double total_weight = (double)h->merged_weight + (double)h->unmerged_weight;
+    // double-precision overflow detected
+    const int overflow_res = _check_td_overflow((double)h->unmerged_weight, (double)total_weight);
+    if (overflow_res != 0)
+        return overflow_res;
+    if (total_weight <= 1)
+        return 0;
+    const double denom = 2 * MM_PI * total_weight * log(total_weight);
+    if (_check_overflow(denom) != 0)
+        return EDOM;
+
+    // Compute the normalizer given compression and number of points.
+    const double normalizer = h->compression / denom;
+    if (_check_overflow(normalizer) != 0)
+        return EDOM;
+    int cur = 0;
+    double weight_so_far = 0;
+
+    for (int i = 1; i < N; i++) {
+        const double proposed_weight = (double)h->nodes_weight[cur] + (double)h->nodes_weight[i];
+        const double z = proposed_weight * normalizer;
+        // quantile up to cur
+        const double q0 = weight_so_far / total_weight;
+        // quantile up to cur + i
+        const double q2 = (weight_so_far + proposed_weight) / total_weight;
+        // Convert  a quantile to the k-scale
+        const bool should_add = (z <= (q0 * (1 - q0))) && (z <= (q2 * (1 - q2)));
+        // next point will fit
+        // so merge into existing centroid
+        if (should_add) {
+            h->nodes_weight[cur] += h->nodes_weight[i];
+            const double delta = h->nodes_mean[i] - h->nodes_mean[cur];
+            const double weighted_delta = (delta * h->nodes_weight[i]) / h->nodes_weight[cur];
+            h->nodes_mean[cur] += weighted_delta;
+        } else {
+            weight_so_far += h->nodes_weight[cur];
+            cur++;
+            h->nodes_weight[cur] = h->nodes_weight[i];
+            h->nodes_mean[cur] = h->nodes_mean[i];
+        }
+        if (cur != i) {
+            h->nodes_weight[i] = 0;
+            h->nodes_mean[i] = 0.0;
+        }
+    }
+    h->merged_nodes = cur + 1;
+    h->merged_weight = total_weight;
+    h->unmerged_nodes = 0;
+    h->unmerged_weight = 0;
+    h->total_compressions++;
+    return 0;
+}
+
+double td_min(td_histogram_t *h) { return h->min; }
+
+double td_max(td_histogram_t *h) { return h->max; }
+
+int td_compression(td_histogram_t *h) { return h->compression; }
+
+const long long *td_centroids_weight(td_histogram_t *h) { return h->nodes_weight; }
+
+const double *td_centroids_mean(td_histogram_t *h) { return h->nodes_mean; }
+
+long long td_centroids_weight_at(td_histogram_t *h, int pos) { return h->nodes_weight[pos]; }
+
+double td_centroids_mean_at(td_histogram_t *h, int pos) {
+    if (pos < 0 || pos > h->merged_nodes) {
+        return NAN;
+    }
+    return h->nodes_mean[pos];
+}
diff --git a/tests/tdigest.h b/tests/tdigest.h
new file mode 100644
index 0000000..c07436c
--- /dev/null
+++ b/tests/tdigest.h
@@ -0,0 +1,258 @@
+#pragma once
+#include <stdlib.h>
+
+/**
+ * Adaptive histogram based on something like streaming k-means crossed with Q-digest.
+ * The implementation is a direct descendent of MergingDigest
+ * https://github.com/tdunning/t-digest/
+ *
+ * Copyright (c) 2021 Redis, All rights reserved.
+ * Copyright (c) 2018 Andrew Werner, All rights reserved.
+ *
+ * The special characteristics of this algorithm are:
+ *
+ * - smaller summaries than Q-digest
+ *
+ * - provides part per million accuracy for extreme quantiles and typically &lt;1000 ppm accuracy
+ * for middle quantiles
+ *
+ * - fast
+ *
+ * - simple
+ *
+ * - easy to adapt for use with map-reduce
+ */
+
+#define MM_PI 3.14159265358979323846
+
+struct td_histogram {
+    // compression is a setting used to configure the size of centroids when merged.
+    double compression;
+
+    double min;
+    double max;
+
+    // cap is the total size of nodes
+    int cap;
+    // merged_nodes is the number of merged nodes at the front of nodes.
+    int merged_nodes;
+    // unmerged_nodes is the number of buffered nodes.
+    int unmerged_nodes;
+
+    // we run the merge in reverse every other merge to avoid left-to-right bias in merging
+    long long total_compressions;
+
+    long long merged_weight;
+    long long unmerged_weight;
+
+    double *nodes_mean;
+    long long *nodes_weight;
+};
+
+typedef struct td_histogram td_histogram_t;
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/**
+ * Allocate the memory, initialise the t-digest, and return the histogram as output parameter.
+ * @param compression The compression parameter.
+ * 100 is a common value for normal uses.
+ * 1000 is extremely large.
+ * The number of centroids retained will be a smallish (usually less than 10) multiple of this
+ * number.
+ * @return the histogram on success, NULL if allocation failed.
+ */
+td_histogram_t *td_new(double compression);
+
+/**
+ * Allocate the memory and initialise the t-digest.
+ *
+ * @param compression The compression parameter.
+ * 100 is a common value for normal uses.
+ * 1000 is extremely large.
+ * The number of centroids retained will be a smallish (usually less than 10) multiple of this
+ * number.
+ * @param result Output parameter to capture allocated histogram.
+ * @return 0 on success, 1 if allocation failed.
+ */
+int td_init(double compression, td_histogram_t **result);
+
+/**
+ * Frees the memory associated with the t-digest.
+ *
+ * @param h The histogram you want to free.
+ */
+void td_free(td_histogram_t *h);
+
+/**
+ * Reset a histogram to zero - empty out a histogram and re-initialise it
+ *
+ * If you want to re-use an existing histogram, but reset everything back to zero, this
+ * is the routine to use.
+ *
+ * @param h The histogram you want to reset to empty.
+ *
+ */
+void td_reset(td_histogram_t *h);
+
+/**
+ * Adds a sample to a histogram.
+ *
+ * @param val The value to add.
+ * @param weight The weight of this point.
+ * @return 0 on success, EDOM if overflow was detected as a consequence of adding the provided
+ * weight.
+ *
+ */
+int td_add(td_histogram_t *h, double val, long long weight);
+
+/**
+ * Re-examines a t-digest to determine whether some centroids are redundant.  If your data are
+ * perversely ordered, this may be a good idea.  Even if not, this may save 20% or so in space.
+ *
+ * The cost is roughly the same as adding as many data points as there are centroids.  This
+ * is typically &lt; 10 * compression, but could be as high as 100 * compression.
+ * This is a destructive operation that is not thread-safe.
+ *
+ * @param h The histogram you want to compress.
+ * @return 0 on success, EDOM if overflow was detected as a consequence of adding the provided
+ * weight. If overflow is detected the histogram is not changed.
+ *
+ */
+int td_compress(td_histogram_t *h);
+
+/**
+ * Merges all of the values from 'from' to 'this' histogram.
+ *
+ * @param h "This" pointer
+ * @param from Histogram to copy values from.
+ * * @return 0 on success, EDOM if overflow was detected as a consequence of merging the the
+ * provided histogram. If overflow is detected the original histogram is not detected.
+ */
+int td_merge(td_histogram_t *h, td_histogram_t *from);
+
+/**
+ * Returns the fraction of all points added which are &le; x.
+ *
+ * @param x The cutoff for the cdf.
+ * @return The fraction of all data which is less or equal to x.
+ */
+double td_cdf(td_histogram_t *h, double x);
+
+/**
+ * Returns an estimate of the cutoff such that a specified fraction of the data
+ * added to this TDigest would be less than or equal to the cutoff.
+ *
+ * @param q The desired fraction
+ * @return The value x such that cdf(x) == q;
+ */
+double td_quantile(td_histogram_t *h, double q);
+
+/**
+ * Returns an estimate of the cutoff such that a specified fraction of the data
+ * added to this TDigest would be less than or equal to the cutoffs.
+ *
+ * @param quantiles The ordered percentiles array to get the values for.
+ * @param values Destination array containing the values at the given quantiles.
+ * The values array should be allocated by the caller.
+ * @return 0 on success, ENOMEM if the provided destination array is null.
+ */
+int td_quantiles(td_histogram_t *h, const double *quantiles, double *values, size_t length);
+
+/**
+ * Returns the trimmed mean ignoring values outside given cutoff upper and lower limits.
+ *
+ * @param leftmost_cut Fraction to cut off of the left tail of the distribution.
+ * @param rightmost_cut Fraction to cut off of the right tail of the distribution.
+ * @return The trimmed mean ignoring values outside given cutoff upper and lower limits;
+ */
+double td_trimmed_mean(td_histogram_t *h, double leftmost_cut, double rightmost_cut);
+
+/**
+ * Returns the trimmed mean ignoring values outside given a symmetric cutoff limits.
+ *
+ * @param proportion_to_cut Fraction to cut off of the left and right tails of the distribution.
+ * @return The trimmed mean ignoring values outside given cutoff upper and lower limits;
+ */
+double td_trimmed_mean_symmetric(td_histogram_t *h, double proportion_to_cut);
+
+/**
+ * Returns the current compression factor.
+ *
+ * @return The compression factor originally used to set up the TDigest.
+ */
+int td_compression(td_histogram_t *h);
+
+/**
+ * Returns the number of points that have been added to this TDigest.
+ *
+ * @return The sum of the weights on all centroids.
+ */
+long long td_size(td_histogram_t *h);
+
+/**
+ * Returns the number of centroids being used by this TDigest.
+ *
+ * @return The number of centroids being used.
+ */
+int td_centroid_count(td_histogram_t *h);
+
+/**
+ * Get minimum value from the histogram.  Will return __DBL_MAX__ if the histogram
+ * is empty.
+ *
+ * @param h "This" pointer
+ */
+double td_min(td_histogram_t *h);
+
+/**
+ * Get maximum value from the histogram.  Will return - __DBL_MAX__ if the histogram
+ * is empty.
+ *
+ * @param h "This" pointer
+ */
+double td_max(td_histogram_t *h);
+
+/**
+ * Get the full centroids weight array for 'this' histogram.
+ *
+ * @param h "This" pointer
+ *
+ * @return The full centroids weight array.
+ */
+const long long *td_centroids_weight(td_histogram_t *h);
+
+/**
+ * Get the full centroids mean array for 'this' histogram.
+ *
+ * @param h "This" pointer
+ *
+ * @return The full centroids mean array.
+ */
+const double *td_centroids_mean(td_histogram_t *h);
+
+/**
+ * Get the centroid weight for 'this' histogram and 'pos'.
+ *
+ * @param h "This" pointer
+ * @param pos centroid position.
+ *
+ * @return The centroid weight.
+ */
+long long td_centroids_weight_at(td_histogram_t *h, int pos);
+
+/**
+ * Get the centroid mean for 'this' histogram and 'pos'.
+ *
+ * @param h "This" pointer
+ * @param pos centroid position.
+ *
+ * @return The centroid mean.
+ */
+double td_centroids_mean_at(td_histogram_t *h, int pos);
+
+#ifdef __cplusplus
+}
+#endif
diff --git a/tests/test.c b/tests/test.c
index a32c1d6..f52d833 100644
--- a/tests/test.c
+++ b/tests/test.c
@@ -10,6 +10,8 @@
 #define MUNIT_NO_FORK (1)
 #define MUNIT_ENABLE_ASSERT_ALIASES (1)
 
+#include <errno.h>
+#include <stdio.h>
 #include <stdlib.h>
 #include <unistd.h>
 
@@ -17,21 +19,27 @@
 #include "common.h"
 #include "munit.h"
 
+#define munit_free free
+
 #if defined(_MSC_VER)
 #pragma warning(disable : 4127)
 #endif
 
+#define SELECT_FALSE
+
 /* !!! Duplicated here for testing purposes. Keep in sync, or suffer. !!! */
 struct sparsemap {
-  uint8_t *m_data;
   size_t m_capacity;
   size_t m_data_used;
+  uint8_t *m_data;
 };
 
 struct user_data {
   int foo;
 };
 
+/* -------------------------- Supporting Functions for Testing */
+
 void
 populate_map(sparsemap_t *map, int size, int max_value)
 {
@@ -42,7 +50,8 @@ populate_map(sparsemap_t *map, int size, int max_value)
   shuffle(array, size);
   for (int i = 0; i < size; i++) {
     sparsemap_set(map, array[i], true);
-    munit_assert_true(sparsemap_is_set(map, array[i]));
+    bool set = sparsemap_is_set(map, array[i]);
+    munit_assert_true(set);
   }
 }
 
@@ -62,11 +71,96 @@ test_api_tear_down(void *fixture)
 {
   sparsemap_t *map = (sparsemap_t *)fixture;
   assert_ptr_not_null(map);
-  free(map);
+  munit_free(map);
+}
+
+/* -------------------------- API Tests */
+
+static MunitResult
+test_api_new(const MunitParameter params[], void *data)
+{
+  sparsemap_t *map = sparsemap(1024);
+  (void)params;
+  (void)data;
+
+  assert_ptr_not_null(map);
+  assert_true(map->m_capacity == 1024);
+  assert_true(map->m_data_used == sizeof(uint32_t));
+
+  munit_free(map);
+
+  return MUNIT_OK;
 }
 
 static MunitResult
-test_api_static_init(const MunitParameter params[], void *data)
+test_api_new_realloc(const MunitParameter params[], void *data)
+{
+  sparsemap_t *map = sparsemap(1024);
+  (void)params;
+  (void)data;
+
+  assert_ptr_not_null(map);
+  assert_true(map->m_capacity == 1024);
+  assert_true(map->m_data_used == sizeof(uint32_t));
+
+  map = sparsemap_set_data_size(map, 2048);
+  assert_true(map->m_capacity == 2048);
+  assert_true(map->m_data_used == sizeof(uint32_t));
+
+  munit_free(map);
+
+  return MUNIT_OK;
+}
+
+static MunitResult
+test_api_new_heap(const MunitParameter params[], void *data)
+{
+  sparsemap_t *map;
+  uint8_t *buf;
+  (void)params;
+  (void)data;
+
+  map = munit_malloc(sizeof(sparsemap_t));
+  assert_ptr_not_null(map);
+  buf = munit_calloc(1024, sizeof(uint8_t));
+  assert_ptr_not_null(buf);
+  sparsemap_init(map, buf, 1024);
+
+  sparsemap_init(map, buf, 1024);
+  assert_ptr_equal(buf, map->m_data);
+  assert_true(map->m_capacity == 1024);
+  assert_true(map->m_data_used == sizeof(uint32_t));
+
+  munit_free(map->m_data);
+  munit_free(map);
+
+  return MUNIT_OK;
+}
+
+static MunitResult
+test_api_new_static(const MunitParameter params[], void *data)
+{
+  sparsemap_t a_map, *map = &a_map;
+  uint8_t *buf;
+  (void)params;
+  (void)data;
+
+  buf = munit_calloc(1024, sizeof(uint8_t));
+  assert_ptr_not_null(buf);
+  sparsemap_init(map, buf, 1024);
+
+  sparsemap_init(map, buf, 1024);
+  assert_ptr_equal(buf, map->m_data);
+  assert_true(map->m_capacity == 1024);
+  assert_true(map->m_data_used == sizeof(uint32_t));
+
+  munit_free(map->m_data);
+
+  return MUNIT_OK;
+}
+
+static MunitResult
+test_api_new_stack(const MunitParameter params[], void *data)
 {
   sparsemap_t a_map, *map = &a_map;
   uint8_t buf[1024] = { 0 };
@@ -74,7 +168,6 @@ test_api_static_init(const MunitParameter params[], void *data)
   (void)params;
   (void)data;
 
-  assert_ptr_not_null(map);
   sparsemap_init(map, buf, 1024);
   assert_ptr_equal(&buf, map->m_data);
   assert_true(map->m_capacity == 1024);
@@ -87,6 +180,7 @@ static void *
 test_api_clear_setup(const MunitParameter params[], void *user_data)
 {
   uint8_t *buf = munit_calloc(1024, sizeof(uint8_t));
+  assert_ptr_not_null(buf);
   sparsemap_t *map = (sparsemap_t *)test_api_setup(params, user_data);
 
   sparsemap_init(map, buf, 1024);
@@ -97,7 +191,8 @@ static void
 test_api_clear_tear_down(void *fixture)
 {
   sparsemap_t *map = (sparsemap_t *)fixture;
-  free(map->m_data);
+  assert_ptr_not_null(map->m_data);
+  munit_free(map->m_data);
   test_api_tear_down(fixture);
 }
 static MunitResult
@@ -120,6 +215,7 @@ static void *
 test_api_open_setup(const MunitParameter params[], void *user_data)
 {
   uint8_t *buf = munit_calloc(1024, sizeof(uint8_t));
+  assert_ptr_not_null(buf);
   sparsemap_t *map = (sparsemap_t *)test_api_setup(params, user_data);
 
   sparsemap_init(map, buf, 1024);
@@ -131,7 +227,8 @@ static void
 test_api_open_tear_down(void *fixture)
 {
   sparsemap_t *map = (sparsemap_t *)fixture;
-  free(map->m_data);
+  assert_ptr_not_null(map->m_data);
+  munit_free(map->m_data);
   test_api_tear_down(fixture);
 }
 static MunitResult
@@ -154,6 +251,7 @@ static void *
 test_api_set_data_size_setup(const MunitParameter params[], void *user_data)
 {
   uint8_t *buf = munit_calloc(1024, sizeof(uint8_t));
+  assert_ptr_not_null(buf);
   sparsemap_t *map = (sparsemap_t *)test_api_setup(params, user_data);
 
   sparsemap_init(map, buf, 1024);
@@ -165,7 +263,8 @@ static void
 test_api_set_data_size_tear_down(void *fixture)
 {
   sparsemap_t *map = (sparsemap_t *)fixture;
-  free(map->m_data);
+  assert_ptr_not_null(map->m_data);
+  munit_free(map->m_data);
   test_api_tear_down(fixture);
 }
 static MunitResult
@@ -187,6 +286,7 @@ static void *
 test_api_remaining_capacity_setup(const MunitParameter params[], void *user_data)
 {
   uint8_t *buf = munit_calloc(1024, sizeof(uint8_t));
+  assert_ptr_not_null(buf);
   sparsemap_t *map = (sparsemap_t *)test_api_setup(params, user_data);
 
   sparsemap_init(map, buf, 1024);
@@ -197,7 +297,8 @@ static void
 test_api_remaining_capacity_tear_down(void *fixture)
 {
   sparsemap_t *map = (sparsemap_t *)fixture;
-  free(map->m_data);
+  assert_ptr_not_null(map->m_data);
+  munit_free(map->m_data);
   test_api_tear_down(fixture);
 }
 static MunitResult
@@ -213,28 +314,32 @@ test_api_remaining_capacity(const MunitParameter params[], void *data)
   do {
     sparsemap_set(map, i++, true);
     cap = sparsemap_capacity_remaining(map);
-  } while (cap > 1.0);
-  // assert_true(i == 169985); when seed is 8675309
-  assert_true(cap <= 1.0);
+  } while (cap > 1.0 && errno != ENOSPC);
+  errno = 0;
+  assert_true(cap <= 2.0);
 
   sparsemap_clear(map);
+  cap = sparsemap_capacity_remaining(map);
+  assert_true(cap > 99);
+
   i = 0;
   do {
     int p = munit_rand_int_range(0, 150000);
     sparsemap_set(map, p, true);
     i++;
     cap = sparsemap_capacity_remaining(map);
-  } while (cap > 2.0);
-  // assert_true(i == 64); when seed is 8675309
+  } while (cap > 1.0 && errno != ENOSPC);
+  errno = 0;
   assert_true(cap <= 2.0);
 
   return MUNIT_OK;
 }
 
 static void *
-test_api_get_range_size_setup(const MunitParameter params[], void *user_data)
+test_api_get_capacity_setup(const MunitParameter params[], void *user_data)
 {
   uint8_t *buf = munit_calloc(1024, sizeof(uint8_t));
+  assert_ptr_not_null(buf);
   sparsemap_t *map = (sparsemap_t *)test_api_setup(params, user_data);
 
   sparsemap_init(map, buf, 1024);
@@ -243,14 +348,15 @@ test_api_get_range_size_setup(const MunitParameter params[], void *user_data)
   return (void *)map;
 }
 static void
-test_api_get_range_size_tear_down(void *fixture)
+test_api_get_capacity_tear_down(void *fixture)
 {
   sparsemap_t *map = (sparsemap_t *)fixture;
-  free(map->m_data);
+  assert_ptr_not_null(map->m_data);
+  munit_free(map->m_data);
   test_api_tear_down(fixture);
 }
 static MunitResult
-test_api_get_range_size(const MunitParameter params[], void *data)
+test_api_get_capacity(const MunitParameter params[], void *data)
 {
   sparsemap_t *map = (sparsemap_t *)data;
   (void)params;
@@ -259,8 +365,7 @@ test_api_get_range_size(const MunitParameter params[], void *data)
 
   sparsemap_set(map, 42, true);
   assert_true(sparsemap_is_set(map, 42));
-  size_t size = sparsemap_get_capacity(map);
-  assert_true(size == 1024);
+  assert_true(sparsemap_get_capacity(map) == 1024);
 
   return MUNIT_OK;
 }
@@ -269,6 +374,7 @@ static void *
 test_api_is_set_setup(const MunitParameter params[], void *user_data)
 {
   uint8_t *buf = munit_calloc(1024, sizeof(uint8_t));
+  assert_ptr_not_null(buf);
   sparsemap_t *map = (sparsemap_t *)test_api_setup(params, user_data);
 
   sparsemap_init(map, buf, 1024);
@@ -280,7 +386,8 @@ static void
 test_api_is_set_tear_down(void *fixture)
 {
   sparsemap_t *map = (sparsemap_t *)fixture;
-  free(map->m_data);
+  assert_ptr_not_null(map->m_data);
+  munit_free(map->m_data);
   test_api_tear_down(fixture);
 }
 static MunitResult
@@ -301,6 +408,7 @@ static void *
 test_api_set_setup(const MunitParameter params[], void *user_data)
 {
   uint8_t *buf = munit_calloc(1024, sizeof(uint8_t));
+  assert_ptr_not_null(buf);
   sparsemap_t *map = (sparsemap_t *)test_api_setup(params, user_data);
 
   sparsemap_init(map, buf, 1024);
@@ -311,7 +419,8 @@ static void
 test_api_set_tear_down(void *fixture)
 {
   sparsemap_t *map = (sparsemap_t *)fixture;
-  free(map->m_data);
+  assert_ptr_not_null(map->m_data);
+  munit_free(map->m_data);
   test_api_tear_down(fixture);
 }
 static MunitResult
@@ -337,9 +446,10 @@ test_api_set(const MunitParameter params[], void *data)
 }
 
 static void *
-test_api_get_start_offset_setup(const MunitParameter params[], void *user_data)
+test_api_get_starting_offset_setup(const MunitParameter params[], void *user_data)
 {
   uint8_t *buf = munit_calloc(1024, sizeof(uint8_t));
+  assert_ptr_not_null(buf);
   sparsemap_t *map = (sparsemap_t *)test_api_setup(params, user_data);
 
   sparsemap_init(map, buf, 1024);
@@ -348,14 +458,15 @@ test_api_get_start_offset_setup(const MunitParameter params[], void *user_data)
   return (void *)map;
 }
 static void
-test_api_get_start_offset_tear_down(void *fixture)
+test_api_get_starting_offset_tear_down(void *fixture)
 {
   sparsemap_t *map = (sparsemap_t *)fixture;
-  free(map->m_data);
+  assert_ptr_not_null(map->m_data);
+  munit_free(map->m_data);
   test_api_tear_down(fixture);
 }
 static MunitResult
-test_api_get_start_offset(const MunitParameter params[], void *data)
+test_api_get_starting_offset(const MunitParameter params[], void *data)
 {
   sparsemap_t *map = (sparsemap_t *)data;
   (void)params;
@@ -364,7 +475,7 @@ test_api_get_start_offset(const MunitParameter params[], void *data)
 
   sparsemap_set(map, 42, true);
   assert_true(sparsemap_is_set(map, 42));
-  size_t offset = sparsemap_get_start_offset(map);
+  size_t offset = sparsemap_get_starting_offset(map);
   assert_true(offset == 0);
 
   return MUNIT_OK;
@@ -374,6 +485,7 @@ static void *
 test_api_get_size_setup(const MunitParameter params[], void *user_data)
 {
   uint8_t *buf = munit_calloc(1024, sizeof(uint8_t));
+  assert_ptr_not_null(buf);
   sparsemap_t *map = (sparsemap_t *)test_api_setup(params, user_data);
 
   sparsemap_init(map, buf, 1024);
@@ -385,7 +497,8 @@ static void
 test_api_get_size_tear_down(void *fixture)
 {
   sparsemap_t *map = (sparsemap_t *)fixture;
-  free(map->m_data);
+  assert_ptr_not_null(map->m_data);
+  munit_free(map->m_data);
   test_api_tear_down(fixture);
 }
 static MunitResult
@@ -406,10 +519,11 @@ static void *
 test_api_scan_setup(const MunitParameter params[], void *user_data)
 {
   uint8_t *buf = munit_calloc(1024, sizeof(uint8_t));
+  assert_ptr_not_null(buf);
   sparsemap_t *map = (sparsemap_t *)test_api_setup(params, user_data);
 
   sparsemap_init(map, buf, 1024);
-  bitmap_from_uint64(map, ((uint64_t)0xfeedface << 32) | 0xbadc0ffee);
+  sm_bitmap_from_uint64(map, ((uint64_t)0xfeedface << 32) | 0xbadc0ffee);
 
   return (void *)map;
 }
@@ -417,7 +531,8 @@ static void
 test_api_scan_tear_down(void *fixture)
 {
   sparsemap_t *map = (sparsemap_t *)fixture;
-  free(map->m_data);
+  assert_ptr_not_null(map->m_data);
+  munit_free(map->m_data);
   test_api_tear_down(fixture);
 }
 void
@@ -436,7 +551,7 @@ test_api_scan(const MunitParameter params[], void *data)
   assert_ptr_not_null(map);
 
   sparsemap_set(map, 4200, true);
-  assert_true(sparsemap_is_set(map, 42));
+  assert_true(sparsemap_is_set(map, 4200));
   sparsemap_scan(map, scan_for_0xfeedfacebadcoffee, 0);
 
   return MUNIT_OK;
@@ -446,6 +561,7 @@ static void *
 test_api_split_setup(const MunitParameter params[], void *user_data)
 {
   uint8_t *buf = munit_calloc(1024, sizeof(uint8_t));
+  assert_ptr_not_null(buf);
   sparsemap_t *map = (sparsemap_t *)test_api_setup(params, user_data);
 
   sparsemap_init(map, buf, 1024);
@@ -458,7 +574,8 @@ static void
 test_api_split_tear_down(void *fixture)
 {
   sparsemap_t *map = (sparsemap_t *)fixture;
-  free(map->m_data);
+  assert_ptr_not_null(map->m_data);
+  munit_free(map->m_data);
   test_api_tear_down(fixture);
 }
 static MunitResult
@@ -489,10 +606,11 @@ static void *
 test_api_select_setup(const MunitParameter params[], void *user_data)
 {
   uint8_t *buf = munit_calloc(1024, sizeof(uint8_t));
+  assert_ptr_not_null(buf);
   sparsemap_t *map = (sparsemap_t *)test_api_setup(params, user_data);
 
   sparsemap_init(map, buf, 1024);
-  bitmap_from_uint64(map, ((uint64_t)0xfeedface << 32) | 0xbadc0ffee);
+  sm_bitmap_from_uint64(map, ((uint64_t)0xfeedface << 32) | 0xbadc0ffee);
 
   return (void *)map;
 }
@@ -500,7 +618,8 @@ static void
 test_api_select_tear_down(void *fixture)
 {
   sparsemap_t *map = (sparsemap_t *)fixture;
-  free(map->m_data);
+  assert_ptr_not_null(map->m_data);
+  munit_free(map->m_data);
   test_api_tear_down(fixture);
 }
 static MunitResult
@@ -513,17 +632,105 @@ test_api_select(const MunitParameter params[], void *data)
 
   /* NOTE: select() is 0-based, to get the bit position of the 1st logical bit set
      call select(map, 0), to get the 18th, select(map, 17), etc. */
-  assert_true(sparsemap_select(map, 0) == 1);
-  assert_true(sparsemap_select(map, 4) == 6);
-  assert_true(sparsemap_select(map, 17) == 26);
+  assert_true(sparsemap_select(map, 0, true) == 1);
+  assert_true(sparsemap_select(map, 4, true) == 6);
+  assert_true(sparsemap_select(map, 17, true) == 26);
 
   return MUNIT_OK;
 }
 
+#ifdef SELECT_FALSE
 static void *
-test_api_rank_setup(const MunitParameter params[], void *user_data)
+test_api_select_false_setup(const MunitParameter params[], void *user_data)
 {
   uint8_t *buf = munit_calloc(1024, sizeof(uint8_t));
+  assert_ptr_not_null(buf);
+  sparsemap_t *map = (sparsemap_t *)test_api_setup(params, user_data);
+
+  sparsemap_init(map, buf, 1024);
+  sm_bitmap_from_uint64(map, ((uint64_t)0xfeedface << 32) | 0xbadc0ffee);
+
+  return (void *)map;
+}
+static void
+test_api_select_false_tear_down(void *fixture)
+{
+  sparsemap_t *map = (sparsemap_t *)fixture;
+  assert_ptr_not_null(map->m_data);
+  munit_free(map->m_data);
+  test_api_tear_down(fixture);
+}
+static MunitResult
+test_api_select_false(const MunitParameter params[], void *data)
+{
+  sparsemap_t *map = (sparsemap_t *)data;
+  (void)params;
+  assert_ptr_not_null(map);
+
+  /* First few 0/off/unset-bits in ((uint64_t)0xfeedface << 32) | 0xbadc0ffee) expressed as an array of offsets. */
+  int off[] = { 0, 4, 16, 17, 18, 19, 20, 21, 25, 28, 30, 36, 37, 40, 42, 49, 52, 56, 64, 65 };
+  for (int i = 0; i < 20; i++) {
+    sparsemap_idx_t f = sparsemap_select(map, i, false);
+    assert_true(f == off[i]);
+    assert_true(sparsemap_is_set(map, f) == false);
+  }
+
+  return MUNIT_OK;
+}
+#endif
+
+#ifdef SELECT_NEG
+static void *
+test_api_select_neg_setup(const MunitParameter params[], void *user_data)
+{
+  uint8_t *buf = munit_calloc(1024, sizeof(uint8_t));
+  assert_ptr_not_null(buf);
+  sparsemap_t *map = (sparsemap_t *)test_api_setup(params, user_data);
+
+  sparsemap_init(map, buf, 1024);
+  sm_bitmap_from_uint64(map, ((uint64_t)0xfeedface << 32) | 0xbadc0ffee);
+
+  return (void *)map;
+}
+static void
+test_api_select_neg_tear_down(void *fixture)
+{
+  sparsemap_t *map = (sparsemap_t *)fixture;
+  assert_ptr_not_null(map->m_data);
+  munit_free(map->m_data);
+  test_api_tear_down(fixture);
+}
+static MunitResult
+test_api_select_neg(const MunitParameter params[], void *data)
+{
+  sparsemap_t *map = (sparsemap_t *)data;
+  (void)params;
+
+  assert_ptr_not_null(map);
+
+  sparsemap_set(map, 42, true);
+  sparsemap_set(map, 420, true);
+  sparsemap_set(map, 4200, true);
+
+  f = sparsemap_select(map, 0, false);
+  assert_true(f == 0);
+
+  f = sparsemap_select(map, -1, true);
+  assert_true(f == 4200);
+  f = sparsemap_select(map, -2, true);
+  assert_true(f == 420);
+  f = sparsemap_select(map, -3, true);
+  assert_true(f == 42);
+
+  return MUNIT_OK;
+}
+#endif
+
+static void *
+test_api_rank_true_setup(const MunitParameter params[], void *user_data)
+{
+  uint8_t *buf = munit_calloc(1024, sizeof(uint8_t));
+  assert_ptr_not_null(buf);
   sparsemap_t *map = (sparsemap_t *)test_api_setup(params, user_data);
 
   sparsemap_init(map, buf, 1024);
@@ -531,14 +738,15 @@ test_api_rank_setup(const MunitParameter params[], void *user_data)
   return (void *)map;
 }
 static void
-test_api_rank_tear_down(void *fixture)
+test_api_rank_true_tear_down(void *fixture)
 {
   sparsemap_t *map = (sparsemap_t *)fixture;
-  free(map->m_data);
+  assert_ptr_not_null(map->m_data);
+  munit_free(map->m_data);
   test_api_tear_down(fixture);
 }
 static MunitResult
-test_api_rank(const MunitParameter params[], void *data)
+test_api_rank_true(const MunitParameter params[], void *data)
 {
   int r1, r2;
   sparsemap_t *map = (sparsemap_t *)data;
@@ -549,36 +757,94 @@ test_api_rank(const MunitParameter params[], void *data)
   for (int i = 0; i < 10; i++) {
     sparsemap_set(map, i, true);
   }
-  for (int i = 0; i < 10; i++) {
-    assert_true(sparsemap_is_set(map, i));
-  }
-  for (int i = 10; i < 1000; i++) {
-    assert_true(!sparsemap_is_set(map, i));
-  }
   /* rank() is also 0-based, for consistency (and confusion sake); consider the
      range as [start, end] of [0, 9] counts the bits set in the first 10
      positions (starting from the LSB) in the index. */
   r1 = rank_uint64((uint64_t)-1, 0, 9);
-  r2 = sparsemap_rank(map, 0, 9);
+  r2 = sparsemap_rank(map, 0, 9, true);
   assert_true(r1 == r2);
-  assert_true(sparsemap_rank(map, 0, 9) == 10);
-  assert_true(sparsemap_rank(map, 1000, 1050) == 0);
+  assert_true(sparsemap_rank(map, 0, 9, true) == 10);
+  assert_true(sparsemap_rank(map, 1000, 1050, true) == 0);
 
-  for (int i = 0; i < 10; i++) {
-    for (int j = i; j < 10; j++) {
-      r1 = rank_uint64((uint64_t)-1, i, j);
-      r2 = sparsemap_rank(map, i, j);
-      assert_true(r1 == r2);
+  sparsemap_clear(map);
+
+  for (int i = 0; i < 10000; i++) {
+    sparsemap_set(map, i, true);
+  }
+  sparsemap_idx_t hole = 4999;
+  sparsemap_set(map, hole, false);
+  for (int i = 0; i < 10000; i++) {
+    for (int j = i; j < 10000; j++) {
+      int amt = (i > j) ? 0 : j - i + 1 - ((hole >= i && j >= hole) ? 1 : 0);
+      int r = sparsemap_rank(map, i, j, true);
+      assert_true(r == amt);
     }
   }
 
   return MUNIT_OK;
 }
 
+static void *
+test_api_rank_false_setup(const MunitParameter params[], void *user_data)
+{
+  uint8_t *buf = munit_calloc(1024, sizeof(uint8_t));
+  assert_ptr_not_null(buf);
+  sparsemap_t *map = (sparsemap_t *)test_api_setup(params, user_data);
+
+  sparsemap_init(map, buf, 1024);
+
+  return (void *)map;
+}
+static void
+test_api_rank_false_tear_down(void *fixture)
+{
+  sparsemap_t *map = (sparsemap_t *)fixture;
+  assert_ptr_not_null(map->m_data);
+  munit_free(map->m_data);
+  test_api_tear_down(fixture);
+}
+static MunitResult
+test_api_rank_false(const MunitParameter params[], void *data)
+{
+  int r;
+  sparsemap_t *map = (sparsemap_t *)data;
+  (void)params;
+
+  assert_ptr_not_null(map);
+
+  // empty map
+  for (int i = 0; i < 10000; i++) {
+    for (int j = i; j < 10000; j++) {
+      r = sparsemap_rank(map, i, j, false);
+      assert_true(r == j - i + 1);
+    }
+  }
+
+  // One chunk means not so empty now!
+  sparsemap_idx_t hole = 4999;
+  sparsemap_set(map, hole, true);
+  for (int i = 0; i < 10000; i++) {
+    for (int j = i; j < 10000; j++) {
+      int amt = (i > j) ? 0 : j - i + 1 - ((hole >= i && j >= hole) ? 1 : 0);
+      r = sparsemap_rank(map, i, j, false);
+      assert_true(r == amt);
+    }
+  }
+
+  sparsemap_clear(map);
+  sparsemap_set(map, 1, true);
+  sparsemap_set(map, 11, true);
+  r = sparsemap_rank(map, 0, 11, false);
+  assert_true(r == 10);
+
+  return MUNIT_OK;
+}
+
 static void *
 test_api_span_setup(const MunitParameter params[], void *user_data)
 {
   uint8_t *buf = munit_calloc(1024, sizeof(uint8_t));
+  assert_ptr_not_null(buf);
   sparsemap_t *map = (sparsemap_t *)test_api_setup(params, user_data);
 
   sparsemap_init(map, buf, 1024);
@@ -589,7 +855,8 @@ static void
 test_api_span_tear_down(void *fixture)
 {
   sparsemap_t *map = (sparsemap_t *)fixture;
-  free(map->m_data);
+  assert_ptr_not_null(map->m_data);
+  munit_free(map->m_data);
   test_api_tear_down(fixture);
 }
 static MunitResult
@@ -602,27 +869,27 @@ test_api_span(const MunitParameter params[], void *data)
 
   int located_at, placed_at, amt = 10000;
 
-  placed_at = create_sequential_set_in_empty_map(map, amt, 1);
-  located_at = sparsemap_span(map, 0, 1);
+  placed_at = sm_add_span(map, amt, 1);
+  located_at = sparsemap_span(map, 0, 1, true);
   assert_true(located_at == placed_at);
 
   sparsemap_clear(map);
 
-  placed_at = create_sequential_set_in_empty_map(map, amt, 50);
-  located_at = sparsemap_span(map, 0, 50);
+  placed_at = sm_add_span(map, amt, 50);
+  located_at = sparsemap_span(map, 0, 50, true);
   assert_true(located_at == placed_at);
 
   sparsemap_clear(map);
 
-  placed_at = create_sequential_set_in_empty_map(map, amt, 50);
-  located_at = sparsemap_span(map, placed_at / 2, 50);
+  placed_at = sm_add_span(map, amt, 50);
+  located_at = sparsemap_span(map, placed_at / 2, 50, true);
   assert_true(located_at == placed_at);
 
   /* TODO
   sparsemap_clear(map);
 
-  placed_at = create_sequential_set_in_empty_map(map, amt, amt - 1);
-  located_at = sparsemap_span(map, 0, amt - 1);
+  placed_at = sm_add_span(map, amt, amt - 1);
+  located_at = sparsemap_span(map, 0, amt - 1, true);
   assert_true(located_at == placed_at);
    */
 
@@ -631,31 +898,286 @@ test_api_span(const MunitParameter params[], void *data)
 
 // clang-format off
 static MunitTest api_test_suite[] = {
-  { (char *)"/static_init", test_api_static_init, NULL, NULL, MUNIT_TEST_OPTION_NONE, NULL },
+  { (char *)"/new", test_api_new, NULL, NULL, MUNIT_TEST_OPTION_NONE, NULL },
+  { (char *)"/new/realloc", test_api_new_realloc, NULL, NULL, MUNIT_TEST_OPTION_NONE, NULL },
+  { (char *)"/new/heap", test_api_new_heap, NULL, NULL, MUNIT_TEST_OPTION_NONE, NULL },
+  { (char *)"/new/static", test_api_new_static, NULL, NULL, MUNIT_TEST_OPTION_NONE, NULL },
+  { (char *)"/new/stack", test_api_new_stack, NULL, NULL, MUNIT_TEST_OPTION_NONE, NULL },
   { (char *)"/clear", test_api_clear, test_api_clear_setup, test_api_clear_tear_down, MUNIT_TEST_OPTION_NONE, NULL },
   { (char *)"/open", test_api_open, test_api_open_setup, test_api_open_tear_down, MUNIT_TEST_OPTION_NONE, NULL },
   { (char *)"/set_data_size", test_api_set_data_size, test_api_set_data_size_setup, test_api_set_data_size_tear_down, MUNIT_TEST_OPTION_NONE, NULL },
   { (char *)"/remaining_capacity", test_api_remaining_capacity, test_api_remaining_capacity_setup, test_api_remaining_capacity_tear_down,
     MUNIT_TEST_OPTION_NONE, NULL },
-  { (char *)"/get_range_size", test_api_get_range_size, test_api_get_range_size_setup, test_api_get_range_size_tear_down, MUNIT_TEST_OPTION_NONE, NULL },
+  { (char *)"/get_capacity", test_api_get_capacity, test_api_get_capacity_setup, test_api_get_capacity_tear_down, MUNIT_TEST_OPTION_NONE, NULL },
   { (char *)"/is_set", test_api_is_set, test_api_is_set_setup, test_api_is_set_tear_down, MUNIT_TEST_OPTION_NONE, NULL },
   { (char *)"/set", test_api_set, test_api_set_setup, test_api_set_tear_down, MUNIT_TEST_OPTION_NONE, NULL },
-  { (char *)"/get_start_offset", test_api_get_start_offset, test_api_get_start_offset_setup, test_api_get_start_offset_tear_down, MUNIT_TEST_OPTION_NONE,
-    NULL },
+  { (char *)"/get_starting_offset", test_api_get_starting_offset, test_api_get_starting_offset_setup, test_api_get_starting_offset_tear_down, MUNIT_TEST_OPTION_NONE, NULL },
+  //TODO { (char *)"/get_ending_offset", test_api_get_ending_offset, test_api_get_ending_offset_setup, test_api_get_ending_offset_tear_down, MUNIT_TEST_OPTION_NONE, NULL },
   { (char *)"/get_size", test_api_get_size, test_api_get_size_setup, test_api_get_size_tear_down, MUNIT_TEST_OPTION_NONE, NULL },
   { (char *)"/scan", test_api_scan, test_api_scan_setup, test_api_scan_tear_down, MUNIT_TEST_OPTION_NONE, NULL },
   { (char *)"/split", test_api_split, test_api_split_setup, test_api_split_tear_down, MUNIT_TEST_OPTION_NONE, NULL },
-  { (char *)"/select", test_api_select, test_api_select_setup, test_api_select_tear_down, MUNIT_TEST_OPTION_NONE, NULL },
-  { (char *)"/rank", test_api_rank, test_api_rank_setup, test_api_rank_tear_down, MUNIT_TEST_OPTION_NONE, NULL },
+  { (char *)"/select/true", test_api_select, test_api_select_setup, test_api_select_tear_down, MUNIT_TEST_OPTION_NONE, NULL },
+#ifdef SELECT_FALSE
+  { (char *)"/select/false", test_api_select_false, test_api_select_false_setup, test_api_select_false_tear_down, MUNIT_TEST_OPTION_NONE, NULL },
+#endif
+#ifdef SELECT_NEG
+  { (char *)"/select/neg", test_api_select_neg, test_api_select_neg_setup, test_api_select_neg_tear_down, MUNIT_TEST_OPTION_NONE, NULL },
+#endif
+  { (char *)"/rank/true", test_api_rank_true, test_api_rank_true_setup, test_api_rank_true_tear_down, MUNIT_TEST_OPTION_NONE, NULL },
+  { (char *)"/rank/false", test_api_rank_false, test_api_rank_false_setup, test_api_rank_false_tear_down, MUNIT_TEST_OPTION_NONE, NULL },
   { (char *)"/span", test_api_span, test_api_span_setup, test_api_span_tear_down, MUNIT_TEST_OPTION_NONE, NULL },
   { NULL, NULL, NULL, NULL, MUNIT_TEST_OPTION_NONE, NULL }
 };
 // clang-format on
 
+/* -------------------------- Scale Tests */
+
+static void *
+test_scale_lots_o_spans_setup(const MunitParameter params[], void *user_data)
+{
+  (void)params;
+  (void)user_data;
+  sparsemap_t *map = sparsemap(10 * 1024);
+  assert_ptr_not_null(map);
+  return (void *)map;
+}
+static void
+test_scale_lots_o_spans_tear_down(void *fixture)
+{
+  sparsemap_t *map = (sparsemap_t *)fixture;
+  assert_ptr_not_null(map);
+  munit_free(map);
+}
+static MunitResult
+test_scale_lots_o_spans(const MunitParameter params[], void *data)
+{
+  size_t amt = 897915; // 268435456
+  sparsemap_t *map = (sparsemap_t *)data;
+  (void)params;
+
+  assert_ptr_not_null(map);
+
+  for (size_t i = 0; i < amt;) {
+    int l = i % 31 + 16;
+    // TODO: sm_add_span(map, amt, l);
+    sm_add_span(map, 10000, l);
+    if (errno == ENOSPC) {
+      map = sparsemap_set_data_size(map, sparsemap_get_capacity(map) * 2);
+      errno = 0;
+    }
+    i += l;
+    /* ANSI esc code to clear line, carrage return, then print on the same line */
+    // printf("\033[2K\r%d", i);
+    // printf("%d\t%d\n", l, i);
+  }
+
+  return MUNIT_OK;
+}
+
+static void *
+test_scale_ondrej_setup(const MunitParameter params[], void *user_data)
+{
+  (void)params;
+  (void)user_data;
+  sparsemap_t *map = sparsemap(10 * 1024);
+  assert_ptr_not_null(map);
+  return (void *)map;
+}
+static void
+test_scale_ondrej_tear_down(void *fixture)
+{
+  sparsemap_t *map = (sparsemap_t *)fixture;
+  assert_ptr_not_null(map);
+  munit_free(map);
+}
+static MunitResult
+test_scale_ondrej(const MunitParameter params[], void *data)
+{
+  sparsemap_t *map = (sparsemap_t *)data;
+  (void)params;
+
+  assert_ptr_not_null(map);
+
+  sparsemap_idx_t stride = 18;
+  sparsemap_idx_t top = 268435456;
+  sparsemap_idx_t needle = munit_rand_int_range(1, top / stride);
+  for (sparsemap_idx_t i = 0; i < top / stride; i += stride) {
+    for (sparsemap_idx_t j = 0; j < stride; j++) {
+      bool set = (i != needle) ? (j < 10) : (j < 9);
+      sparsemap_set(map, i, set);
+      if (errno == ENOSPC) {
+        map = sparsemap_set_data_size(map, sparsemap_get_capacity(map) * 2);
+        errno = 0;
+      }
+    }
+    assert_true(sm_is_span(map, i + ((i != needle) ? 10 : 9), (i != needle) ? 8 : 9, false));
+  }
+  sparsemap_idx_t a = sparsemap_span(map, 0, 9, false);
+  sparsemap_idx_t l = a / stride;
+  printf("%ld\t%ld\n", a, l);
+  assert_true(l == needle);
+  return MUNIT_OK;
+}
+
+static void *
+test_scale_spans_come_spans_go_setup(const MunitParameter params[], void *user_data)
+{
+  (void)params;
+  (void)user_data;
+  sparsemap_t *map = sparsemap(1024);
+  assert_ptr_not_null(map);
+  return (void *)map;
+}
+static void
+test_scale_spans_come_spans_go_tear_down(void *fixture)
+{
+  sparsemap_t *map = (sparsemap_t *)fixture;
+  assert_ptr_not_null(map);
+  munit_free(map);
+}
+static MunitResult
+test_scale_spans_come_spans_go(const MunitParameter params[], void *data)
+{
+  size_t amt = 8192; // 268435456, ~5e7 interations due to 2e9 / avg(l)
+  sparsemap_t *map = (sparsemap_t *)data;
+  (void)params;
+
+  assert_ptr_not_null(map);
+
+  for (size_t i = 0; i < amt;) {
+    int l = i % 31 + 16;
+    sm_add_span(map, amt, l);
+    if (errno == ENOSPC) {
+      map = sparsemap_set_data_size(map, sparsemap_get_capacity(map) * 2);
+      errno = 0;
+    }
+
+    /* After 1,000 spans are in the map start consuming a span every iteration. */
+    if (l > 1000) {
+      do {
+        int s = munit_rand_int_range(1, 30);
+        int o = munit_rand_int_range(1, 268435456 - s - 1);
+        size_t b = sparsemap_span(map, o, s, true);
+        if (b == SPARSEMAP_IDX_MAX) {
+          continue;
+        }
+        for (int j = b; j < s; j++) {
+          assert_true(sparsemap_is_set(map, j) == true);
+        }
+        for (int j = b; j < s; j++) {
+          sparsemap_set(map, j, false);
+        }
+        for (int j = b; j < s; j++) {
+          assert_true(sparsemap_is_set(map, j) == false);
+        }
+        break;
+      } while (true);
+    }
+    i += l;
+  }
+
+  return MUNIT_OK;
+}
+
+static void *
+test_scale_best_case_setup(const MunitParameter params[], void *user_data)
+{
+  uint8_t *buf = munit_calloc(1024, sizeof(uint8_t));
+  assert_ptr_not_null(buf);
+  sparsemap_t *map = (sparsemap_t *)test_api_setup(params, user_data);
+
+  sparsemap_init(map, buf, 1024);
+
+  return (void *)map;
+}
+static void
+test_scale_best_case_tear_down(void *fixture)
+{
+  sparsemap_t *map = (sparsemap_t *)fixture;
+  assert_ptr_not_null(map->m_data);
+  munit_free(map->m_data);
+  test_api_tear_down(fixture);
+}
+static MunitResult
+test_scale_best_case(const MunitParameter params[], void *data)
+{
+  sparsemap_t *map = (sparsemap_t *)data;
+  (void)params;
+
+  assert_ptr_not_null(map);
+
+  /* Best case a map can contain 2048 bits in 8 bytes.
+     So, in a 1KiB buffer you have:
+       (1024 KiB / 8 bytes) * 2048 = 268,435,456 bits
+     or 1.09 TiB of 4KiB pages. Let's investigate, and find out if that's the case.
+
+     TODO: Actually, 172032 are stored before SEGV, or 706 MiB of 4KiB pages.
+  */
+
+  /* Set every bit on, that should be the best case. */
+  // for (int i = 0; i < 268435456; i++) {
+  for (int i = 0; i < 172032; i++) {
+    /* ANSI esc code to clear line, carrage return, then print on the same line */
+    //    printf("\033[2K\r%d", i);
+    sparsemap_set(map, i, true);
+  }
+
+  return MUNIT_OK;
+}
+
+static void *
+test_scale_worst_case_setup(const MunitParameter params[], void *user_data)
+{
+  uint8_t *buf = munit_calloc(1024, sizeof(uint8_t));
+  assert_ptr_not_null(buf);
+  sparsemap_t *map = (sparsemap_t *)test_api_setup(params, user_data);
+
+  sparsemap_init(map, buf, 1024);
+
+  return (void *)map;
+}
+static void
+test_scale_worst_case_tear_down(void *fixture)
+{
+  sparsemap_t *map = (sparsemap_t *)fixture;
+  assert_ptr_not_null(map->m_data);
+  munit_free(map->m_data);
+  test_api_tear_down(fixture);
+}
+static MunitResult
+test_scale_worst_case(const MunitParameter params[], void *data)
+{
+  sparsemap_t *map = (sparsemap_t *)data;
+  (void)params;
+
+  assert_ptr_not_null(map);
+
+  /* Worst case a map can contain 2048 bits in 265 + 8 = 264 bytes.
+     So, in a 1KiB buffer you have:
+       (1024 KiB / 264 bytes) * 2048 = 8,134,407.75758 bits
+     or 33.3 GiB of 4KiB pages. Let's investigate, and find out if that's the case.
+
+     TODO: actually 7744 are stored before SEGV, or 31MiB of 4KiB pages.
+  */
+
+  /* Set every other bit, that has to be the "worst case" for this index. */
+  // for (int i = 0; i < 8134407; i += 2) {
+  for (int i = 0; i < 7744; i += 2) {
+    /* ANSI esc code to clear line, carrage return, then print on the same line */
+    //    printf("\033[2K\r%d", i);
+    sparsemap_set(map, i, true);
+  }
+
+  return MUNIT_OK;
+}
+
+/* -------------------------- Performance Tests */
+
 static void *
 test_perf_span_solo_setup(const MunitParameter params[], void *user_data)
 {
   uint8_t *buf = munit_calloc(1024 * 3, sizeof(uint8_t));
+  assert_ptr_not_null(buf);
   sparsemap_t *map = (sparsemap_t *)test_api_setup(params, user_data);
 
   sparsemap_init(map, buf, 3 * 1024);
@@ -666,40 +1188,38 @@ static void
 test_perf_span_solo_tear_down(void *fixture)
 {
   sparsemap_t *map = (sparsemap_t *)fixture;
-  free(map->m_data);
+  assert_ptr_not_null(map->m_data);
+  munit_free(map->m_data);
   test_api_tear_down(fixture);
 }
-EST_MEDIAN_DECL(solo, 10000)
 static MunitResult
 test_perf_span_solo(const MunitParameter params[], void *data)
 {
   sparsemap_t *map = (sparsemap_t *)data;
-  uint64_t stop, start;
+  // double stop, start;
   (void)params;
-  int located_at, placed_at, amt = 1000;
+  int located_at, placed_at, amt = 500;
 
   assert_ptr_not_null(map);
+  return MUNIT_OK; // TODO
 
   for (int i = 1; i < amt; i++) {
-    for (int j = 1; j < amt / 10; j++) {
+    for (int length = 1; length <= 100; length++) {
       sparsemap_clear(map);
-      placed_at = create_sequential_set_in_empty_map(map, amt, j);
-      // logf("i = %d, j = %d\tplaced_at %d\n", i, j, placed_at);
-      // whats_set(map, 5000);
-      start = tsc();
-      located_at = sparsemap_span(map, 0, j);
-      stop = tsc();
-      //      fprintf(stdout, "%ll - %ll = %ll\n", stop, start, stop - start);
-      EST_MEDIAN_ADD(solo, stop - start);
+      placed_at = sm_add_span(map, amt, length);
+      // logf("i = %d, length = %d\tplaced_at %d\n", i, length, placed_at);
+      // sm_whats_set(map, 5000);
+      // start = nsts();
+      located_at = sparsemap_span(map, 0, length, true);
+      // stop = nsts();
+      // double amt = (stop - start) * 1e6;
+      // if (amt > 0) {
+      //   fprintf(stdout, "%0.8f\n", amt);
+      // }
       if (placed_at != located_at)
-        logf("a: i = %d, j = %d\tplaced_at %d located_at %d\n", i, j, placed_at, located_at);
+        logf("a: i = %d, length = %d\tplaced_at %d located_at %d\n", i, length, placed_at, located_at);
     }
   }
-  uint64_t est = EST_MEDIAN_GET(solo);
-  // fprintf(stdout, "median time %zu or %f ns\n", est, tsc_ticks_to_ns(est)); // measured 228
-  assert_true(est < 500);
-  fflush(stdout);
-
   return MUNIT_OK;
 }
 
@@ -707,6 +1227,7 @@ static void *
 test_perf_span_tainted_setup(const MunitParameter params[], void *user_data)
 {
   uint8_t *buf = munit_calloc(1024 * 3, sizeof(uint8_t));
+  assert_ptr_not_null(buf);
   sparsemap_t *map = (sparsemap_t *)test_api_setup(params, user_data);
 
   sparsemap_init(map, buf, 3 * 1024);
@@ -717,48 +1238,52 @@ static void
 test_perf_span_tainted_tear_down(void *fixture)
 {
   sparsemap_t *map = (sparsemap_t *)fixture;
-  free(map->m_data);
+  assert_ptr_not_null(map->m_data);
+  munit_free(map->m_data);
   test_api_tear_down(fixture);
 }
-EST_MEDIAN_DECL(tainted, 10000)
 static MunitResult
 test_perf_span_tainted(const MunitParameter params[], void *data)
 {
   sparsemap_t *map = (sparsemap_t *)data;
-  uint64_t stop, start;
+  // double stop, start;
   (void)params;
 
   assert_ptr_not_null(map);
 
-  int located_at, placed_at, amt = 1000;
+  int located_at, placed_at, amt = 500;
   for (int i = 1; i < amt; i++) {
-    for (int j = 1; j < amt / 10; j++) {
+    for (int j = 100; j <= 10; j++) {
       sparsemap_clear(map);
       populate_map(map, 1024, 1 * 1024);
-      placed_at = create_sequential_set_in_empty_map(map, amt, j);
-      start = tsc();
-      located_at = sparsemap_span(map, 0, j);
-      stop = tsc();
-      EST_MEDIAN_ADD(tainted, stop - start);
+      placed_at = sm_add_span(map, amt, j);
+      // start = nsts();
+      located_at = sparsemap_span(map, 0, j, true);
+      // stop = nsts();
+      // double amt = (stop - start) * 1e6;
+      // if (amt > 0) {
+      // fprintf(stdout, "%0.8f\n", amt);
+      // }
       if (located_at >= placed_at)
         logf("b: i = %d, j = %d\tplaced_at %d located_at %d\n", i, j, placed_at, located_at);
-      // assert_true(located_at >= placed_at);
-      // start = tsc();
-      // located_at = sparsemap_span(map, (placed_at < j ? 0 : placed_at / 2), i);
-      // stop = tsc();
-      // EST_MEDIAN_ADD(solo, stop - start);
-      // assert_true(placed_at == located_at);
     }
   }
-  uint64_t est = EST_MEDIAN_GET(tainted);
-  // fprintf(stdout, "median time %zu or %f ns\n", est, tsc_ticks_to_ns(est)); // measured 228
-  assert_true(est < 500);
 
   return MUNIT_OK;
 }
 
 // clang-format off
-static MunitTest performance_test_suite[] = {
+static MunitTest scale_test_suite[] = {
+  { (char *)"/lots-o-spans", test_scale_lots_o_spans, test_scale_lots_o_spans_setup, test_scale_lots_o_spans_tear_down, MUNIT_TEST_OPTION_NONE, NULL },
+  { (char *)"/ondrej", test_scale_ondrej, test_scale_ondrej_setup, test_scale_ondrej_tear_down, MUNIT_TEST_OPTION_NONE, NULL },
+  { (char *)"/spans_come_spans_go", test_scale_spans_come_spans_go, test_scale_spans_come_spans_go_setup, test_scale_spans_come_spans_go_tear_down, MUNIT_TEST_OPTION_NONE, NULL },
+  { (char *)"/best-case", test_scale_best_case, test_scale_best_case_setup, test_scale_best_case_tear_down, MUNIT_TEST_OPTION_NONE, NULL },
+  { (char *)"/worst-case", test_scale_worst_case, test_scale_worst_case_setup, test_scale_worst_case_tear_down, MUNIT_TEST_OPTION_NONE, NULL },
+  { NULL, NULL, NULL, NULL, MUNIT_TEST_OPTION_NONE, NULL } };
+// clang-format on
+
+// clang-format off
+static MunitTest perf_test_suite[] = {
   { (char *)"/span/solo", test_perf_span_solo, test_perf_span_solo_setup, test_perf_span_solo_tear_down, MUNIT_TEST_OPTION_NONE, NULL },
   { (char *)"/span/tainted", test_perf_span_tainted, test_perf_span_tainted_setup, test_perf_span_tainted_tear_down, MUNIT_TEST_OPTION_NONE, NULL },
   { NULL, NULL, NULL, NULL, MUNIT_TEST_OPTION_NONE, NULL } };
@@ -766,16 +1291,29 @@ static MunitTest performance_test_suite[] = {
 
 // clang-format off
 static MunitSuite other_test_suite[] = {
-  { "/performance", performance_test_suite, NULL, 1, MUNIT_SUITE_OPTION_NONE },
+  { "/api", api_test_suite, NULL, 1, MUNIT_SUITE_OPTION_NONE },
+  { "/perf", perf_test_suite, NULL, 1, MUNIT_SUITE_OPTION_NONE },
+  { "/scale", scale_test_suite, NULL, 1, MUNIT_SUITE_OPTION_NONE },
   { NULL, NULL, NULL, 0, MUNIT_SUITE_OPTION_NONE } };
 // clang-format on
 
-static const MunitSuite main_test_suite = { (char *)"/api", api_test_suite, other_test_suite, 1, MUNIT_SUITE_OPTION_NONE };
+// clang-format off
+static MunitTest sparsemap_test_suite[] = {
+  { NULL, NULL, NULL, NULL, MUNIT_TEST_OPTION_NONE, NULL }
+};
+// clang-format on
+
+static const MunitSuite main_test_suite = { (char *)"/sparsemap", sparsemap_test_suite, other_test_suite, 1, MUNIT_SUITE_OPTION_NONE };
 
 int
 main(int argc, char *argv[MUNIT_ARRAY_PARAM(argc + 1)])
 {
   struct user_data info;
+
+  /* Disable buffering on std{out,err} to avoid having to call fflush(). */
+  setvbuf(stdout, NULL, _IONBF, 0);
+  setvbuf(stderr, NULL, _IONBF, 0);
+
   return munit_suite_main(&main_test_suite, (void *)&info, argc, argv);
 }