Branch data Line data Source code
1 : : /**
2 : : * @file v2_bitvec.c
3 : : * @brief Ranked bit vector and LOUDS navigation — format v2 §5 and §6.
4 : : *
5 : : * Copyright (c) 2026 M. A. Chatterjee <deftio at deftio dot com>
6 : : * BSD-2-Clause — see LICENSE.txt
7 : : */
8 : :
9 : : #include "triepack/v2_bitvec.h"
10 : :
11 : : #include <stdlib.h>
12 : : #include <string.h>
13 : :
14 : : /* §5.1 layout constants. Conventional rather than measured; see the open
15 : : questions in the format spec. */
16 : : #define TP2_SUPER_BITS 2048
17 : : #define TP2_BLOCK_BITS 256
18 : :
19 : : struct tp2_bitvec {
20 : : uint8_t *bits;
21 : : uint64_t nbits;
22 : : uint64_t *super; /* rank1 of everything before this superblock */
23 : : uint16_t *block; /* rank1 within the superblock, before this block */
24 : : uint64_t nsuper, nblock;
25 : : };
26 : :
27 : 502844 : static int popcount64(uint64_t x)
28 : : {
29 : : #if defined(__GNUC__) || defined(__clang__)
30 : 502844 : return __builtin_popcountll(x);
31 : : #else
32 : : x = x - ((x >> 1) & 0x5555555555555555ull);
33 : : x = (x & 0x3333333333333333ull) + ((x >> 2) & 0x3333333333333333ull);
34 : : x = (x + (x >> 4)) & 0x0F0F0F0F0F0F0F0Full;
35 : : return (int)((x * 0x0101010101010101ull) >> 56);
36 : : #endif
37 : : }
38 : :
39 : 43 : tp2_bitvec *tp2_bv_create(const uint8_t *bits, uint64_t nbits)
40 : : {
41 : 43 : tp2_bitvec *bv = calloc(1, sizeof(*bv));
42 [ - + ]: 43 : if (!bv)
43 : 0 : return NULL;
44 : :
45 : 43 : bv->nbits = nbits;
46 : 43 : size_t nbytes = (size_t)((nbits + 7) / 8);
47 : : /* One spare byte so a word-at-a-time read near the end stays in bounds. */
48 : 43 : bv->bits = calloc(nbytes + 8, 1);
49 [ - + ]: 43 : if (!bv->bits) {
50 : 0 : free(bv);
51 : 0 : return NULL;
52 : : }
53 [ + + + - ]: 43 : if (bits && nbytes)
54 : 42 : memcpy(bv->bits, bits, nbytes);
55 : :
56 : : /* Trailing bits of the final byte are not part of the vector and must
57 : : not be counted; zero them so rank never sees them. */
58 [ + + ]: 43 : if (nbits % 8)
59 : 27 : bv->bits[nbytes - 1] &= (uint8_t)(0xFF << (8 - (nbits % 8)));
60 : :
61 : 43 : bv->nsuper = (nbits + TP2_SUPER_BITS - 1) / TP2_SUPER_BITS + 1;
62 : 43 : bv->nblock = (nbits + TP2_BLOCK_BITS - 1) / TP2_BLOCK_BITS + 1;
63 : 43 : bv->super = calloc((size_t)bv->nsuper, sizeof(*bv->super));
64 : 43 : bv->block = calloc((size_t)bv->nblock, sizeof(*bv->block));
65 [ + - - + ]: 43 : if (!bv->super || !bv->block) {
66 : 0 : tp2_bv_destroy(&bv);
67 : 0 : return NULL;
68 : : }
69 : :
70 : : /* Run to <= nbits so the final entries are produced by the same rule as
71 : : every other one. Writing them afterwards as separate "terminators"
72 : : gives them different semantics, and where a superblock boundary lands
73 : : exactly at the end both entries hold the same count and rank adds it
74 : : twice -- at n = 2048, rank1(2048) returned 1366 instead of 683. */
75 : 43 : uint64_t total = 0, within = 0;
76 [ + + ]: 37893 : for (uint64_t i = 0; i <= nbits; i++) {
77 [ + + ]: 37850 : if (i % TP2_SUPER_BITS == 0) {
78 : 57 : bv->super[i / TP2_SUPER_BITS] = total;
79 : 57 : within = 0;
80 : : }
81 [ + + ]: 37850 : if (i % TP2_BLOCK_BITS == 0)
82 : 183 : bv->block[i / TP2_BLOCK_BITS] = (uint16_t)within;
83 [ + + + + ]: 37850 : if (i < nbits && tp2_bv_get(bv, i)) {
84 : 22890 : total++;
85 : 22890 : within++;
86 : : }
87 : : }
88 : 43 : return bv;
89 : : }
90 : :
91 : 45 : void tp2_bv_destroy(tp2_bitvec **bv)
92 : : {
93 [ + + + + ]: 45 : if (!bv || !*bv)
94 : 2 : return;
95 : 43 : free((*bv)->bits);
96 : 43 : free((*bv)->super);
97 : 43 : free((*bv)->block);
98 : 43 : free(*bv);
99 : 43 : *bv = NULL;
100 : : }
101 : :
102 : 6 : uint64_t tp2_bv_size(const tp2_bitvec *bv)
103 : : {
104 [ + - ]: 6 : return bv ? bv->nbits : 0;
105 : : }
106 : :
107 : 394145 : int tp2_bv_get(const tp2_bitvec *bv, uint64_t i)
108 : : {
109 : 394145 : return (bv->bits[i >> 3] >> (7 - (i & 7))) & 1;
110 : : }
111 : :
112 : 103772 : uint64_t tp2_bv_rank1(const tp2_bitvec *bv, uint64_t i)
113 : : {
114 [ + - + + ]: 103772 : if (!bv || i == 0)
115 : 39 : return 0;
116 [ - + ]: 103733 : if (i > bv->nbits)
117 : 0 : i = bv->nbits;
118 : :
119 : 103733 : uint64_t c = bv->super[i / TP2_SUPER_BITS] + bv->block[i / TP2_BLOCK_BITS];
120 : 103733 : uint64_t j = (i / TP2_BLOCK_BITS) * TP2_BLOCK_BITS;
121 : :
122 : : /* Whole 64-bit words, then whole bytes, then the ragged tail. */
123 [ + + ]: 268057 : for (; j + 64 <= i; j += 64) {
124 : 164324 : uint64_t w = 0;
125 [ + + ]: 1478916 : for (int k = 0; k < 8; k++)
126 : 1314592 : w = (w << 8) | bv->bits[(j >> 3) + (uint64_t)k];
127 : 164324 : c += (uint64_t)popcount64(w);
128 : : }
129 [ + + ]: 442253 : for (; j + 8 <= i; j += 8)
130 : 338520 : c += (uint64_t)popcount64((uint64_t)bv->bits[j >> 3]);
131 [ + + ]: 460044 : for (; j < i; j++)
132 : 356311 : c += (uint64_t)tp2_bv_get(bv, j);
133 : 103733 : return c;
134 : : }
135 : :
136 : 56373 : uint64_t tp2_bv_rank0(const tp2_bitvec *bv, uint64_t i)
137 : : {
138 [ - + ]: 56373 : if (!bv)
139 : 0 : return 0;
140 [ - + ]: 56373 : if (i > bv->nbits)
141 : 0 : i = bv->nbits;
142 : 56373 : return i - tp2_bv_rank1(bv, i);
143 : : }
144 : :
145 : : /* §5.3: binary search over the rank index, no extra stored structure. The
146 : : predicate is monotonic in i, which is what makes the search valid. */
147 : 5501 : static uint64_t bv_select(const tp2_bitvec *bv, uint64_t k, int want_one)
148 : : {
149 [ + - + + ]: 5501 : if (!bv || bv->nbits == 0)
150 : 1 : return 0;
151 [ + + ]: 5500 : uint64_t available = want_one ? tp2_bv_rank1(bv, bv->nbits) : tp2_bv_rank0(bv, bv->nbits);
152 [ + + ]: 5500 : if (k >= available)
153 : 3 : return bv->nbits; /* defined "not found" */
154 : :
155 : 5497 : uint64_t lo = 0, hi = bv->nbits;
156 [ + + ]: 71170 : while (lo < hi) {
157 : 65673 : uint64_t mid = lo + (hi - lo) / 2;
158 [ + + ]: 65673 : uint64_t seen = want_one ? tp2_bv_rank1(bv, mid + 1) : tp2_bv_rank0(bv, mid + 1);
159 [ + + ]: 65673 : if (seen <= k)
160 : 31858 : lo = mid + 1;
161 : : else
162 : 33815 : hi = mid;
163 : : }
164 : 5497 : return lo;
165 : : }
166 : :
167 : 1121 : uint64_t tp2_bv_select1(const tp2_bitvec *bv, uint64_t k)
168 : : {
169 : 1121 : return bv_select(bv, k, 1);
170 : : }
171 : :
172 : 4380 : uint64_t tp2_bv_select0(const tp2_bitvec *bv, uint64_t k)
173 : : {
174 : 4380 : return bv_select(bv, k, 0);
175 : : }
176 : :
177 : 1 : uint64_t tp2_bv_index_bytes(const tp2_bitvec *bv)
178 : : {
179 [ - + ]: 1 : if (!bv)
180 : 0 : return 0;
181 : 1 : return bv->nsuper * sizeof(uint64_t) + bv->nblock * sizeof(uint16_t);
182 : : }
183 : :
184 : : /* ── LOUDS (§6.2) ────────────────────────────────────────────────────── */
185 : :
186 : 462 : uint64_t tp2_louds_child_begin(const tp2_bitvec *bv, uint64_t v)
187 : : {
188 [ + + ]: 462 : return (v == 0) ? 0 : tp2_bv_select0(bv, v - 1) + 1;
189 : : }
190 : :
191 : 233 : uint64_t tp2_louds_child_count(const tp2_bitvec *bv, uint64_t v)
192 : : {
193 : 233 : uint64_t begin = tp2_louds_child_begin(bv, v);
194 : 233 : uint64_t end = tp2_bv_select0(bv, v);
195 [ + + ]: 233 : return (end > begin) ? end - begin : 0;
196 : : }
197 : :
198 : 225 : uint64_t tp2_louds_child(const tp2_bitvec *bv, uint64_t v, uint64_t k)
199 : : {
200 : 225 : return tp2_bv_rank1(bv, tp2_louds_child_begin(bv, v) + k) + 1;
201 : : }
|
