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
90namespace std _GLIBCXX_VISIBILITY(default)
91{
92_GLIBCXX_BEGIN_NAMESPACE_VERSION
93
94#if __cplusplus >= 201402L
95namespace __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 */
750

source code of include/c++/11/numeric