Negate sparse arrays into run containers; SIMD in-place array union

include/roaring/containers/containers.h

@@ -2018,9 +2018,8 @@ static inline container_t *container_not(const container_t *c, uint8_t type,
                     : ARRAY_CONTAINER_TYPE;
             return result;
         case ARRAY_CONTAINER_TYPE:
-            result = bitset_container_create();
-            *result_type = BITSET_CONTAINER_TYPE;
-            array_container_negation(const_CAST_array(c), CAST_bitset(result));
+            *result_type = (uint8_t)array_container_negation_range_to_any(
+                const_CAST_array(c), 0, 1 << 16, &result);
             return result;
         case RUN_CONTAINER_TYPE:
             *result_type =
@@ -2052,11 +2051,8 @@ static inline container_t *container_not_range(const container_t *c,
                     : ARRAY_CONTAINER_TYPE;
             return result;
         case ARRAY_CONTAINER_TYPE:
-            *result_type =
-                array_container_negation_range(const_CAST_array(c), range_start,
-                                               range_end, &result)
-                    ? BITSET_CONTAINER_TYPE
-                    : ARRAY_CONTAINER_TYPE;
+            *result_type = (uint8_t)array_container_negation_range_to_any(
+                const_CAST_array(c), range_start, range_end, &result);
             return result;
         case RUN_CONTAINER_TYPE:
             *result_type = (uint8_t)run_container_negation_range(
@@ -2085,9 +2081,8 @@ static inline container_t *container_inot(container_t *c, uint8_t type,
             return result;
         case ARRAY_CONTAINER_TYPE:
             // will never be inplace
-            result = bitset_container_create();
-            *result_type = BITSET_CONTAINER_TYPE;
-            array_container_negation(CAST_array(c), CAST_bitset(result));
+            *result_type = (uint8_t)array_container_negation_range_to_any(
+                CAST_array(c), 0, 1 << 16, &result);
             array_container_free(CAST_array(c));
             return result;
         case RUN_CONTAINER_TYPE:
@@ -2118,10 +2113,9 @@ static inline container_t *container_inot_range(container_t *c, uint8_t type,
                                : ARRAY_CONTAINER_TYPE;
             return result;
         case ARRAY_CONTAINER_TYPE:
-            *result_type = array_container_negation_range_inplace(
-                               CAST_array(c), range_start, range_end, &result)
-                               ? BITSET_CONTAINER_TYPE
-                               : ARRAY_CONTAINER_TYPE;
+            *result_type = (uint8_t)array_container_negation_range_to_any(
+                CAST_array(c), range_start, range_end, &result);
+            array_container_free(CAST_array(c));
             return result;
         case RUN_CONTAINER_TYPE:
             *result_type = (uint8_t)run_container_negation_range_inplace(

include/roaring/containers/mixed_negation.h

@@ -93,6 +93,17 @@ bool array_container_negation_range_inplace(array_container_t *src,
                                             const int range_end,
                                             container_t **dst);
 
+/* Negation across a range of the container.
+ * Compute the negation of src and write the result to *dst, which is not
+ * preallocated. Unlike array_container_negation_range, the result may be a
+ * run container: the negation of a sparse array over a large range has few
+ * runs. Returns the typecode of the result.
+ */
+int array_container_negation_range_to_any(const array_container_t *src,
+                                          const int range_start,
+                                          const int range_end,
+                                          container_t **dst);
+
 /* Negation across a range of the container
  * Compute the  negation of src  and write the result
  * to *dst.  A true return value indicates a bitset result,

src/containers/mixed_negation.c

@@ -172,6 +172,71 @@ bool array_container_negation_range_inplace(array_container_t *src,
     return ans;
 }
 
+// Appends the run [start, end] (inclusive) after the n runs, merging it with
+// the last one when they touch. Returns the new number of runs.
+static inline int32_t append_run(rle16_t *runs, int32_t n, uint32_t start,
+                                 uint32_t end) {
+    if (n > 0 &&
+        (uint32_t)runs[n - 1].value + runs[n - 1].length + 1 == start) {
+        runs[n - 1].length = (uint16_t)(end - runs[n - 1].value);
+        return n;
+    }
+    runs[n].value = (uint16_t)start;
+    runs[n].length = (uint16_t)(end - start);
+    return n + 1;
+}
+
+int array_container_negation_range_to_any(const array_container_t *src,
+                                          const int range_start,
+                                          const int range_end,
+                                          container_t **dst) {
+    // The result has more than DEFAULT_MAX_SIZE values, so it would be a
+    // bitset, but it has at most cardinality + 1 runs: when these runs take
+    // less room than a bitset, we build them directly. In the range, the
+    // runs are the gaps between the values; outside, the values themselves.
+    if (range_end - range_start - src->cardinality > DEFAULT_MAX_SIZE &&
+        run_container_serialized_size_in_bytes(src->cardinality + 1) <
+            BITSET_CONTAINER_SIZE_IN_WORDS * (int32_t)sizeof(uint64_t)) {
+        run_container_t *ans =
+            run_container_create_given_capacity(src->cardinality + 1);
+        if (ans == NULL) {
+            *dst = NULL;
+            return RUN_CONTAINER_TYPE;
+        }
+        const uint16_t *array = src->array;
+        const int32_t card = src->cardinality;
+        rle16_t *runs = ans->runs;
+        int32_t n = 0, k = 0;
+        for (; k < card && array[k] < range_start; k++) {
+            n = append_run(runs, n, array[k], array[k]);
+        }
+        uint32_t next = (uint32_t)range_start;  // first value not yet decided
+        for (; k < card && array[k] < range_end; k++) {
+            if (array[k] > next) n = append_run(runs, n, next, array[k] - 1u);
+            next = array[k] + 1u;
+        }
+        if (next < (uint32_t)range_end) {
+            n = append_run(runs, n, next, (uint32_t)range_end - 1);
+        }
+        for (; k < card; k++) {
+            n = append_run(runs, n, array[k], array[k]);
+        }
+        ans->n_runs = n;
+        *dst = ans;
+        return RUN_CONTAINER_TYPE;
+    }
+    if (range_start == 0 && range_end == (1 << 16)) {
+        // The negation of a whole array is always a bitset, built in one pass.
+        bitset_container_t *ans = bitset_container_create();
+        if (ans != NULL) array_container_negation(src, ans);
+        *dst = ans;
+        return BITSET_CONTAINER_TYPE;
+    }
+    return array_container_negation_range(src, range_start, range_end, dst)
+               ? BITSET_CONTAINER_TYPE
+               : ARRAY_CONTAINER_TYPE;
+}
+
 /* Negation across a range of the container
  * Compute the  negation of src  and write the result
  * to *dst.  A true return value indicates a bitset result,

src/containers/mixed_union.c

@@ -206,15 +206,17 @@ bool array_array_container_inplace_union(array_container_t *src_1,
             }
             return false;  // not a bitset
         } else {
-            memmove(src_1->array + src_2->cardinality, src_1->array,
-                    src_1->cardinality * sizeof(uint16_t));
-            // In theory, we could use fast_union_uint16, but it is unsafe. It
-            // fails with Intel compilers in particular.
+            // fast_union_uint16 cannot work in place: its vectorized stores
+            // could overwrite values not yet read.
             // https://github.com/RoaringBitmap/CRoaring/pull/452
             // See report https://github.com/RoaringBitmap/CRoaring/issues/476
-            src_1->cardinality = (int32_t)union_uint16(
-                src_1->array + src_2->cardinality, src_1->cardinality,
-                src_2->array, src_2->cardinality, src_1->array);
+            // So we copy src_1 aside first, which costs about as much as
+            // moving it to the end of its array for the scalar union.
+            uint16_t buffer[DEFAULT_MAX_SIZE];
+            memcpy(buffer, src_1->array, src_1->cardinality * sizeof(uint16_t));
+            src_1->cardinality = (int32_t)fast_union_uint16(
+                buffer, src_1->cardinality, src_2->array, src_2->cardinality,
+                src_1->array);
             return false;  // not a bitset
         }
     }

tests/mixed_container_unit.c

@@ -1991,11 +1991,59 @@ DEFINE_TEST(run_negation_range_test9) {
                              false, false);
 }
 
+// In-place union of two arrays when the first one has room for the result,
+// compared with the union into a new container.
+DEFINE_TEST(array_array_inplace_union_spare_capacity_test) {
+    uint32_t seed = 12345;
+    for (int trial = 0; trial < 300; trial++) {
+        int n1 = trial % 41 * 50, n2 = (trial * 7) % 37 * 50;
+        array_container_t* a = array_container_create_given_capacity(4096);
+        array_container_t* b = array_container_create();
+        assert_non_null(a);
+        assert_non_null(b);
+        for (int i = 0; i < n1; i++) {
+            seed = seed * 1664525 + 1013904223;
+            array_container_add(
+                a, (uint16_t)(seed >> 16) % (trial % 3 ? 65535 : 5000));
+        }
+        for (int i = 0; i < n2; i++) {
+            seed = seed * 1664525 + 1013904223;
+            array_container_add(
+                b, (uint16_t)(seed >> 16) % (trial % 3 ? 65535 : 5000));
+        }
+        array_container_t* want = array_container_create();
+        assert_non_null(want);
+        array_container_union(a, b, want);
+        container_t* dst = NULL;
+        bool is_bitset = array_array_container_inplace_union(a, b, &dst);
+        if (a->cardinality + b->cardinality <= DEFAULT_MAX_SIZE &&
+            dst == NULL) {
+            assert_false(is_bitset);
+            assert_true(array_container_equals(a, want));
+        } else {
+            uint8_t type =
+                is_bitset ? BITSET_CONTAINER_TYPE : ARRAY_CONTAINER_TYPE;
+            container_t* result = dst != NULL ? dst : (container_t*)a;
+            int card = container_get_cardinality(result, type);
+            assert_int_equal(card, want->cardinality);
+            for (int i = 0; i < want->cardinality; i++) {
+                assert_true(container_contains(result, want->array[i], type));
+            }
+            if (dst != NULL && dst != (container_t*)a)
+                container_free(dst, type);
+        }
+        array_container_free(a);
+        array_container_free(b);
+        array_container_free(want);
+    }
+}
+
 int main() {
     tellmeall();
 
     const struct CMUnitTest tests[] = {
         cmocka_unit_test(array_bitset_and_or_xor_andnot_test),
+        cmocka_unit_test(array_array_inplace_union_spare_capacity_test),
         cmocka_unit_test(array_bitset_run_lazy_xor_test),
         cmocka_unit_test(run_xor_test),
         cmocka_unit_test(run_ixor_test),

tests/toplevel_unit.c

@@ -3545,6 +3545,82 @@ DEFINE_TEST(test_flip_run_container_removal2) {
     roaring_bitmap_free(bm);
 }
 
+// Flipping a sparse array container over a large range gives few runs.
+DEFINE_TEST(test_flip_sparse_array_container) {
+    const uint32_t cards[] = {1, 2, 100, 2045, 2046, 2047, 3000, 4096};
+    const uint32_t ranges[][2] = {
+        {0, 65536}, {0, 1 << 20}, {100, 60000}, {7, 5000}, {65530, 131072}};
+    for (size_t c = 0; c < sizeof(cards) / sizeof(cards[0]); c++) {
+        for (size_t g = 0; g < sizeof(ranges) / sizeof(ranges[0]); g++) {
+            roaring_bitmap_t *bm = roaring_bitmap_create();
+            // Spread the values over the first chunk, with some neighbors.
+            for (uint32_t k = 0; k < cards[c]; k++) {
+                roaring_bitmap_add(bm, (k * 15u) % 65536 + (k & 1));
+            }
+            uint32_t start = ranges[g][0], end = ranges[g][1];
+            roaring_bitmap_t *flipped = roaring_bitmap_flip(bm, start, end);
+            roaring_bitmap_t *inplace = roaring_bitmap_copy(bm);
+            roaring_bitmap_flip_inplace(inplace, start, end);
+            assert_true(roaring_bitmap_internal_validate(flipped, NULL));
+            assert_true(roaring_bitmap_internal_validate(inplace, NULL));
+            assert_true(roaring_bitmap_equals(flipped, inplace));
+            uint64_t expected =
+                roaring_bitmap_get_cardinality(bm) + end - start -
+                2 * roaring_bitmap_range_cardinality(bm, start, end);
+            assert_int_equal(roaring_bitmap_get_cardinality(flipped), expected);
+            for (uint32_t v = 0; v < 70000; v++) {
+                bool in_range = v >= start && v < end;
+                assert_true(roaring_bitmap_contains(flipped, v) ==
+                            (roaring_bitmap_contains(bm, v) != in_range));
+            }
+            if (cards[c] <= 2046 && start == 0 && end % 65536 == 0) {
+                // The flipped chunks are all runs, whole or partial flips.
+                roaring_statistics_t stats;
+                roaring_bitmap_statistics(flipped, &stats);
+                assert_int_equal(stats.n_run_containers, end / 65536);
+                roaring_bitmap_statistics(inplace, &stats);
+                assert_int_equal(stats.n_run_containers, end / 65536);
+            }
+            roaring_bitmap_free(bm);
+            roaring_bitmap_free(flipped);
+            roaring_bitmap_free(inplace);
+        }
+    }
+    // Clustered values around the ends of the range, so that runs touch.
+    uint32_t seed = 4321;
+    for (int trial = 0; trial < 300; trial++) {
+        seed = seed * 1664525 + 1013904223;
+        uint32_t start = (seed >> 8) % 3000;
+        seed = seed * 1664525 + 1013904223;
+        uint32_t end = 65536 - (seed >> 8) % 3000;
+        if (trial % 3 == 0) end = 65536;
+        if (trial % 5 == 0) start = 0;
+        roaring_bitmap_t *bm = roaring_bitmap_create();
+        const uint32_t centers[] = {0, start, end, 30000, 65535};
+        for (size_t c = 0; c < 5; c++) {
+            for (int k = 0; k < 40; k++) {
+                seed = seed * 1664525 + 1013904223;
+                int64_t v =
+                    (int64_t)centers[c] + (int64_t)((seed >> 8) % 9) - 4;
+                if (v >= 0 && v < 65536) roaring_bitmap_add(bm, (uint32_t)v);
+            }
+        }
+        roaring_bitmap_t *flipped = roaring_bitmap_flip(bm, start, end);
+        roaring_bitmap_t *inplace = roaring_bitmap_copy(bm);
+        roaring_bitmap_flip_inplace(inplace, start, end);
+        assert_true(roaring_bitmap_internal_validate(flipped, NULL));
+        assert_true(roaring_bitmap_equals(flipped, inplace));
+        for (uint32_t v = 0; v < 65536; v++) {
+            bool in_range = v >= start && v < end;
+            assert_true(roaring_bitmap_contains(flipped, v) ==
+                        (roaring_bitmap_contains(bm, v) != in_range));
+        }
+        roaring_bitmap_free(bm);
+        roaring_bitmap_free(flipped);
+        roaring_bitmap_free(inplace);
+    }
+}
+
 // randomized test for rank query
 DEFINE_TEST(select_test) {
     const int min_runs = 1;
@@ -5863,6 +5939,7 @@ int main() {
         cmocka_unit_test(test_flip_bitset_container_removal),
         cmocka_unit_test(test_flip_run_container_removal),
         cmocka_unit_test(test_flip_run_container_removal2),
+        cmocka_unit_test(test_flip_sparse_array_container),
         cmocka_unit_test(select_test),
         cmocka_unit_test(test_subset),
         cmocka_unit_test(test_or_many_memory_leak),