libstdc++
multiseq_selection.h
Go to the documentation of this file.
1 // -*- C++ -*-
2 
3 // Copyright (C) 2007-2025 Free Software Foundation, Inc.
4 //
5 // This file is part of the GNU ISO C++ Library. This library is free
6 // software; you can redistribute it and/or modify it under the terms
7 // of the GNU General Public License as published by the Free Software
8 // Foundation; either version 3, or (at your option) any later
9 // version.
10 
11 // This library is distributed in the hope that it will be useful, but
12 // WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 // General Public License for more details.
15 
16 // Under Section 7 of GPL version 3, you are granted additional
17 // permissions described in the GCC Runtime Library Exception, version
18 // 3.1, as published by the Free Software Foundation.
19 
20 // You should have received a copy of the GNU General Public License and
21 // a copy of the GCC Runtime Library Exception along with this program;
22 // see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
23 // <http://www.gnu.org/licenses/>.
24 
25 /** @file parallel/multiseq_selection.h
26  * @brief Functions to find elements of a certain global __rank in
27  * multiple sorted sequences. Also serves for splitting such
28  * sequence sets.
29  *
30  * The algorithm description can be found in
31  *
32  * P. J. Varman, S. D. Scheufler, B. R. Iyer, and G. R. Ricard.
33  * Merging Multiple Lists on Hierarchical-Memory Multiprocessors.
34  * Journal of Parallel and Distributed Computing, 12(2):171-177, 1991.
35  *
36  * This file is a GNU parallel extension to the Standard C++ Library.
37  */
38 
39 // Written by Johannes Singler.
40 
41 #ifndef _GLIBCXX_PARALLEL_MULTISEQ_SELECTION_H
42 #define _GLIBCXX_PARALLEL_MULTISEQ_SELECTION_H 1
43 
44 #include <vector>
45 #include <queue>
46 
47 #include <bits/stl_algo.h>
48 
49 namespace __gnu_parallel
50 {
51 
52 #pragma GCC diagnostic push
53 #pragma GCC diagnostic ignored "-Wdeprecated-declarations" // *nary_function
54 
55  /** @brief Compare __a pair of types lexicographically, ascending. */
56  template<typename _T1, typename _T2, typename _Compare>
58  : public std::binary_function<std::pair<_T1, _T2>,
59  std::pair<_T1, _T2>, bool>
60  {
61  private:
62  _Compare& _M_comp;
63 
64  public:
65  _Lexicographic(_Compare& __comp) : _M_comp(__comp) { }
66 
67  bool
68  operator()(const std::pair<_T1, _T2>& __p1,
69  const std::pair<_T1, _T2>& __p2) const
70  {
71  if (_M_comp(__p1.first, __p2.first))
72  return true;
73 
74  if (_M_comp(__p2.first, __p1.first))
75  return false;
76 
77  // Firsts are equal.
78  return __p1.second < __p2.second;
79  }
80  };
81 
82  /** @brief Compare __a pair of types lexicographically, descending. */
83  template<typename _T1, typename _T2, typename _Compare>
84  class _LexicographicReverse : public std::binary_function<_T1, _T2, bool>
85  {
86  private:
87  _Compare& _M_comp;
88 
89  public:
90  _LexicographicReverse(_Compare& __comp) : _M_comp(__comp) { }
91 
92  bool
93  operator()(const std::pair<_T1, _T2>& __p1,
94  const std::pair<_T1, _T2>& __p2) const
95  {
96  if (_M_comp(__p2.first, __p1.first))
97  return true;
98 
99  if (_M_comp(__p1.first, __p2.first))
100  return false;
101 
102  // Firsts are equal.
103  return __p2.second < __p1.second;
104  }
105  };
106 
107 #pragma GCC diagnostic pop // -Wdeprecated-declarations
108 
109  /**
110  * @brief Splits several sorted sequences at a certain global __rank,
111  * resulting in a splitting point for each sequence.
112  * The sequences are passed via a sequence of random-access
113  * iterator pairs, none of the sequences may be empty. If there
114  * are several equal elements across the split, the ones on the
115  * __left side will be chosen from sequences with smaller number.
116  * @param __begin_seqs Begin of the sequence of iterator pairs.
117  * @param __end_seqs End of the sequence of iterator pairs.
118  * @param __rank The global rank to partition at.
119  * @param __begin_offsets A random-access __sequence __begin where the
120  * __result will be stored in. Each element of the sequence is an
121  * iterator that points to the first element on the greater part of
122  * the respective __sequence.
123  * @param __comp The ordering functor, defaults to std::less<_Tp>.
124  */
125  template<typename _RanSeqs, typename _RankType, typename _RankIterator,
126  typename _Compare>
127  void
128  multiseq_partition(_RanSeqs __begin_seqs, _RanSeqs __end_seqs,
129  _RankType __rank,
130  _RankIterator __begin_offsets,
131  _Compare __comp = std::less<
132  typename std::iterator_traits<typename
134  first_type>::value_type>()) // std::less<_Tp>
135  {
136  _GLIBCXX_CALL(__end_seqs - __begin_seqs)
137 
139  _It;
141  _SeqNumber;
143  _DifferenceType;
144  typedef typename std::iterator_traits<_It>::value_type _ValueType;
145 
148 
149  // Number of sequences, number of elements in total (possibly
150  // including padding).
151  _DifferenceType __m = std::distance(__begin_seqs, __end_seqs), __nn = 0,
152  __nmax, __n, __r;
153 
154  for (_SeqNumber __i = 0; __i < __m; __i++)
155  {
156  __nn += std::distance(__begin_seqs[__i].first,
157  __begin_seqs[__i].second);
158  _GLIBCXX_PARALLEL_ASSERT(
159  std::distance(__begin_seqs[__i].first,
160  __begin_seqs[__i].second) > 0);
161  }
162 
163  if (__rank == __nn)
164  {
165  for (_SeqNumber __i = 0; __i < __m; __i++)
166  __begin_offsets[__i] = __begin_seqs[__i].second; // Very end.
167  // Return __m - 1;
168  return;
169  }
170 
171  _GLIBCXX_PARALLEL_ASSERT(__m != 0);
172  _GLIBCXX_PARALLEL_ASSERT(__nn != 0);
173  _GLIBCXX_PARALLEL_ASSERT(__rank >= 0);
174  _GLIBCXX_PARALLEL_ASSERT(__rank < __nn);
175 
176  _DifferenceType* __ns = new _DifferenceType[__m];
177  _DifferenceType* __a = new _DifferenceType[__m];
178  _DifferenceType* __b = new _DifferenceType[__m];
179  _DifferenceType __l;
180 
181  __ns[0] = std::distance(__begin_seqs[0].first, __begin_seqs[0].second);
182  __nmax = __ns[0];
183  for (_SeqNumber __i = 0; __i < __m; __i++)
184  {
185  __ns[__i] = std::distance(__begin_seqs[__i].first,
186  __begin_seqs[__i].second);
187  __nmax = std::max(__nmax, __ns[__i]);
188  }
189 
190  __r = __rd_log2(__nmax) + 1;
191 
192 #pragma GCC diagnostic push
193 #pragma GCC diagnostic ignored "-Wlong-long" // LL literal
194  // Pad all lists to this length, at least as long as any ns[__i],
195  // equality iff __nmax = 2^__k - 1.
196  __l = (1ULL << __r) - 1;
197 #pragma GCC diagnostic pop
198 
199  for (_SeqNumber __i = 0; __i < __m; __i++)
200  {
201  __a[__i] = 0;
202  __b[__i] = __l;
203  }
204  __n = __l / 2;
205 
206  // Invariants:
207  // 0 <= __a[__i] <= __ns[__i], 0 <= __b[__i] <= __l
208 
209 #define __S(__i) (__begin_seqs[__i].first)
210 
211  // Initial partition.
213 
214  for (_SeqNumber __i = 0; __i < __m; __i++)
215  if (__n < __ns[__i]) //__sequence long enough
216  __sample.push_back(std::make_pair(__S(__i)[__n], __i));
217  __gnu_sequential::sort(__sample.begin(), __sample.end(), __lcomp);
218 
219  for (_SeqNumber __i = 0; __i < __m; __i++) //conceptual infinity
220  if (__n >= __ns[__i]) //__sequence too short, conceptual infinity
221  __sample.push_back(
222  std::make_pair(__S(__i)[0] /*__dummy element*/, __i));
223 
224  _DifferenceType __localrank = __rank / __l;
225 
226  _SeqNumber __j;
227  for (__j = 0;
228  __j < __localrank && ((__n + 1) <= __ns[__sample[__j].second]);
229  ++__j)
230  __a[__sample[__j].second] += __n + 1;
231  for (; __j < __m; __j++)
232  __b[__sample[__j].second] -= __n + 1;
233 
234  // Further refinement.
235  while (__n > 0)
236  {
237  __n /= 2;
238 
239  _SeqNumber __lmax_seq = -1; // to avoid warning
240  const _ValueType* __lmax = 0; // impossible to avoid the warning?
241  for (_SeqNumber __i = 0; __i < __m; __i++)
242  {
243  if (__a[__i] > 0)
244  {
245  if (!__lmax)
246  {
247  __lmax = &(__S(__i)[__a[__i] - 1]);
248  __lmax_seq = __i;
249  }
250  else
251  {
252  // Max, favor rear sequences.
253  if (!__comp(__S(__i)[__a[__i] - 1], *__lmax))
254  {
255  __lmax = &(__S(__i)[__a[__i] - 1]);
256  __lmax_seq = __i;
257  }
258  }
259  }
260  }
261 
262  _SeqNumber __i;
263  for (__i = 0; __i < __m; __i++)
264  {
265  _DifferenceType __middle = (__b[__i] + __a[__i]) / 2;
266  if (__lmax && __middle < __ns[__i] &&
267  __lcomp(std::make_pair(__S(__i)[__middle], __i),
268  std::make_pair(*__lmax, __lmax_seq)))
269  __a[__i] = std::min(__a[__i] + __n + 1, __ns[__i]);
270  else
271  __b[__i] -= __n + 1;
272  }
273 
274  _DifferenceType __leftsize = 0;
275  for (_SeqNumber __i = 0; __i < __m; __i++)
276  __leftsize += __a[__i] / (__n + 1);
277 
278  _DifferenceType __skew = __rank / (__n + 1) - __leftsize;
279 
280  if (__skew > 0)
281  {
282  // Move to the left, find smallest.
286  __pq(__lrcomp);
287 
288  for (_SeqNumber __i = 0; __i < __m; __i++)
289  if (__b[__i] < __ns[__i])
290  __pq.push(std::make_pair(__S(__i)[__b[__i]], __i));
291 
292  for (; __skew != 0 && !__pq.empty(); --__skew)
293  {
294  _SeqNumber __source = __pq.top().second;
295  __pq.pop();
296 
297  __a[__source]
298  = std::min(__a[__source] + __n + 1, __ns[__source]);
299  __b[__source] += __n + 1;
300 
301  if (__b[__source] < __ns[__source])
302  __pq.push(
303  std::make_pair(__S(__source)[__b[__source]], __source));
304  }
305  }
306  else if (__skew < 0)
307  {
308  // Move to the right, find greatest.
312  __pq(__lcomp);
313 
314  for (_SeqNumber __i = 0; __i < __m; __i++)
315  if (__a[__i] > 0)
316  __pq.push(std::make_pair(__S(__i)[__a[__i] - 1], __i));
317 
318  for (; __skew != 0; ++__skew)
319  {
320  _SeqNumber __source = __pq.top().second;
321  __pq.pop();
322 
323  __a[__source] -= __n + 1;
324  __b[__source] -= __n + 1;
325 
326  if (__a[__source] > 0)
327  __pq.push(std::make_pair(
328  __S(__source)[__a[__source] - 1], __source));
329  }
330  }
331  }
332 
333  // Postconditions:
334  // __a[__i] == __b[__i] in most cases, except when __a[__i] has been
335  // clamped because of having reached the boundary
336 
337  // Now return the result, calculate the offset.
338 
339  // Compare the keys on both edges of the border.
340 
341  // Maximum of left edge, minimum of right edge.
342  _ValueType* __maxleft = 0;
343  _ValueType* __minright = 0;
344  for (_SeqNumber __i = 0; __i < __m; __i++)
345  {
346  if (__a[__i] > 0)
347  {
348  if (!__maxleft)
349  __maxleft = &(__S(__i)[__a[__i] - 1]);
350  else
351  {
352  // Max, favor rear sequences.
353  if (!__comp(__S(__i)[__a[__i] - 1], *__maxleft))
354  __maxleft = &(__S(__i)[__a[__i] - 1]);
355  }
356  }
357  if (__b[__i] < __ns[__i])
358  {
359  if (!__minright)
360  __minright = &(__S(__i)[__b[__i]]);
361  else
362  {
363  // Min, favor fore sequences.
364  if (__comp(__S(__i)[__b[__i]], *__minright))
365  __minright = &(__S(__i)[__b[__i]]);
366  }
367  }
368  }
369 
370  _SeqNumber __seq = 0;
371  for (_SeqNumber __i = 0; __i < __m; __i++)
372  __begin_offsets[__i] = __S(__i) + __a[__i];
373 
374  delete[] __ns;
375  delete[] __a;
376  delete[] __b;
377  }
378 
379 
380  /**
381  * @brief Selects the element at a certain global __rank from several
382  * sorted sequences.
383  *
384  * The sequences are passed via a sequence of random-access
385  * iterator pairs, none of the sequences may be empty.
386  * @param __begin_seqs Begin of the sequence of iterator pairs.
387  * @param __end_seqs End of the sequence of iterator pairs.
388  * @param __rank The global rank to partition at.
389  * @param __offset The rank of the selected element in the global
390  * subsequence of elements equal to the selected element. If the
391  * selected element is unique, this number is 0.
392  * @param __comp The ordering functor, defaults to std::less.
393  */
394  template<typename _Tp, typename _RanSeqs, typename _RankType,
395  typename _Compare>
396  _Tp
397  multiseq_selection(_RanSeqs __begin_seqs, _RanSeqs __end_seqs,
398  _RankType __rank,
399  _RankType& __offset, _Compare __comp = std::less<_Tp>())
400  {
401  _GLIBCXX_CALL(__end_seqs - __begin_seqs)
402 
404  _It;
406  _SeqNumber;
408  _DifferenceType;
409 
412 
413  // Number of sequences, number of elements in total (possibly
414  // including padding).
415  _DifferenceType __m = std::distance(__begin_seqs, __end_seqs);
416  _DifferenceType __nn = 0;
417  _DifferenceType __nmax, __n, __r;
418 
419  for (_SeqNumber __i = 0; __i < __m; __i++)
420  __nn += std::distance(__begin_seqs[__i].first,
421  __begin_seqs[__i].second);
422 
423  if (__m == 0 || __nn == 0 || __rank < 0 || __rank >= __nn)
424  {
425  // result undefined if there is no data or __rank is outside bounds
426  throw std::exception();
427  }
428 
429 
430  _DifferenceType* __ns = new _DifferenceType[__m];
431  _DifferenceType* __a = new _DifferenceType[__m];
432  _DifferenceType* __b = new _DifferenceType[__m];
433  _DifferenceType __l;
434 
435  __ns[0] = std::distance(__begin_seqs[0].first, __begin_seqs[0].second);
436  __nmax = __ns[0];
437  for (_SeqNumber __i = 0; __i < __m; ++__i)
438  {
439  __ns[__i] = std::distance(__begin_seqs[__i].first,
440  __begin_seqs[__i].second);
441  __nmax = std::max(__nmax, __ns[__i]);
442  }
443 
444  __r = __rd_log2(__nmax) + 1;
445 
446  // Pad all lists to this length, at least as long as any ns[__i],
447  // equality iff __nmax = 2^__k - 1
448  __l = __round_up_to_pow2(__r) - 1;
449 
450  for (_SeqNumber __i = 0; __i < __m; ++__i)
451  {
452  __a[__i] = 0;
453  __b[__i] = __l;
454  }
455  __n = __l / 2;
456 
457  // Invariants:
458  // 0 <= __a[__i] <= __ns[__i], 0 <= __b[__i] <= __l
459 
460 #define __S(__i) (__begin_seqs[__i].first)
461 
462  // Initial partition.
464 
465  for (_SeqNumber __i = 0; __i < __m; __i++)
466  if (__n < __ns[__i])
467  __sample.push_back(std::make_pair(__S(__i)[__n], __i));
468  __gnu_sequential::sort(__sample.begin(), __sample.end(),
469  __lcomp, sequential_tag());
470 
471  // Conceptual infinity.
472  for (_SeqNumber __i = 0; __i < __m; __i++)
473  if (__n >= __ns[__i])
474  __sample.push_back(
475  std::make_pair(__S(__i)[0] /*__dummy element*/, __i));
476 
477  _DifferenceType __localrank = __rank / __l;
478 
479  _SeqNumber __j;
480  for (__j = 0;
481  __j < __localrank && ((__n + 1) <= __ns[__sample[__j].second]);
482  ++__j)
483  __a[__sample[__j].second] += __n + 1;
484  for (; __j < __m; ++__j)
485  __b[__sample[__j].second] -= __n + 1;
486 
487  // Further refinement.
488  while (__n > 0)
489  {
490  __n /= 2;
491 
492  const _Tp* __lmax = 0;
493  for (_SeqNumber __i = 0; __i < __m; ++__i)
494  {
495  if (__a[__i] > 0)
496  {
497  if (!__lmax)
498  __lmax = &(__S(__i)[__a[__i] - 1]);
499  else
500  {
501  if (__comp(*__lmax, __S(__i)[__a[__i] - 1])) //max
502  __lmax = &(__S(__i)[__a[__i] - 1]);
503  }
504  }
505  }
506 
507  _SeqNumber __i;
508  for (__i = 0; __i < __m; __i++)
509  {
510  _DifferenceType __middle = (__b[__i] + __a[__i]) / 2;
511  if (__lmax && __middle < __ns[__i]
512  && __comp(__S(__i)[__middle], *__lmax))
513  __a[__i] = std::min(__a[__i] + __n + 1, __ns[__i]);
514  else
515  __b[__i] -= __n + 1;
516  }
517 
518  _DifferenceType __leftsize = 0;
519  for (_SeqNumber __i = 0; __i < __m; ++__i)
520  __leftsize += __a[__i] / (__n + 1);
521 
522  _DifferenceType __skew = __rank / (__n + 1) - __leftsize;
523 
524  if (__skew > 0)
525  {
526  // Move to the left, find smallest.
530  __pq(__lrcomp);
531 
532  for (_SeqNumber __i = 0; __i < __m; ++__i)
533  if (__b[__i] < __ns[__i])
534  __pq.push(std::make_pair(__S(__i)[__b[__i]], __i));
535 
536  for (; __skew != 0 && !__pq.empty(); --__skew)
537  {
538  _SeqNumber __source = __pq.top().second;
539  __pq.pop();
540 
541  __a[__source]
542  = std::min(__a[__source] + __n + 1, __ns[__source]);
543  __b[__source] += __n + 1;
544 
545  if (__b[__source] < __ns[__source])
546  __pq.push(
547  std::make_pair(__S(__source)[__b[__source]], __source));
548  }
549  }
550  else if (__skew < 0)
551  {
552  // Move to the right, find greatest.
556 
557  for (_SeqNumber __i = 0; __i < __m; ++__i)
558  if (__a[__i] > 0)
559  __pq.push(std::make_pair(__S(__i)[__a[__i] - 1], __i));
560 
561  for (; __skew != 0; ++__skew)
562  {
563  _SeqNumber __source = __pq.top().second;
564  __pq.pop();
565 
566  __a[__source] -= __n + 1;
567  __b[__source] -= __n + 1;
568 
569  if (__a[__source] > 0)
570  __pq.push(std::make_pair(
571  __S(__source)[__a[__source] - 1], __source));
572  }
573  }
574  }
575 
576  // Postconditions:
577  // __a[__i] == __b[__i] in most cases, except when __a[__i] has been
578  // clamped because of having reached the boundary
579 
580  // Now return the result, calculate the offset.
581 
582  // Compare the keys on both edges of the border.
583 
584  // Maximum of left edge, minimum of right edge.
585  bool __maxleftset = false, __minrightset = false;
586 
587  // Impossible to avoid the warning?
588  _Tp __maxleft, __minright;
589  for (_SeqNumber __i = 0; __i < __m; ++__i)
590  {
591  if (__a[__i] > 0)
592  {
593  if (!__maxleftset)
594  {
595  __maxleft = __S(__i)[__a[__i] - 1];
596  __maxleftset = true;
597  }
598  else
599  {
600  // Max.
601  if (__comp(__maxleft, __S(__i)[__a[__i] - 1]))
602  __maxleft = __S(__i)[__a[__i] - 1];
603  }
604  }
605  if (__b[__i] < __ns[__i])
606  {
607  if (!__minrightset)
608  {
609  __minright = __S(__i)[__b[__i]];
610  __minrightset = true;
611  }
612  else
613  {
614  // Min.
615  if (__comp(__S(__i)[__b[__i]], __minright))
616  __minright = __S(__i)[__b[__i]];
617  }
618  }
619  }
620 
621  // Minright is the __splitter, in any case.
622 
623  if (!__maxleftset || __comp(__minright, __maxleft))
624  {
625  // Good luck, everything is split unambiguously.
626  __offset = 0;
627  }
628  else
629  {
630  // We have to calculate an offset.
631  __offset = 0;
632 
633  for (_SeqNumber __i = 0; __i < __m; ++__i)
634  {
635  _DifferenceType lb
636  = std::lower_bound(__S(__i), __S(__i) + __ns[__i],
637  __minright,
638  __comp) - __S(__i);
639  __offset += __a[__i] - lb;
640  }
641  }
642 
643  delete[] __ns;
644  delete[] __a;
645  delete[] __b;
646 
647  return __minright;
648  }
649 }
650 
651 #undef __S
652 
653 #endif /* _GLIBCXX_PARALLEL_MULTISEQ_SELECTION_H */
#define _GLIBCXX_CALL(__n)
Macro to produce log message when entering a function.
constexpr const _Tp & max(const _Tp &, const _Tp &)
This does what you think it does.
Definition: stl_algobase.h:258
constexpr const _Tp & min(const _Tp &, const _Tp &)
This does what you think it does.
Definition: stl_algobase.h:234
_RandomAccessIterator __sample(_InputIterator __first, _InputIterator __last, input_iterator_tag, _RandomAccessIterator __out, random_access_iterator_tag, _Size __n, _UniformRandomBitGenerator &&__g)
Reservoir sampling algorithm.
Definition: stl_algo.h:5807
constexpr iterator_traits< _InputIterator >::difference_type distance(_InputIterator __first, _InputIterator __last)
A generalization of pointer arithmetic.
GNU parallel code for public use.
_Tp multiseq_selection(_RanSeqs __begin_seqs, _RanSeqs __end_seqs, _RankType __rank, _RankType &__offset, _Compare __comp=std::less< _Tp >())
Selects the element at a certain global __rank from several sorted sequences.
_Tp __round_up_to_pow2(_Tp __x)
Round up to the next greater power of 2.
void multiseq_partition(_RanSeqs __begin_seqs, _RanSeqs __end_seqs, _RankType __rank, _RankIterator __begin_offsets, _Compare __comp=std::less< typename std::iterator_traits< typename std::iterator_traits< _RanSeqs >::value_type::first_type >::value_type >())
Splits several sorted sequences at a certain global __rank, resulting in a splitting point for each s...
_Size __rd_log2(_Size __n)
Calculates the rounded-down logarithm of __n for base 2.
Definition: base.h:102
Base class for all library exceptions.
Definition: exception.h:62
One of the comparison functors.
Definition: stl_function.h:401
Struct holding two objects of arbitrary type.
Definition: stl_pair.h:304
_T1 first
The first member.
Definition: stl_pair.h:308
_T2 second
The second member.
Definition: stl_pair.h:309
Traits class for iterators.
A standard container automatically sorting its contents.
Definition: stl_queue.h:553
bool empty() const
Definition: stl_queue.h:808
void pop()
Removes first element.
Definition: stl_queue.h:886
const_reference top() const
Definition: stl_queue.h:823
void push(const value_type &__x)
Add data to the queue.
Definition: stl_queue.h:838
A standard container which offers fixed time access to individual elements in any order.
Definition: stl_vector.h:459
Compare __a pair of types lexicographically, ascending.
Compare __a pair of types lexicographically, descending.
Forces sequential execution at compile time.
Definition: tags.h:42