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