RapidFuzz
Loading...
Searching...
No Matches
OSA_impl.hpp
1
2/* SPDX-License-Identifier: MIT */
3/* Copyright © 2022-present Max Bachmann */
4
5#pragma once
6#include <cstdint>
7#include <rapidfuzz/details/PatternMatchVector.hpp>
8#include <rapidfuzz/details/Range.hpp>
9#include <rapidfuzz/details/common.hpp>
10#include <rapidfuzz/details/distance.hpp>
11#include <rapidfuzz/details/simd.hpp>
12
13namespace rapidfuzz {
14namespace detail {
15
34template <typename PM_Vec, typename InputIt1, typename InputIt2>
35size_t osa_hyrroe2003(const PM_Vec& PM, const Range<InputIt1>& s1, const Range<InputIt2>& s2, size_t max)
36{
37 /* VP is set to 1^m. Shifting by bitwidth would be undefined behavior */
38 uint64_t VP = ~UINT64_C(0);
39 uint64_t VN = 0;
40 uint64_t D0 = 0;
41 uint64_t PM_j_old = 0;
42 size_t currDist = s1.size();
43 assert(s1.size() != 0);
44
45 /* mask used when computing D[m,j] in the paper 10^(m-1) */
46 uint64_t mask = UINT64_C(1) << (s1.size() - 1);
47
48 /* Searching */
49 for (const auto& ch : s2) {
50 /* Step 1: Computing D0 */
51 uint64_t PM_j = PM.get(0, ch);
52 uint64_t TR = (((~D0) & PM_j) << 1) & PM_j_old;
53 D0 = (((PM_j & VP) + VP) ^ VP) | PM_j | VN;
54 D0 = D0 | TR;
55
56 /* Step 2: Computing HP and HN */
57 uint64_t HP = VN | ~(D0 | VP);
58 uint64_t HN = D0 & VP;
59
60 /* Step 3: Computing the value D[m,j] */
61 currDist += bool(HP & mask);
62 currDist -= bool(HN & mask);
63
64 /* Step 4: Computing Vp and VN */
65 HP = (HP << 1) | 1;
66 HN = (HN << 1);
67
68 VP = HN | ~(D0 | HP);
69 VN = HP & D0;
70 PM_j_old = PM_j;
71 }
72
73 return (currDist <= max) ? currDist : max + 1;
74}
75
76#ifdef RAPIDFUZZ_SIMD
77template <typename VecType, typename InputIt, int _lto_hack = RAPIDFUZZ_LTO_HACK>
78void osa_hyrroe2003_simd(Range<size_t*> scores, const detail::BlockPatternMatchVector& block,
79 const std::vector<size_t>& s1_lengths, const Range<InputIt>& s2,
80 size_t score_cutoff) noexcept
81{
82# ifdef RAPIDFUZZ_AVX2
83 using namespace simd_avx2;
84# else
85 using namespace simd_sse2;
86# endif
87 static constexpr size_t alignment = native_simd<VecType>::alignment;
88 static constexpr size_t vec_width = native_simd<VecType>::size;
89 static constexpr size_t vecs = native_simd<uint64_t>::size;
90 assert(block.size() % vecs == 0);
91
92 native_simd<VecType> zero(VecType(0));
93 native_simd<VecType> one(1);
94 size_t result_index = 0;
95
96 for (size_t cur_vec = 0; cur_vec < block.size(); cur_vec += vecs) {
97 /* VP is set to 1^m */
98 native_simd<VecType> VP(static_cast<VecType>(-1));
99 native_simd<VecType> VN(VecType(0));
100 native_simd<VecType> D0(VecType(0));
101 native_simd<VecType> PM_j_old(VecType(0));
102
103 alignas(alignment) std::array<VecType, vec_width> currDist_;
104 unroll<size_t, vec_width>(
105 [&](size_t i) { currDist_[i] = static_cast<VecType>(s1_lengths[result_index + i]); });
106 native_simd<VecType> currDist(reinterpret_cast<uint64_t*>(currDist_.data()));
107 /* mask used when computing D[m,j] in the paper 10^(m-1) */
108 alignas(alignment) std::array<VecType, vec_width> mask_;
109 unroll<size_t, vec_width>([&](size_t i) {
110 if (s1_lengths[result_index + i] == 0)
111 mask_[i] = 0;
112 else
113 mask_[i] = static_cast<VecType>(UINT64_C(1) << (s1_lengths[result_index + i] - 1));
114 });
115 native_simd<VecType> mask(reinterpret_cast<uint64_t*>(mask_.data()));
116
117 for (const auto& ch : s2) {
118 /* Step 1: Computing D0 */
119 alignas(alignment) std::array<uint64_t, vecs> stored;
120 unroll<size_t, vecs>([&](size_t i) { stored[i] = block.get(cur_vec + i, ch); });
121
122 native_simd<VecType> PM_j(stored.data());
123 auto TR = (andnot(PM_j, D0) << 1) & PM_j_old;
124 D0 = (((PM_j & VP) + VP) ^ VP) | PM_j | VN;
125 D0 = D0 | TR;
126
127 /* Step 2: Computing HP and HN */
128 auto HP = VN | ~(D0 | VP);
129 auto HN = D0 & VP;
130
131 /* Step 3: Computing the value D[m,j] */
132 currDist += andnot(one, (HP & mask) == zero);
133 currDist -= andnot(one, (HN & mask) == zero);
134
135 /* Step 4: Computing Vp and VN */
136 HP = (HP << 1) | one;
137 HN = (HN << 1);
138
139 VP = HN | ~(D0 | HP);
140 VN = HP & D0;
141 PM_j_old = PM_j;
142 }
143
144 alignas(alignment) std::array<VecType, vec_width> distances;
145 currDist.store(distances.data());
146
147 unroll<size_t, vec_width>([&](size_t i) {
148 size_t score = 0;
149 /* strings of length 0 are not handled correctly */
150 if (s1_lengths[result_index] == 0) {
151 score = s2.size();
152 }
153 /* calculate score under consideration of wraparounds in parallel counter */
154 else {
155 RAPIDFUZZ_IF_CONSTEXPR (std::numeric_limits<VecType>::max() <
156 std::numeric_limits<size_t>::max())
157 {
158 size_t min_dist = abs_diff(s1_lengths[result_index], s2.size());
159 size_t wraparound_score = static_cast<size_t>(std::numeric_limits<VecType>::max()) + 1;
160
161 score = (min_dist / wraparound_score) * wraparound_score;
162 VecType remainder = static_cast<VecType>(min_dist % wraparound_score);
163
164 if (distances[i] < remainder) score += wraparound_score;
165 }
166
167 score += distances[i];
168 }
169 scores[result_index] = (score <= score_cutoff) ? score : score_cutoff + 1;
170 result_index++;
171 });
172 }
173}
174#endif
175
176template <typename InputIt1, typename InputIt2>
177size_t osa_hyrroe2003_block(const BlockPatternMatchVector& PM, const Range<InputIt1>& s1,
178 const Range<InputIt2>& s2, size_t max = std::numeric_limits<size_t>::max())
179{
180 struct Row {
181 uint64_t VP;
182 uint64_t VN;
183 uint64_t D0;
184 uint64_t PM;
185
186 Row() : VP(~UINT64_C(0)), VN(0), D0(0), PM(0)
187 {}
188 };
189
190 size_t word_size = sizeof(uint64_t) * 8;
191 size_t words = PM.size();
192 uint64_t Last = UINT64_C(1) << ((s1.size() - 1) % word_size);
193
194 size_t currDist = s1.size();
195 std::vector<Row> old_vecs(words + 1);
196 std::vector<Row> new_vecs(words + 1);
197
198 /* Searching */
199 auto iter_s2 = s2.begin();
200 for (size_t row = 0; row < s2.size(); ++iter_s2, ++row) {
201 uint64_t HP_carry = 1;
202 uint64_t HN_carry = 0;
203
204 for (size_t word = 0; word < words; word++) {
205 /* retrieve bit vectors from last iterations */
206 uint64_t VN = old_vecs[word + 1].VN;
207 uint64_t VP = old_vecs[word + 1].VP;
208 uint64_t D0 = old_vecs[word + 1].D0;
209 /* D0 last word */
210 uint64_t D0_last = old_vecs[word].D0;
211
212 /* PM of last char same word */
213 uint64_t PM_j_old = old_vecs[word + 1].PM;
214 /* PM of last word */
215 uint64_t PM_last = new_vecs[word].PM;
216
217 uint64_t PM_j = PM.get(word, *iter_s2);
218 uint64_t X = PM_j;
219 uint64_t TR = ((((~D0) & X) << 1) | (((~D0_last) & PM_last) >> 63)) & PM_j_old;
220
221 X |= HN_carry;
222 D0 = (((X & VP) + VP) ^ VP) | X | VN | TR;
223
224 uint64_t HP = VN | ~(D0 | VP);
225 uint64_t HN = D0 & VP;
226
227 if (word == words - 1) {
228 currDist += bool(HP & Last);
229 currDist -= bool(HN & Last);
230 }
231
232 uint64_t HP_carry_temp = HP_carry;
233 HP_carry = HP >> 63;
234 HP = (HP << 1) | HP_carry_temp;
235 uint64_t HN_carry_temp = HN_carry;
236 HN_carry = HN >> 63;
237 HN = (HN << 1) | HN_carry_temp;
238
239 new_vecs[word + 1].VP = HN | ~(D0 | HP);
240 new_vecs[word + 1].VN = HP & D0;
241 new_vecs[word + 1].D0 = D0;
242 new_vecs[word + 1].PM = PM_j;
243 }
244
245 std::swap(new_vecs, old_vecs);
246 }
247
248 return (currDist <= max) ? currDist : max + 1;
249}
250
251class OSA : public DistanceBase<OSA, size_t, 0, std::numeric_limits<int64_t>::max()> {
252 friend DistanceBase<OSA, size_t, 0, std::numeric_limits<int64_t>::max()>;
253 friend NormalizedMetricBase<OSA>;
254
255 template <typename InputIt1, typename InputIt2>
256 static size_t maximum(const Range<InputIt1>& s1, const Range<InputIt2>& s2)
257 {
258 return std::max(s1.size(), s2.size());
259 }
260
261 template <typename InputIt1, typename InputIt2>
262 static size_t _distance(Range<InputIt1> s1, Range<InputIt2> s2, size_t score_cutoff, size_t score_hint)
263 {
264 if (s2.size() < s1.size()) return _distance(s2, s1, score_cutoff, score_hint);
265
266 remove_common_affix(s1, s2);
267 if (s1.empty())
268 return (s2.size() <= score_cutoff) ? s2.size() : score_cutoff + 1;
269 else if (s1.size() < 64)
270 return osa_hyrroe2003(PatternMatchVector(s1), s1, s2, score_cutoff);
271 else
272 return osa_hyrroe2003_block(BlockPatternMatchVector(s1), s1, s2, score_cutoff);
273 }
274};
275
276} // namespace detail
277} // namespace rapidfuzz
Generated by doxygen 1.9.8