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