libstdc++
numeric
Go to the documentation of this file.
1 // <numeric> -*- C++ -*-
2 
3 // Copyright (C) 2001-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
7 // terms of the GNU General Public License as published by the
8 // Free Software Foundation; either version 3, or (at your option)
9 // any later version.
10 
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 // GNU 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 /*
26  *
27  * Copyright (c) 1994
28  * Hewlett-Packard Company
29  *
30  * Permission to use, copy, modify, distribute and sell this software
31  * and its documentation for any purpose is hereby granted without fee,
32  * provided that the above copyright notice appear in all copies and
33  * that both that copyright notice and this permission notice appear
34  * in supporting documentation. Hewlett-Packard Company makes no
35  * representations about the suitability of this software for any
36  * purpose. It is provided "as is" without express or implied warranty.
37  *
38  *
39  * Copyright (c) 1996,1997
40  * Silicon Graphics Computer Systems, Inc.
41  *
42  * Permission to use, copy, modify, distribute and sell this software
43  * and its documentation for any purpose is hereby granted without fee,
44  * provided that the above copyright notice appear in all copies and
45  * that both that copyright notice and this permission notice appear
46  * in supporting documentation. Silicon Graphics makes no
47  * representations about the suitability of this software for any
48  * purpose. It is provided "as is" without express or implied warranty.
49  */
50 
51 /** @file include/numeric
52  * This is a Standard C++ Library header.
53  */
54 
55 #ifndef _GLIBCXX_NUMERIC
56 #define _GLIBCXX_NUMERIC 1
57 
58 #ifdef _GLIBCXX_SYSHDR
59 #pragma GCC system_header
60 #endif
61 
62 #include <bits/c++config.h>
63 #include <bits/stl_iterator_base_types.h>
64 #include <bits/stl_numeric.h>
65 
66 #ifdef _GLIBCXX_PARALLEL
67 # include <parallel/numeric>
68 #endif
69 
70 #if __cplusplus >= 201402L
71 # include <type_traits>
72 # include <bit>
73 # include <ext/numeric_traits.h>
74 #endif
75 
76 #if __cplusplus >= 201703L
77 # include <bits/stl_function.h>
78 #endif
79 
80 #if __cplusplus > 201703L
81 # include <limits>
82 #endif
83 
84 #define __glibcxx_want_constexpr_numeric
85 #define __glibcxx_want_gcd
86 #define __glibcxx_want_gcd_lcm
87 #define __glibcxx_want_interpolate
88 #define __glibcxx_want_lcm
89 #define __glibcxx_want_parallel_algorithm
90 #define __glibcxx_want_ranges_iota
91 #define __glibcxx_want_saturation_arithmetic
92 #include <bits/version.h>
93 
94 #if __glibcxx_ranges_iota >= 202202L // C++ >= 23
95 # include <bits/ranges_algobase.h> // for ranges::out_value_result
96 #endif
97 
98 #ifdef __glibcxx_saturation_arithmetic // C++ >= 26
99 # include <bits/sat_arith.h>
100 #endif
101 
102 /**
103  * @defgroup numerics Numerics
104  *
105  * Components for performing numeric operations. Includes support for
106  * complex number types, random number generation, numeric (n-at-a-time)
107  * arrays, generalized numeric algorithms, and mathematical special functions.
108  */
109 
110 namespace std _GLIBCXX_VISIBILITY(default)
111 {
112 _GLIBCXX_BEGIN_NAMESPACE_VERSION
113 
114 #if __cplusplus >= 201402L
115 namespace __detail
116 {
117  // Like std::abs, but supports unsigned types and returns the specified type,
118  // so |std::numeric_limits<_Tp>::min()| is OK if representable in _Res.
119  template<typename _Res, typename _Tp>
120  constexpr _Res
121  __abs_r(_Tp __val)
122  {
123  static_assert(sizeof(_Res) >= sizeof(_Tp),
124  "result type must be at least as wide as the input type");
125 
126  if (__val >= 0)
127  return __val;
128 #ifdef _GLIBCXX_ASSERTIONS
129  if (!__is_constant_evaluated()) // overflow already detected in constexpr
130  __glibcxx_assert(__val != __gnu_cxx::__int_traits<_Res>::__min);
131 #endif
132  return -static_cast<_Res>(__val);
133  }
134 
135  template<typename> void __abs_r(bool) = delete;
136 
137  // GCD implementation, using Stein's algorithm
138  template<typename _Tp>
139  constexpr _Tp
140  __gcd(_Tp __m, _Tp __n)
141  {
142  static_assert(is_unsigned<_Tp>::value, "type must be unsigned");
143 
144  if (__m == 0)
145  return __n;
146  if (__n == 0)
147  return __m;
148 
149  const int __i = std::__countr_zero(__m);
150  __m >>= __i;
151  const int __j = std::__countr_zero(__n);
152  __n >>= __j;
153  const int __k = __i < __j ? __i : __j; // min(i, j)
154 
155  while (true)
156  {
157  if (__m > __n)
158  {
159  _Tp __tmp = __m;
160  __m = __n;
161  __n = __tmp;
162  }
163 
164  __n -= __m;
165 
166  if (__n == 0)
167  return __m << __k;
168 
169  __n >>= std::__countr_zero(__n);
170  }
171  }
172 } // namespace __detail
173 #endif // C++14
174 
175 #ifdef __cpp_lib_gcd_lcm // C++ >= 17
176  /// Greatest common divisor
177  template<typename _Mn, typename _Nn>
178  constexpr common_type_t<_Mn, _Nn>
179  gcd(_Mn __m, _Nn __n) noexcept
180  {
181  static_assert(is_integral_v<_Mn> && is_integral_v<_Nn>,
182  "std::gcd arguments must be integers");
183  static_assert(_Mn(2) == 2 && _Nn(2) == 2,
184  "std::gcd arguments must not be bool");
185  using _Ct = common_type_t<_Mn, _Nn>;
186  const _Ct __m2 = __detail::__abs_r<_Ct>(__m);
187  const _Ct __n2 = __detail::__abs_r<_Ct>(__n);
188  return __detail::__gcd<make_unsigned_t<_Ct>>(__m2, __n2);
189  }
190 
191  /// Least common multiple
192  template<typename _Mn, typename _Nn>
193  constexpr common_type_t<_Mn, _Nn>
194  lcm(_Mn __m, _Nn __n) noexcept
195  {
196  static_assert(is_integral_v<_Mn> && is_integral_v<_Nn>,
197  "std::lcm arguments must be integers");
198  static_assert(_Mn(2) == 2 && _Nn(2) == 2,
199  "std::lcm arguments must not be bool");
200  using _Ct = common_type_t<_Mn, _Nn>;
201  const _Ct __m2 = __detail::__abs_r<_Ct>(__m);
202  const _Ct __n2 = __detail::__abs_r<_Ct>(__n);
203  if (__m2 == 0 || __n2 == 0)
204  return 0;
205  _Ct __r = __m2 / __detail::__gcd<make_unsigned_t<_Ct>>(__m2, __n2);
206 
207  if constexpr (is_signed_v<_Ct>)
208  if (__is_constant_evaluated())
209  return __r * __n2; // constant evaluation can detect overflow here.
210 
211  bool __overflow = __builtin_mul_overflow(__r, __n2, &__r);
212  __glibcxx_assert(!__overflow);
213  return __r;
214  }
215 
216 #endif // __cpp_lib_gcd_lcm
217 
218  // midpoint
219 #ifdef __cpp_lib_interpolate // C++ >= 20
220  template<typename _Tp>
221  constexpr
222  enable_if_t<__and_v<is_arithmetic<_Tp>, is_same<remove_cv_t<_Tp>, _Tp>,
223  __not_<is_same<_Tp, bool>>>,
224  _Tp>
225  midpoint(_Tp __a, _Tp __b) noexcept
226  {
227  if constexpr (is_integral_v<_Tp>)
228  {
229  using _Up = make_unsigned_t<_Tp>;
230 
231  int __k = 1;
232  _Up __m = __a;
233  _Up __M = __b;
234  if (__a > __b)
235  {
236  __k = -1;
237  __m = __b;
238  __M = __a;
239  }
240  return __a + __k * _Tp(_Up(__M - __m) / 2);
241  }
242  else // is_floating
243  {
244  constexpr _Tp __lo = numeric_limits<_Tp>::min() * 2;
245  constexpr _Tp __hi = numeric_limits<_Tp>::max() / 2;
246  const _Tp __abs_a = __a < 0 ? -__a : __a;
247  const _Tp __abs_b = __b < 0 ? -__b : __b;
248  if (__abs_a <= __hi && __abs_b <= __hi) [[likely]]
249  return (__a + __b) / 2; // always correctly rounded
250  if (__abs_a < __lo) // not safe to halve __a
251  return __a + __b/2;
252  if (__abs_b < __lo) // not safe to halve __b
253  return __a/2 + __b;
254  return __a/2 + __b/2; // otherwise correctly rounded
255  }
256  }
257 
258  template<typename _Tp>
259  constexpr enable_if_t<is_object_v<_Tp>, _Tp*>
260  midpoint(_Tp* __a, _Tp* __b) noexcept
261  {
262  static_assert( sizeof(_Tp) != 0, "type must be complete" );
263  return __a + (__b - __a) / 2;
264  }
265 #endif // __cpp_lib_interpolate
266 
267 #if __cplusplus >= 201703L
268  /// @addtogroup numeric_ops
269  /// @{
270 
271  /**
272  * @brief Calculate reduction of values in a range.
273  *
274  * @param __first Start of range.
275  * @param __last End of range.
276  * @param __init Starting value to add other values to.
277  * @param __binary_op A binary function object.
278  * @return The final sum.
279  *
280  * Reduce the values in the range `[first,last)` using a binary operation.
281  * The initial value is `init`. The values are not necessarily processed
282  * in order.
283  *
284  * This algorithm is similar to `std::accumulate` but is not required to
285  * perform the operations in order from first to last. For operations
286  * that are commutative and associative the result will be the same as
287  * for `std::accumulate`, but for other operations (such as floating point
288  * arithmetic) the result can be different.
289  */
290  template<typename _InputIterator, typename _Tp, typename _BinaryOperation>
291  _GLIBCXX20_CONSTEXPR
292  _Tp
293  reduce(_InputIterator __first, _InputIterator __last, _Tp __init,
294  _BinaryOperation __binary_op)
295  {
296  using __ref = typename iterator_traits<_InputIterator>::reference;
297  static_assert(is_invocable_r_v<_Tp, _BinaryOperation&, _Tp&, __ref>);
298  static_assert(is_invocable_r_v<_Tp, _BinaryOperation&, __ref, _Tp&>);
299  static_assert(is_invocable_r_v<_Tp, _BinaryOperation&, _Tp&, _Tp&>);
300  static_assert(is_invocable_r_v<_Tp, _BinaryOperation&, __ref, __ref>);
301  if constexpr (__is_random_access_iter<_InputIterator>::value)
302  {
303  while ((__last - __first) >= 4)
304  {
305  _Tp __v1 = __binary_op(__first[0], __first[1]);
306  _Tp __v2 = __binary_op(__first[2], __first[3]);
307  _Tp __v3 = __binary_op(__v1, __v2);
308  __init = __binary_op(__init, __v3);
309  __first += 4;
310  }
311  }
312  for (; __first != __last; ++__first)
313  __init = __binary_op(__init, *__first);
314  return __init;
315  }
316 
317  /**
318  * @brief Calculate reduction of values in a range.
319  *
320  * @param __first Start of range.
321  * @param __last End of range.
322  * @param __init Starting value to add other values to.
323  * @return The final sum.
324  *
325  * Reduce the values in the range `[first,last)` using addition.
326  * Equivalent to calling `std::reduce(first, last, init, std::plus<>())`.
327  */
328  template<typename _InputIterator, typename _Tp>
329  _GLIBCXX20_CONSTEXPR
330  inline _Tp
331  reduce(_InputIterator __first, _InputIterator __last, _Tp __init)
332  { return std::reduce(__first, __last, std::move(__init), plus<>()); }
333 
334  /**
335  * @brief Calculate reduction of values in a range.
336  *
337  * @param __first Start of range.
338  * @param __last End of range.
339  * @return The final sum.
340  *
341  * Reduce the values in the range `[first,last)` using addition, with
342  * an initial value of `T{}`, where `T` is the iterator's value type.
343  * Equivalent to calling `std::reduce(first, last, T{}, std::plus<>())`.
344  */
345  template<typename _InputIterator>
346  _GLIBCXX20_CONSTEXPR
347  inline typename iterator_traits<_InputIterator>::value_type
348  reduce(_InputIterator __first, _InputIterator __last)
349  {
350  using value_type = typename iterator_traits<_InputIterator>::value_type;
351  return std::reduce(__first, __last, value_type{}, plus<>());
352  }
353 
354  /**
355  * @brief Combine elements from two ranges and reduce
356  *
357  * @param __first1 Start of first range.
358  * @param __last1 End of first range.
359  * @param __first2 Start of second range.
360  * @param __init Starting value to add other values to.
361  * @param __binary_op1 The function used to perform reduction.
362  * @param __binary_op2 The function used to combine values from the ranges.
363  * @return The final sum.
364  *
365  * Call `binary_op2(first1[n],first2[n])` for each `n` in `[0,last1-first1)`
366  * and then use `binary_op1` to reduce the values returned by `binary_op2`
367  * to a single value of type `T`.
368  *
369  * The range beginning at `first2` must contain at least `last1-first1`
370  * elements.
371  */
372  template<typename _InputIterator1, typename _InputIterator2, typename _Tp,
373  typename _BinaryOperation1, typename _BinaryOperation2>
374  _GLIBCXX20_CONSTEXPR
375  _Tp
376  transform_reduce(_InputIterator1 __first1, _InputIterator1 __last1,
377  _InputIterator2 __first2, _Tp __init,
378  _BinaryOperation1 __binary_op1,
379  _BinaryOperation2 __binary_op2)
380  {
381  if constexpr (__and_v<__is_random_access_iter<_InputIterator1>,
382  __is_random_access_iter<_InputIterator2>>)
383  {
384  while ((__last1 - __first1) >= 4)
385  {
386  _Tp __v1 = __binary_op1(__binary_op2(__first1[0], __first2[0]),
387  __binary_op2(__first1[1], __first2[1]));
388  _Tp __v2 = __binary_op1(__binary_op2(__first1[2], __first2[2]),
389  __binary_op2(__first1[3], __first2[3]));
390  _Tp __v3 = __binary_op1(__v1, __v2);
391  __init = __binary_op1(__init, __v3);
392  __first1 += 4;
393  __first2 += 4;
394  }
395  }
396  for (; __first1 != __last1; ++__first1, (void) ++__first2)
397  __init = __binary_op1(__init, __binary_op2(*__first1, *__first2));
398  return __init;
399  }
400 
401  /**
402  * @brief Combine elements from two ranges and reduce
403  *
404  * @param __first1 Start of first range.
405  * @param __last1 End of first range.
406  * @param __first2 Start of second range.
407  * @param __init Starting value to add other values to.
408  * @return The final sum.
409  *
410  * Call `first1[n]*first2[n]` for each `n` in `[0,last1-first1)` and then
411  * use addition to sum those products to a single value of type `T`.
412  *
413  * The range beginning at `first2` must contain at least `last1-first1`
414  * elements.
415  */
416  template<typename _InputIterator1, typename _InputIterator2, typename _Tp>
417  _GLIBCXX20_CONSTEXPR
418  inline _Tp
419  transform_reduce(_InputIterator1 __first1, _InputIterator1 __last1,
420  _InputIterator2 __first2, _Tp __init)
421  {
422  return std::transform_reduce(__first1, __last1, __first2,
423  std::move(__init),
424  plus<>(), multiplies<>());
425  }
426 
427  /**
428  * @brief Transform the elements of a range and reduce
429  *
430  * @param __first Start of range.
431  * @param __last End of range.
432  * @param __init Starting value to add other values to.
433  * @param __binary_op The function used to perform reduction.
434  * @param __unary_op The function used to transform values from the range.
435  * @return The final sum.
436  *
437  * Call `unary_op(first[n])` for each `n` in `[0,last-first)` and then
438  * use `binary_op` to reduce the values returned by `unary_op`
439  * to a single value of type `T`.
440  */
441  template<typename _InputIterator, typename _Tp,
442  typename _BinaryOperation, typename _UnaryOperation>
443  _GLIBCXX20_CONSTEXPR
444  _Tp
445  transform_reduce(_InputIterator __first, _InputIterator __last, _Tp __init,
446  _BinaryOperation __binary_op, _UnaryOperation __unary_op)
447  {
448  if constexpr (__is_random_access_iter<_InputIterator>::value)
449  {
450  while ((__last - __first) >= 4)
451  {
452  _Tp __v1 = __binary_op(__unary_op(__first[0]),
453  __unary_op(__first[1]));
454  _Tp __v2 = __binary_op(__unary_op(__first[2]),
455  __unary_op(__first[3]));
456  _Tp __v3 = __binary_op(__v1, __v2);
457  __init = __binary_op(__init, __v3);
458  __first += 4;
459  }
460  }
461  for (; __first != __last; ++__first)
462  __init = __binary_op(__init, __unary_op(*__first));
463  return __init;
464  }
465 
466  /** @brief Output the cumulative sum of one range to a second range
467  *
468  * @param __first Start of input range.
469  * @param __last End of input range.
470  * @param __result Start of output range.
471  * @param __init Initial value.
472  * @param __binary_op Function to perform summation.
473  * @return The end of the output range.
474  *
475  * Write the cumulative sum (aka prefix sum, aka scan) of the input range
476  * to the output range. Each element of the output range contains the
477  * running total of all earlier elements (and the initial value),
478  * using `binary_op` for summation.
479  *
480  * This function generates an "exclusive" scan, meaning the Nth element
481  * of the output range is the sum of the first N-1 input elements,
482  * so the Nth input element is not included.
483  */
484  template<typename _InputIterator, typename _OutputIterator, typename _Tp,
485  typename _BinaryOperation>
486  _GLIBCXX20_CONSTEXPR
487  _OutputIterator
488  exclusive_scan(_InputIterator __first, _InputIterator __last,
489  _OutputIterator __result, _Tp __init,
490  _BinaryOperation __binary_op)
491  {
492  while (__first != __last)
493  {
494  _Tp __v = std::move(__init);
495  __init = __binary_op(__v, *__first);
496  ++__first;
497  *__result++ = std::move(__v);
498  }
499  return __result;
500  }
501 
502  /** @brief Output the cumulative sum of one range to a second range
503  *
504  * @param __first Start of input range.
505  * @param __last End of input range.
506  * @param __result Start of output range.
507  * @param __init Initial value.
508  * @return The end of the output range.
509  *
510  * Write the cumulative sum (aka prefix sum, aka scan) of the input range
511  * to the output range. Each element of the output range contains the
512  * running total of all earlier elements (and the initial value),
513  * using `std::plus<>` for summation.
514  *
515  * This function generates an "exclusive" scan, meaning the Nth element
516  * of the output range is the sum of the first N-1 input elements,
517  * so the Nth input element is not included.
518  */
519  template<typename _InputIterator, typename _OutputIterator, typename _Tp>
520  _GLIBCXX20_CONSTEXPR
521  inline _OutputIterator
522  exclusive_scan(_InputIterator __first, _InputIterator __last,
523  _OutputIterator __result, _Tp __init)
524  {
525  return std::exclusive_scan(__first, __last, __result, std::move(__init),
526  plus<>());
527  }
528 
529  /** @brief Output the cumulative sum of one range to a second range
530  *
531  * @param __first Start of input range.
532  * @param __last End of input range.
533  * @param __result Start of output range.
534  * @param __binary_op Function to perform summation.
535  * @param __init Initial value.
536  * @return The end of the output range.
537  *
538  * Write the cumulative sum (aka prefix sum, aka scan) of the input range
539  * to the output range. Each element of the output range contains the
540  * running total of all earlier elements (and the initial value),
541  * using `binary_op` for summation.
542  *
543  * This function generates an "inclusive" scan, meaning the Nth element
544  * of the output range is the sum of the first N input elements,
545  * so the Nth input element is included.
546  */
547  template<typename _InputIterator, typename _OutputIterator,
548  typename _BinaryOperation, typename _Tp>
549  _GLIBCXX20_CONSTEXPR
550  _OutputIterator
551  inclusive_scan(_InputIterator __first, _InputIterator __last,
552  _OutputIterator __result, _BinaryOperation __binary_op,
553  _Tp __init)
554  {
555  for (; __first != __last; ++__first)
556  *__result++ = __init = __binary_op(__init, *__first);
557  return __result;
558  }
559 
560  /** @brief Output the cumulative sum of one range to a second range
561  *
562  * @param __first Start of input range.
563  * @param __last End of input range.
564  * @param __result Start of output range.
565  * @param __binary_op Function to perform summation.
566  * @return The end of the output range.
567  *
568  * Write the cumulative sum (aka prefix sum, aka scan) of the input range
569  * to the output range. Each element of the output range contains the
570  * running total of all earlier elements, using `binary_op` for summation.
571  *
572  * This function generates an "inclusive" scan, meaning the Nth element
573  * of the output range is the sum of the first N input elements,
574  * so the Nth input element is included.
575  */
576  template<typename _InputIterator, typename _OutputIterator,
577  typename _BinaryOperation>
578  _GLIBCXX20_CONSTEXPR
579  _OutputIterator
580  inclusive_scan(_InputIterator __first, _InputIterator __last,
581  _OutputIterator __result, _BinaryOperation __binary_op)
582  {
583  if (__first != __last)
584  {
585  auto __init = *__first;
586  *__result++ = __init;
587  ++__first;
588  if (__first != __last)
589  __result = std::inclusive_scan(__first, __last, __result,
590  __binary_op, std::move(__init));
591  }
592  return __result;
593  }
594 
595  /** @brief Output the cumulative sum of one range to a second range
596  *
597  * @param __first Start of input range.
598  * @param __last End of input range.
599  * @param __result Start of output range.
600  * @return The end of the output range.
601  *
602  * Write the cumulative sum (aka prefix sum, aka scan) of the input range
603  * to the output range. Each element of the output range contains the
604  * running total of all earlier elements, using `std::plus<>` for summation.
605  *
606  * This function generates an "inclusive" scan, meaning the Nth element
607  * of the output range is the sum of the first N input elements,
608  * so the Nth input element is included.
609  */
610  template<typename _InputIterator, typename _OutputIterator>
611  _GLIBCXX20_CONSTEXPR
612  inline _OutputIterator
613  inclusive_scan(_InputIterator __first, _InputIterator __last,
614  _OutputIterator __result)
615  { return std::inclusive_scan(__first, __last, __result, plus<>()); }
616 
617  /** @brief Output the cumulative sum of one range to a second range
618  *
619  * @param __first Start of input range.
620  * @param __last End of input range.
621  * @param __result Start of output range.
622  * @param __init Initial value.
623  * @param __binary_op Function to perform summation.
624  * @param __unary_op Function to transform elements of the input range.
625  * @return The end of the output range.
626  *
627  * Write the cumulative sum (aka prefix sum, aka scan) of the input range
628  * to the output range. Each element of the output range contains the
629  * running total of all earlier elements (and the initial value),
630  * using `__unary_op` to transform the input elements
631  * and using `__binary_op` for summation.
632  *
633  * This function generates an "exclusive" scan, meaning the Nth element
634  * of the output range is the sum of the first N-1 input elements,
635  * so the Nth input element is not included.
636  */
637  template<typename _InputIterator, typename _OutputIterator, typename _Tp,
638  typename _BinaryOperation, typename _UnaryOperation>
639  _GLIBCXX20_CONSTEXPR
640  _OutputIterator
641  transform_exclusive_scan(_InputIterator __first, _InputIterator __last,
642  _OutputIterator __result, _Tp __init,
643  _BinaryOperation __binary_op,
644  _UnaryOperation __unary_op)
645  {
646  while (__first != __last)
647  {
648  auto __v = __init;
649  __init = __binary_op(__init, __unary_op(*__first));
650  ++__first;
651  *__result++ = std::move(__v);
652  }
653  return __result;
654  }
655 
656  /** @brief Output the cumulative sum of one range to a second range
657  *
658  * @param __first Start of input range.
659  * @param __last End of input range.
660  * @param __result Start of output range.
661  * @param __binary_op Function to perform summation.
662  * @param __unary_op Function to transform elements of the input range.
663  * @param __init Initial value.
664  * @return The end of the output range.
665  *
666  * Write the cumulative sum (aka prefix sum, aka scan) of the input range
667  * to the output range. Each element of the output range contains the
668  * running total of all earlier elements (and the initial value),
669  * using `__unary_op` to transform the input elements
670  * and using `__binary_op` for summation.
671  *
672  * This function generates an "inclusive" scan, meaning the Nth element
673  * of the output range is the sum of the first N input elements,
674  * so the Nth input element is included.
675  */
676  template<typename _InputIterator, typename _OutputIterator,
677  typename _BinaryOperation, typename _UnaryOperation, typename _Tp>
678  _GLIBCXX20_CONSTEXPR
679  _OutputIterator
680  transform_inclusive_scan(_InputIterator __first, _InputIterator __last,
681  _OutputIterator __result,
682  _BinaryOperation __binary_op,
683  _UnaryOperation __unary_op,
684  _Tp __init)
685  {
686  for (; __first != __last; ++__first)
687  *__result++ = __init = __binary_op(__init, __unary_op(*__first));
688  return __result;
689  }
690 
691  /** @brief Output the cumulative sum of one range to a second range
692  *
693  * @param __first Start of input range.
694  * @param __last End of input range.
695  * @param __result Start of output range.
696  * @param __binary_op Function to perform summation.
697  * @param __unary_op Function to transform elements of the input range.
698  * @return The end of the output range.
699  *
700  * Write the cumulative sum (aka prefix sum, aka scan) of the input range
701  * to the output range. Each element of the output range contains the
702  * running total of all earlier elements,
703  * using `__unary_op` to transform the input elements
704  * and using `__binary_op` for summation.
705  *
706  * This function generates an "inclusive" scan, meaning the Nth element
707  * of the output range is the sum of the first N input elements,
708  * so the Nth input element is included.
709  */
710  template<typename _InputIterator, typename _OutputIterator,
711  typename _BinaryOperation, typename _UnaryOperation>
712  _GLIBCXX20_CONSTEXPR
713  _OutputIterator
714  transform_inclusive_scan(_InputIterator __first, _InputIterator __last,
715  _OutputIterator __result,
716  _BinaryOperation __binary_op,
717  _UnaryOperation __unary_op)
718  {
719  if (__first != __last)
720  {
721  auto __init = __unary_op(*__first);
722  *__result++ = __init;
723  ++__first;
724  if (__first != __last)
725  __result = std::transform_inclusive_scan(__first, __last, __result,
726  __binary_op, __unary_op,
727  std::move(__init));
728  }
729  return __result;
730  }
731 
732  /// @} group numeric_ops
733 #endif // C++17
734 
735 #if __glibcxx_ranges_iota >= 202202L // C++ >= 23
736 namespace ranges
737 {
738  template<typename _Out, typename _Tp>
739  using iota_result = out_value_result<_Out, _Tp>;
740 
741  struct __iota_fn
742  {
743  template<input_or_output_iterator _Out, sentinel_for<_Out> _Sent, weakly_incrementable _Tp>
744  requires indirectly_writable<_Out, const _Tp&>
745  constexpr iota_result<_Out, _Tp>
746  operator()(_Out __first, _Sent __last, _Tp __value) const
747  {
748  while (__first != __last)
749  {
750  *__first = static_cast<const _Tp&>(__value);
751  ++__first;
752  ++__value;
753  }
754  return {std::move(__first), std::move(__value)};
755  }
756 
757  template<weakly_incrementable _Tp, output_range<const _Tp&> _Range>
758  constexpr iota_result<borrowed_iterator_t<_Range>, _Tp>
759  operator()(_Range&& __r, _Tp __value) const
760  { return (*this)(ranges::begin(__r), ranges::end(__r), std::move(__value)); }
761  };
762 
763  inline constexpr __iota_fn iota{};
764 } // namespace ranges
765 #endif // __glibcxx_ranges_iota
766 
767 _GLIBCXX_END_NAMESPACE_VERSION
768 } // namespace std
769 
770 #if __cplusplus >= 201703L && _GLIBCXX_HOSTED
771 // Parallel STL algorithms
772 # if _PSTL_EXECUTION_POLICIES_DEFINED
773 // If <execution> has already been included, pull in implementations
774 # include <pstl/glue_numeric_impl.h>
775 # else
776 // Otherwise just pull in forward declarations
777 # include <pstl/glue_numeric_defs.h>
778 # define _PSTL_NUMERIC_FORWARD_DECLARED 1
779 # endif
780 #endif // C++17
781 
782 #endif /* _GLIBCXX_NUMERIC */