simple_sds_sbwt/
bits.rs

1//! Low-level functions for bit manipulation.
2//!
3//! Most operations are built around storing bits in integer arrays of type `u64` using the least significant bits first.
4
5use std::ops::{Index, IndexMut};
6
7//-----------------------------------------------------------------------------
8
9/// Number of bytes in [`u64`].
10pub const WORD_BYTES: usize = 8;
11
12/// Number of bits in [`u64`].
13pub const WORD_BITS: usize = 64;
14
15// Bit shift for transforming a bit offset into an array index.
16const INDEX_SHIFT: usize = 6;
17
18// Bit mask for transforming a bit offset into an offset in [`u64`].
19const OFFSET_MASK: usize = 0b111111;
20
21//-----------------------------------------------------------------------------
22
23const LOW_SET: [u64; 65] = [
24    0x0000_0000_0000_0000,
25
26    0x0000_0000_0000_0001, 0x0000_0000_0000_0003, 0x0000_0000_0000_0007, 0x0000_0000_0000_000F,
27    0x0000_0000_0000_001F, 0x0000_0000_0000_003F, 0x0000_0000_0000_007F, 0x0000_0000_0000_00FF,
28    0x0000_0000_0000_01FF, 0x0000_0000_0000_03FF, 0x0000_0000_0000_07FF, 0x0000_0000_0000_0FFF,
29    0x0000_0000_0000_1FFF, 0x0000_0000_0000_3FFF, 0x0000_0000_0000_7FFF, 0x0000_0000_0000_FFFF,
30
31    0x0000_0000_0001_FFFF, 0x0000_0000_0003_FFFF, 0x0000_0000_0007_FFFF, 0x0000_0000_000F_FFFF,
32    0x0000_0000_001F_FFFF, 0x0000_0000_003F_FFFF, 0x0000_0000_007F_FFFF, 0x0000_0000_00FF_FFFF,
33    0x0000_0000_01FF_FFFF, 0x0000_0000_03FF_FFFF, 0x0000_0000_07FF_FFFF, 0x0000_0000_0FFF_FFFF,
34    0x0000_0000_1FFF_FFFF, 0x0000_0000_3FFF_FFFF, 0x0000_0000_7FFF_FFFF, 0x0000_0000_FFFF_FFFF,
35
36    0x0000_0001_FFFF_FFFF, 0x0000_0003_FFFF_FFFF, 0x0000_0007_FFFF_FFFF, 0x0000_000F_FFFF_FFFF,
37    0x0000_001F_FFFF_FFFF, 0x0000_003F_FFFF_FFFF, 0x0000_007F_FFFF_FFFF, 0x0000_00FF_FFFF_FFFF,
38    0x0000_01FF_FFFF_FFFF, 0x0000_03FF_FFFF_FFFF, 0x0000_07FF_FFFF_FFFF, 0x0000_0FFF_FFFF_FFFF,
39    0x0000_1FFF_FFFF_FFFF, 0x0000_3FFF_FFFF_FFFF, 0x0000_7FFF_FFFF_FFFF, 0x0000_FFFF_FFFF_FFFF,
40
41    0x0001_FFFF_FFFF_FFFF, 0x0003_FFFF_FFFF_FFFF, 0x0007_FFFF_FFFF_FFFF, 0x000F_FFFF_FFFF_FFFF,
42    0x001F_FFFF_FFFF_FFFF, 0x003F_FFFF_FFFF_FFFF, 0x007F_FFFF_FFFF_FFFF, 0x00FF_FFFF_FFFF_FFFF,
43    0x01FF_FFFF_FFFF_FFFF, 0x03FF_FFFF_FFFF_FFFF, 0x07FF_FFFF_FFFF_FFFF, 0x0FFF_FFFF_FFFF_FFFF,
44    0x1FFF_FFFF_FFFF_FFFF, 0x3FFF_FFFF_FFFF_FFFF, 0x7FFF_FFFF_FFFF_FFFF, 0xFFFF_FFFF_FFFF_FFFF,
45];
46
47const HIGH_SET: [u64; 65] = [
48    0x0000_0000_0000_0000,
49
50    0x8000_0000_0000_0000, 0xC000_0000_0000_0000, 0xE000_0000_0000_0000, 0xF000_0000_0000_0000,
51    0xF800_0000_0000_0000, 0xFC00_0000_0000_0000, 0xFE00_0000_0000_0000, 0xFF00_0000_0000_0000,
52    0xFF80_0000_0000_0000, 0xFFC0_0000_0000_0000, 0xFFE0_0000_0000_0000, 0xFFF0_0000_0000_0000,
53    0xFFF8_0000_0000_0000, 0xFFFC_0000_0000_0000, 0xFFFE_0000_0000_0000, 0xFFFF_0000_0000_0000,
54
55    0xFFFF_8000_0000_0000, 0xFFFF_C000_0000_0000, 0xFFFF_E000_0000_0000, 0xFFFF_F000_0000_0000,
56    0xFFFF_F800_0000_0000, 0xFFFF_FC00_0000_0000, 0xFFFF_FE00_0000_0000, 0xFFFF_FF00_0000_0000,
57    0xFFFF_FF80_0000_0000, 0xFFFF_FFC0_0000_0000, 0xFFFF_FFE0_0000_0000, 0xFFFF_FFF0_0000_0000,
58    0xFFFF_FFF8_0000_0000, 0xFFFF_FFFC_0000_0000, 0xFFFF_FFFE_0000_0000, 0xFFFF_FFFF_0000_0000,
59
60    0xFFFF_FFFF_8000_0000, 0xFFFF_FFFF_C000_0000, 0xFFFF_FFFF_E000_0000, 0xFFFF_FFFF_F000_0000,
61    0xFFFF_FFFF_F800_0000, 0xFFFF_FFFF_FC00_0000, 0xFFFF_FFFF_FE00_0000, 0xFFFF_FFFF_FF00_0000,
62    0xFFFF_FFFF_FF80_0000, 0xFFFF_FFFF_FFC0_0000, 0xFFFF_FFFF_FFE0_0000, 0xFFFF_FFFF_FFF0_0000,
63    0xFFFF_FFFF_FFF8_0000, 0xFFFF_FFFF_FFFC_0000, 0xFFFF_FFFF_FFFE_0000, 0xFFFF_FFFF_FFFF_0000,
64
65    0xFFFF_FFFF_FFFF_8000, 0xFFFF_FFFF_FFFF_C000, 0xFFFF_FFFF_FFFF_E000, 0xFFFF_FFFF_FFFF_F000,
66    0xFFFF_FFFF_FFFF_F800, 0xFFFF_FFFF_FFFF_FC00, 0xFFFF_FFFF_FFFF_FE00, 0xFFFF_FFFF_FFFF_FF00,
67    0xFFFF_FFFF_FFFF_FF80, 0xFFFF_FFFF_FFFF_FFC0, 0xFFFF_FFFF_FFFF_FFE0, 0xFFFF_FFFF_FFFF_FFF0,
68    0xFFFF_FFFF_FFFF_FFF8, 0xFFFF_FFFF_FFFF_FFFC, 0xFFFF_FFFF_FFFF_FFFE, 0xFFFF_FFFF_FFFF_FFFF,
69];
70
71//-----------------------------------------------------------------------------
72
73// Lookup tables for the `select` implementation without the PDEP instruction.
74
75// Each byte in `_PS_OVERFLOW[i]` contains `128 - i`.
76const _PS_OVERFLOW: [u64; 65] = [
77    0x8080_8080_8080_8080,
78
79    0x7F7F_7F7F_7F7F_7F7F, 0x7E7E_7E7E_7E7E_7E7E, 0x7D7D_7D7D_7D7D_7D7D, 0x7C7C_7C7C_7C7C_7C7C,
80    0x7B7B_7B7B_7B7B_7B7B, 0x7A7A_7A7A_7A7A_7A7A, 0x7979_7979_7979_7979, 0x7878_7878_7878_7878,
81    0x7777_7777_7777_7777, 0x7676_7676_7676_7676, 0x7575_7575_7575_7575, 0x7474_7474_7474_7474,
82    0x7373_7373_7373_7373, 0x7272_7272_7272_7272, 0x7171_7171_7171_7171, 0x7070_7070_7070_7070,
83
84    0x6F6F_6F6F_6F6F_6F6F, 0x6E6E_6E6E_6E6E_6E6E, 0x6D6D_6D6D_6D6D_6D6D, 0x6C6C_6C6C_6C6C_6C6C,
85    0x6B6B_6B6B_6B6B_6B6B, 0x6A6A_6A6A_6A6A_6A6A, 0x6969_6969_6969_6969, 0x6868_6868_6868_6868,
86    0x6767_6767_6767_6767, 0x6666_6666_6666_6666, 0x6565_6565_6565_6565, 0x6464_6464_6464_6464,
87    0x6363_6363_6363_6363, 0x6262_6262_6262_6262, 0x6161_6161_6161_6161, 0x6060_6060_6060_6060,
88
89    0x5F5F_5F5F_5F5F_5F5F, 0x5E5E_5E5E_5E5E_5E5E, 0x5D5D_5D5D_5D5D_5D5D, 0x5C5C_5C5C_5C5C_5C5C,
90    0x5B5B_5B5B_5B5B_5B5B, 0x5A5A_5A5A_5A5A_5A5A, 0x5959_5959_5959_5959, 0x5858_5858_5858_5858,
91    0x5757_5757_5757_5757, 0x5656_5656_5656_5656, 0x5555_5555_5555_5555, 0x5454_5454_5454_5454,
92    0x5353_5353_5353_5353, 0x5252_5252_5252_5252, 0x5151_5151_5151_5151, 0x5050_5050_5050_5050,
93
94    0x4F4F_4F4F_4F4F_4F4F, 0x4E4E_4E4E_4E4E_4E4E, 0x4D4D_4D4D_4D4D_4D4D, 0x4C4C_4C4C_4C4C_4C4C,
95    0x4B4B_4B4B_4B4B_4B4B, 0x4A4A_4A4A_4A4A_4A4A, 0x4949_4949_4949_4949, 0x4848_4848_4848_4848,
96    0x4747_4747_4747_4747, 0x4646_4646_4646_4646, 0x4545_4545_4545_4545, 0x4444_4444_4444_4444,
97    0x4343_4343_4343_4343, 0x4242_4242_4242_4242, 0x4141_4141_4141_4141, 0x4040_4040_4040_4040,
98];
99
100// `_SELECT_IN_BYTE[256 * i + x]` is the position of the `i`th set bit in `x`.
101const _SELECT_IN_BYTE: [u8; 2048] = [
102    0, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
103    4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
104    5, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
105    4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
106    6, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
107    4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
108    5, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
109    4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
110    7, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
111    4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
112    5, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
113    4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
114    6, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
115    4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
116    5, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
117    4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
118
119    0, 0, 0, 1, 0, 2, 2, 1, 0, 3, 3, 1, 3, 2, 2, 1,
120    0, 4, 4, 1, 4, 2, 2, 1, 4, 3, 3, 1, 3, 2, 2, 1,
121    0, 5, 5, 1, 5, 2, 2, 1, 5, 3, 3, 1, 3, 2, 2, 1,
122    5, 4, 4, 1, 4, 2, 2, 1, 4, 3, 3, 1, 3, 2, 2, 1,
123    0, 6, 6, 1, 6, 2, 2, 1, 6, 3, 3, 1, 3, 2, 2, 1,
124    6, 4, 4, 1, 4, 2, 2, 1, 4, 3, 3, 1, 3, 2, 2, 1,
125    6, 5, 5, 1, 5, 2, 2, 1, 5, 3, 3, 1, 3, 2, 2, 1,
126    5, 4, 4, 1, 4, 2, 2, 1, 4, 3, 3, 1, 3, 2, 2, 1,
127    0, 7, 7, 1, 7, 2, 2, 1, 7, 3, 3, 1, 3, 2, 2, 1,
128    7, 4, 4, 1, 4, 2, 2, 1, 4, 3, 3, 1, 3, 2, 2, 1,
129    7, 5, 5, 1, 5, 2, 2, 1, 5, 3, 3, 1, 3, 2, 2, 1,
130    5, 4, 4, 1, 4, 2, 2, 1, 4, 3, 3, 1, 3, 2, 2, 1,
131    7, 6, 6, 1, 6, 2, 2, 1, 6, 3, 3, 1, 3, 2, 2, 1,
132    6, 4, 4, 1, 4, 2, 2, 1, 4, 3, 3, 1, 3, 2, 2, 1,
133    6, 5, 5, 1, 5, 2, 2, 1, 5, 3, 3, 1, 3, 2, 2, 1,
134    5, 4, 4, 1, 4, 2, 2, 1, 4, 3, 3, 1, 3, 2, 2, 1,
135
136    0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 3, 0, 3, 3, 2,
137    0, 0, 0, 4, 0, 4, 4, 2, 0, 4, 4, 3, 4, 3, 3, 2,
138    0, 0, 0, 5, 0, 5, 5, 2, 0, 5, 5, 3, 5, 3, 3, 2,
139    0, 5, 5, 4, 5, 4, 4, 2, 5, 4, 4, 3, 4, 3, 3, 2,
140    0, 0, 0, 6, 0, 6, 6, 2, 0, 6, 6, 3, 6, 3, 3, 2,
141    0, 6, 6, 4, 6, 4, 4, 2, 6, 4, 4, 3, 4, 3, 3, 2,
142    0, 6, 6, 5, 6, 5, 5, 2, 6, 5, 5, 3, 5, 3, 3, 2,
143    6, 5, 5, 4, 5, 4, 4, 2, 5, 4, 4, 3, 4, 3, 3, 2,
144    0, 0, 0, 7, 0, 7, 7, 2, 0, 7, 7, 3, 7, 3, 3, 2,
145    0, 7, 7, 4, 7, 4, 4, 2, 7, 4, 4, 3, 4, 3, 3, 2,
146    0, 7, 7, 5, 7, 5, 5, 2, 7, 5, 5, 3, 5, 3, 3, 2,
147    7, 5, 5, 4, 5, 4, 4, 2, 5, 4, 4, 3, 4, 3, 3, 2,
148    0, 7, 7, 6, 7, 6, 6, 2, 7, 6, 6, 3, 6, 3, 3, 2,
149    7, 6, 6, 4, 6, 4, 4, 2, 6, 4, 4, 3, 4, 3, 3, 2,
150    7, 6, 6, 5, 6, 5, 5, 2, 6, 5, 5, 3, 5, 3, 3, 2,
151    6, 5, 5, 4, 5, 4, 4, 2, 5, 4, 4, 3, 4, 3, 3, 2,
152
153    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3,
154    0, 0, 0, 0, 0, 0, 0, 4, 0, 0, 0, 4, 0, 4, 4, 3,
155    0, 0, 0, 0, 0, 0, 0, 5, 0, 0, 0, 5, 0, 5, 5, 3,
156    0, 0, 0, 5, 0, 5, 5, 4, 0, 5, 5, 4, 5, 4, 4, 3,
157    0, 0, 0, 0, 0, 0, 0, 6, 0, 0, 0, 6, 0, 6, 6, 3,
158    0, 0, 0, 6, 0, 6, 6, 4, 0, 6, 6, 4, 6, 4, 4, 3,
159    0, 0, 0, 6, 0, 6, 6, 5, 0, 6, 6, 5, 6, 5, 5, 3,
160    0, 6, 6, 5, 6, 5, 5, 4, 6, 5, 5, 4, 5, 4, 4, 3,
161    0, 0, 0, 0, 0, 0, 0, 7, 0, 0, 0, 7, 0, 7, 7, 3,
162    0, 0, 0, 7, 0, 7, 7, 4, 0, 7, 7, 4, 7, 4, 4, 3,
163    0, 0, 0, 7, 0, 7, 7, 5, 0, 7, 7, 5, 7, 5, 5, 3,
164    0, 7, 7, 5, 7, 5, 5, 4, 7, 5, 5, 4, 5, 4, 4, 3,
165    0, 0, 0, 7, 0, 7, 7, 6, 0, 7, 7, 6, 7, 6, 6, 3,
166    0, 7, 7, 6, 7, 6, 6, 4, 7, 6, 6, 4, 6, 4, 4, 3,
167    0, 7, 7, 6, 7, 6, 6, 5, 7, 6, 6, 5, 6, 5, 5, 3,
168    7, 6, 6, 5, 6, 5, 5, 4, 6, 5, 5, 4, 5, 4, 4, 3,
169
170    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
171    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4,
172    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5,
173    0, 0, 0, 0, 0, 0, 0, 5, 0, 0, 0, 5, 0, 5, 5, 4,
174    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6,
175    0, 0, 0, 0, 0, 0, 0, 6, 0, 0, 0, 6, 0, 6, 6, 4,
176    0, 0, 0, 0, 0, 0, 0, 6, 0, 0, 0, 6, 0, 6, 6, 5,
177    0, 0, 0, 6, 0, 6, 6, 5, 0, 6, 6, 5, 6, 5, 5, 4,
178    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 7,
179    0, 0, 0, 0, 0, 0, 0, 7, 0, 0, 0, 7, 0, 7, 7, 4,
180    0, 0, 0, 0, 0, 0, 0, 7, 0, 0, 0, 7, 0, 7, 7, 5,
181    0, 0, 0, 7, 0, 7, 7, 5, 0, 7, 7, 5, 7, 5, 5, 4,
182    0, 0, 0, 0, 0, 0, 0, 7, 0, 0, 0, 7, 0, 7, 7, 6,
183    0, 0, 0, 7, 0, 7, 7, 6, 0, 7, 7, 6, 7, 6, 6, 4,
184    0, 0, 0, 7, 0, 7, 7, 6, 0, 7, 7, 6, 7, 6, 6, 5,
185    0, 7, 7, 6, 7, 6, 6, 5, 7, 6, 6, 5, 6, 5, 5, 4,
186
187    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
188    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
189    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
190    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5,
191    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
192    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6,
193    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6,
194    0, 0, 0, 0, 0, 0, 0, 6, 0, 0, 0, 6, 0, 6, 6, 5,
195    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
196    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 7,
197    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 7,
198    0, 0, 0, 0, 0, 0, 0, 7, 0, 0, 0, 7, 0, 7, 7, 5,
199    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 7,
200    0, 0, 0, 0, 0, 0, 0, 7, 0, 0, 0, 7, 0, 7, 7, 6,
201    0, 0, 0, 0, 0, 0, 0, 7, 0, 0, 0, 7, 0, 7, 7, 6,
202    0, 0, 0, 7, 0, 7, 7, 6, 0, 7, 7, 6, 7, 6, 6, 5,
203
204    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
205    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
206    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
207    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
208    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
209    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
210    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
211    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6,
212    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
213    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
214    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
215    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 7,
216    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
217    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 7,
218    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 7,
219    0, 0, 0, 0, 0, 0, 0, 7, 0, 0, 0, 7, 0, 7, 7, 6,
220
221    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
222    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
223    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
224    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
225    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
226    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
227    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
228    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
229    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
230    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
231    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
232    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
233    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
234    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
235    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
236    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 7,
237];
238
239//-----------------------------------------------------------------------------
240
241/// Returns an integer with the lowest `n` bits set.
242///
243/// # Examples
244///
245/// ```
246/// use simple_sds_sbwt::bits;
247///
248/// assert_eq!(bits::low_set(13), 0x1FFF);
249/// ```
250///
251/// # Panics
252///
253/// May panic if `n > 64`.
254#[inline]
255pub fn low_set(n: usize) -> u64 {
256    LOW_SET[n]
257}
258
259/// Unsafe version of [`low_set`].
260///
261/// # Safety
262///
263/// Behavior is undefined if `n > 64`.
264#[inline]
265pub unsafe fn low_set_unchecked(n: usize) -> u64 {
266    *LOW_SET.get_unchecked(n)
267}
268
269/// Returns an integer with the highest `n` bits set.
270///
271/// # Examples
272///
273/// ```
274/// use simple_sds_sbwt::bits;
275///
276/// assert_eq!(bits::high_set(13), 0xFFF8_0000_0000_0000);
277/// ```
278///
279/// # Panics
280///
281/// May panic if `n > 64`.
282#[inline]
283pub fn high_set(n: usize) -> u64 {
284    HIGH_SET[n]
285}
286
287/// Unsafe version of [`high_set`].
288///
289/// # Safety
290///
291/// Behavior is undefined if `n > 64`.
292#[inline]
293pub unsafe fn high_set_unchecked(n: usize) -> u64 {
294    *HIGH_SET.get_unchecked(n)
295}
296
297/// Returns the length of the binary representation of integer `n`.
298///
299/// # Examples
300///
301/// ```
302/// use simple_sds_sbwt::bits;
303///
304/// assert_eq!(bits::bit_len(0), 1);
305/// assert_eq!(bits::bit_len(0x1FFF), 13);
306/// ```
307#[inline]
308pub fn bit_len(n: u64) -> usize {
309    // The `| 1` to avoid a branch was obvious in hindsight.
310    WORD_BITS - ((n | 1).leading_zeros() as usize)
311}
312
313/// Returns the bit offset of the set bit of specified rank.
314///
315/// # Arguments
316///
317/// * `n`: An integer.
318/// * `rank`: Rank of the set bit we want to find.
319///
320/// # Examples
321///
322/// ```
323/// use simple_sds_sbwt::bits;
324///
325/// unsafe {
326///     assert_eq!(bits::select(0b00100001_00010000, 0), 4);
327///     assert_eq!(bits::select(0b00100001_00010000, 1), 8);
328///     assert_eq!(bits::select(0b00100001_00010000, 2), 13);
329/// }
330/// ```
331///
332/// # Safety
333///
334/// Behavior is undefined if `rank >= n.count_ones()`.
335#[inline]
336pub unsafe fn select(n: u64, rank: usize) -> usize {
337    // The first argument to `__pdep_u64` has a single 1 at bit offset `rank`. The
338    // number `n` we are interested in is used as a mask. PDEP takes low-order bits
339    // from the value and places them to the offsets specified by the mask. In
340    // particular, the only 1 is placed to the bit offset of the set bit of rank
341    // `rank` in `n`. We get the offset by counting the trailing zeros.
342    #[cfg(all(target_arch = "x86_64", target_feature = "bmi2"))]
343    {
344        let pos = core::arch::x86_64::_pdep_u64(1u64 << rank, n);
345        pos.trailing_zeros() as usize
346    }
347
348    // This is borrowed from SDSL.
349    #[cfg(not(all(target_arch = "x86_64", target_feature = "bmi2")))]
350    {
351        // Each byte in `cumulative` will contain the cumulative number of set bits
352        // in bytes up to and including that byte in `n`.
353        let cumulative = n - ((n >> 1) & 0x5555_5555_5555_5555);
354        let cumulative = (cumulative & 0x3333_3333_3333_3333) + ((cumulative >> 2) & 0x3333_3333_3333_3333);
355        let cumulative = (cumulative + (cumulative >> 4)) & 0x0F0F_0F0F_0F0F_0F0F;
356        let (cumulative, _) = cumulative.overflowing_mul(0x0101_0101_0101_0101);
357
358        // We add `128 - rank - 1` to each byte and mask out all bits except `128`. We get
359        // the bit offset for the byte containing the answer by counting trailing zeros.
360        let mask = (cumulative + *_PS_OVERFLOW.get_unchecked(rank + 1)) & 0x8080_8080_8080_8080;
361        let offset = ((mask.trailing_zeros() >> 3) << 3) as usize;
362
363        // Subtract the number of set bits in the previous bytes from the rank.
364        let relative_rank = rank - (((cumulative << 8) >> offset) as usize & 0xFF);
365
366        offset + (*_SELECT_IN_BYTE.get_unchecked((relative_rank << 8) + ((n >> offset) as usize & 0xFF)) as usize)
367    }
368}
369
370//-----------------------------------------------------------------------------
371
372/// Returns the number of bytes that can be stored in `n` integers of type [`u64`].
373///
374/// # Examples
375///
376/// ```
377/// use simple_sds_sbwt::bits;
378///
379/// assert_eq!(bits::words_to_bytes(3), 24);
380/// ```
381///
382/// # Panics
383///
384/// May panic if `n * 8 > usize::MAX`.
385#[inline]
386pub fn words_to_bytes(n: usize) -> usize {
387    n * WORD_BYTES
388}
389
390/// Returns the number of integers of type [`u64`] required to store `n` bytes.
391///
392/// # Examples
393///
394/// ```
395/// use simple_sds_sbwt::bits;
396///
397/// assert_eq!(bits::bytes_to_words(8), 1);
398/// assert_eq!(bits::bytes_to_words(9), 2);
399/// ```
400///
401/// # Panics
402///
403/// May panic if `n + 7 > usize::MAX`.
404#[inline]
405pub fn bytes_to_words(n: usize) -> usize {
406    (n + WORD_BYTES - 1) / WORD_BYTES
407}
408
409/// Rounds `n` up to the next multiple of 8.
410///
411/// # Examples
412///
413/// ```
414/// use simple_sds_sbwt::bits;
415///
416/// assert_eq!(bits::round_up_to_word_bytes(0), 0);
417/// assert_eq!(bits::round_up_to_word_bytes(8), 8);
418/// assert_eq!(bits::round_up_to_word_bytes(9), 16);
419/// ```
420///
421/// # Panics
422///
423/// May panic if `n + 7 > usize::MAX`.
424#[inline]
425pub fn round_up_to_word_bytes(n: usize) -> usize {
426    words_to_bytes(bytes_to_words(n))
427}
428
429//-----------------------------------------------------------------------------
430
431/// Returns the number of bits that can be stored in `n` integers of type [`u64`].
432///
433/// # Examples
434///
435/// ```
436/// use simple_sds_sbwt::bits;
437///
438/// assert_eq!(bits::words_to_bits(3), 192);
439/// ```
440///
441/// # Panics
442///
443/// May panic if `n * 64 > usize::MAX`.
444#[inline]
445pub fn words_to_bits(n: usize) -> usize {
446    n * WORD_BITS
447}
448
449/// Returns the number of integers of type [`u64`] required to store `n` bits.
450///
451/// # Examples
452///
453/// ```
454/// use simple_sds_sbwt::bits;
455///
456/// assert_eq!(bits::bits_to_words(64), 1);
457/// assert_eq!(bits::bits_to_words(65), 2);
458/// ```
459///
460/// # Panics
461///
462/// May panic if `n + 63 > usize::MAX`.
463#[inline]
464pub fn bits_to_words(n: usize) -> usize {
465    (n + WORD_BITS - 1) / WORD_BITS
466}
467
468/// Rounds `n` up to the next multiple of 64.
469///
470/// # Examples
471///
472/// ```
473/// use simple_sds_sbwt::bits;
474///
475/// assert_eq!(bits::round_up_to_word_bits(0), 0);
476/// assert_eq!(bits::round_up_to_word_bits(64), 64);
477/// assert_eq!(bits::round_up_to_word_bits(65), 128);
478/// ```
479///
480/// # Panics
481///
482/// May panic if `n + 63 > usize::MAX`.
483#[inline]
484pub fn round_up_to_word_bits(n: usize) -> usize {
485    words_to_bits(bits_to_words(n))
486}
487
488//-----------------------------------------------------------------------------
489
490/// Divides `value` by `n` and rounds the result up.
491///
492/// # Examples
493///
494/// ```
495/// use simple_sds_sbwt::bits;
496///
497/// assert_eq!(bits::div_round_up(129, 13), 10);
498/// assert_eq!(bits::div_round_up(130, 13), 10);
499/// assert_eq!(bits::div_round_up(131, 13), 11);
500/// ```
501///
502/// # Panics
503///
504/// May panic if `value + n > usize::MAX` or `n == 0`.
505#[inline]
506pub fn div_round_up(value: usize, n: usize) -> usize {
507    (value + n - 1) / n
508}
509
510/// Returns a [`u64`] value consisting entirely of bit `value`.
511///
512/// # Examples
513///
514/// ```
515/// use simple_sds_sbwt::bits;
516///
517/// assert_eq!(bits::filler_value(false), 0);
518/// assert_eq!(bits::filler_value(true), !0u64);
519/// ```
520#[inline]
521pub fn filler_value(value: bool) -> u64 {
522    match value {
523        true  => !0u64,
524        false => 0u64,
525    }
526}
527
528/// Splits a bit offset into an index in an array of [`u64`] and an offset within the integer.
529///
530/// # Examples
531///
532/// ```
533/// use simple_sds_sbwt::bits;
534///
535/// assert_eq!(bits::split_offset(123), (1, 59));
536/// ```
537#[inline]
538pub fn split_offset(bit_offset: usize) -> (usize, usize) {
539    (bit_offset >> INDEX_SHIFT, bit_offset & OFFSET_MASK)
540}
541
542/// Combines an index in an array of [`u64`] and an offset within the integer into a bit offset.
543///
544/// # Arguments
545///
546/// * `index`: Array index.
547/// * `offset`: Offset within the integer.
548///
549/// # Examples
550///
551/// ```
552/// use simple_sds_sbwt::bits;
553///
554/// assert_eq!(bits::bit_offset(1, 59), 123);
555/// ```
556///
557/// # Panics
558///
559/// May panic if the result would be greater than [`usize::MAX`].
560#[inline]
561pub fn bit_offset(index: usize, offset: usize) -> usize {
562    (index << INDEX_SHIFT) + offset
563}
564
565//-----------------------------------------------------------------------------
566
567/// Writes an integer into a bit array implemented as an array of [`u64`] values.
568///
569/// # Arguments
570///
571/// * `bit_offset`: Starting offset in the bit array.
572/// * `value`: The integer to be written.
573/// * `width`: The width of the integer in bits.
574///
575/// # Examples
576///
577/// ```
578/// use simple_sds_sbwt::bits;
579///
580/// let mut array: Vec<u64> = vec![0];
581/// unsafe {
582///     bits::write_int(&mut array, 0, 7, 4);
583///     bits::write_int(&mut array, 4, 3, 4);
584/// }
585/// assert_eq!(array[0], 0x37);
586/// ```
587///
588/// # Safety
589///
590/// Behavior is undefined if `width > 64`.
591///
592/// # Panics
593///
594/// May panic if `(bit_offset + width - 1) / 64` is not a valid index in the array.
595pub unsafe fn write_int<T: IndexMut<usize, Output = u64>>(array: &mut T, bit_offset: usize, value: u64, width: usize) {
596    let value = value & low_set(width);
597    let (index, offset) = split_offset(bit_offset);
598
599    if offset + width <= WORD_BITS {
600        array[index] &= high_set_unchecked(WORD_BITS - width - offset) | low_set_unchecked(offset);
601        array[index] |= value << offset;
602    } else {
603        array[index] &= low_set_unchecked(offset);
604        array[index] |= value << offset;
605        array[index + 1] &= high_set(2 * WORD_BITS - width - offset);
606        array[index + 1] |= value >> (WORD_BITS - offset);
607    }
608}
609
610/// Reads an integer from a bit array implemented as an array of [`u64`] values.
611///
612/// # Arguments
613///
614/// * `bit_offset`: Starting offset in the bit array.
615/// * `width`: The width of the integer in bits.
616///
617/// # Examples
618///
619/// ```
620/// use simple_sds_sbwt::bits;
621///
622/// let array: Vec<u64> = vec![0x37];
623/// unsafe {
624///     assert_eq!(bits::read_int(&array, 0, 4), 7);
625///     assert_eq!(bits::read_int(&array, 4, 4), 3);
626/// }
627/// ```
628///
629/// # Safety
630///
631/// Behavior is undefined if `width > 64`.
632///
633/// # Panics
634///
635/// May panic if `(bit_offset + width - 1) / 64` is not a valid index in the array.
636pub unsafe fn read_int<T: Index<usize, Output = u64>>(array: &T, bit_offset: usize, width: usize) -> u64 {
637    let (index, offset) = split_offset(bit_offset);
638    let first = array[index] >> offset;
639
640    if offset + width <= WORD_BITS {
641        first & low_set_unchecked(width)
642    } else {
643        first | ((array[index + 1] & low_set_unchecked((offset + width) & OFFSET_MASK)) << (WORD_BITS - offset))
644    }
645}
646
647//-----------------------------------------------------------------------------
648
649#[cfg(test)]
650mod tests {
651    use super::*;
652    use std::mem;
653    use rand::Rng;
654
655    #[test]
656    fn low_set_test() {
657        assert_eq!(low_set(0), 0u64, "low_set(0) failed");
658        assert_eq!(low_set(13), 0x1FFF, "low_set(13) failed");
659        assert_eq!(low_set(64), !0u64, "low_set(64) failed");
660
661        unsafe {
662            assert_eq!(low_set_unchecked(0), 0u64, "low_set_unchecked(0) failed");
663            assert_eq!(low_set_unchecked(13), 0x1FFF, "low_set_unchecked(13) failed");
664            assert_eq!(low_set_unchecked(64), !0u64, "low_set_unchecked(64) failed")
665        }
666    }
667
668    #[test]
669    fn high_set_test() {
670        assert_eq!(high_set(0), 0, "high_set(0) failed");
671        assert_eq!(high_set(13), 0xFFF8_0000_0000_0000, "high_set(13) failed");
672        assert_eq!(high_set(64), !0u64, "high_set(64) failed");
673
674        unsafe {
675            assert_eq!(high_set_unchecked(0), 0, "high_set_unchecked(0) failed");
676            assert_eq!(high_set_unchecked(13), 0xFFF8_0000_0000_0000, "high_set_unchecked(13) failed");
677            assert_eq!(high_set_unchecked(64), !0u64, "high_set_unchecked(64) failed")
678        }
679    }
680
681    #[test]
682    fn bit_len_test() {
683        assert_eq!(bit_len(0), 1, "bit_len(0) failed");
684        assert_eq!(bit_len(1), 1, "bit_len(1) failed");
685        assert_eq!(bit_len(0x1FFF), 13, "bit_len(0x1FFF) failed");
686        assert_eq!(bit_len(0x7FFF_FFFF_FFFF_FFFF), 63, "bit_len(0x7FFF_FFFF_FFFF_FFFF) failed");
687        assert_eq!(bit_len(0xFFFF_FFFF_FFFF_FFFF), 64, "bit_len(0xFFFF_FFFF_FFFF_FFFF) failed");
688    }
689
690    #[test]
691    fn select_test() {
692        let values: Vec<(u64, usize, usize)> = vec![
693        (0x0000_0000_0000_0001, 0, 0),
694        (0x8000_0000_0000_0000, 0, 63),
695        (0x8000_0000_0000_0001, 0, 0),
696        (0x8000_0000_0000_0001, 1, 63),
697        (0x8000_0010_0000_0001, 0, 0),
698        (0x8000_0010_0000_0001, 1, 36),
699        (0x8000_0010_0000_0001, 2, 63),
700        ];
701        for (value, rank, result) in values.iter() {
702            unsafe { assert_eq!(select(*value, *rank), *result, "select({:X}, {}) failed", value, rank); }
703        }
704    }
705
706    #[test]
707    fn read_write() {
708        let mut correct: Vec<(u64, u64, u64, u64)> = Vec::new();
709        let mut rng = rand::thread_rng();
710        for _ in 0..64 {
711            let mut tuple: (u64, u64, u64, u64) = rng.gen();
712            tuple.0 &= low_set(31); tuple.1 &= low_set(64); tuple.2 &= low_set(35); tuple.3 &= low_set(63);
713            correct.push(tuple);
714        }
715
716        let mut array: Vec<u64> = vec![0; 256];
717        let mut bit_offset = 0;
718        for i in 0..64 {
719            unsafe {
720                write_int(&mut array, bit_offset, correct[i].0, 31); bit_offset += 31;
721                write_int(&mut array, bit_offset, correct[i].1, 64); bit_offset += 64;
722                write_int(&mut array, bit_offset, correct[i].2, 35); bit_offset += 35;
723                write_int(&mut array, bit_offset, correct[i].3, 63); bit_offset += 63;
724            }
725        }
726
727        bit_offset = 0;
728        for i in 0..64 {
729            unsafe {
730                assert_eq!(read_int(&array, bit_offset, 31), correct[i].0, "Invalid value at [{}].0", i); bit_offset += 31;
731                assert_eq!(read_int(&array, bit_offset, 64), correct[i].1, "Invalid value at [{}].1", i); bit_offset += 64;
732                assert_eq!(read_int(&array, bit_offset, 35), correct[i].2, "Invalid value at [{}].2", i); bit_offset += 35;
733                assert_eq!(read_int(&array, bit_offset, 63), correct[i].3, "Invalid value at [{}].3", i); bit_offset += 63;
734            }
735        }
736    }
737
738    #[test]
739    fn no_extra_bits() {
740        let mut array: Vec<u64> = vec![0; 2];
741        unsafe {
742            write_int(&mut array, 16, 2, 16);
743            write_int(&mut array, 48, 2, 16);
744            write_int(&mut array, 32, !0u64, 16); // This should not overwrite the integer at offset 48.
745            write_int(&mut array, 0, !0u64, 16); // This should not overwrite the other integers.
746            assert_eq!(read_int(&array, 0, 16), 0xFFFF, "Incorrect 16-bit integer at offset 0");
747            assert_eq!(read_int(&array, 16, 16), 2, "Incorrect 16-bit integer at offset 16");
748            assert_eq!(read_int(&array, 32, 16), 0xFFFF, "Incorrect 16-bit integer at offset 32");
749            assert_eq!(read_int(&array, 48, 16), 2, "Incorrect 16-bit integer at offset 48");
750        }
751    }
752
753    #[test]
754    #[ignore]
755    fn environment() {
756        assert_eq!(mem::size_of::<usize>(), 8, "Things may not work if usize is not 64 bits");
757        assert_eq!(mem::align_of_val(&1usize), 8, "Things may not work if the minimum alignment of usize is not 8 bytes");
758        assert_eq!(mem::align_of_val(&1u64), 8, "Things may not work if the minimum alignment of u64 is not 8 bytes");
759        assert!(cfg!(target_endian = "little"), "Things may not work on a big-endian system");
760        assert!(cfg!(any(target_arch = "x86_64", target_arch = "aarch64")), "Target architecture should be x86_64 or aarch64");
761        assert!(cfg!(unix), "Memory mapping requires Unix-like OS");
762
763        #[cfg(target_arch = "x86_64")]
764        {
765            assert!(cfg!(target_feature = "popcnt"), "Performance may be worse without the POPCNT instruction");
766            assert!(cfg!(target_feature = "lzcnt"), "Performance may be worse without the LZCNT instruction");
767            assert!(cfg!(target_feature = "bmi1"), "Performance may be worse without the TZCNT instruction from BMI1");
768            assert!(cfg!(target_feature = "bmi2"), "Performance may be worse without the PDEP instruction from BMI2");
769        }
770    }
771}
772
773//-----------------------------------------------------------------------------