1// <numeric> -*- C++ -*-
2
3// Copyright (C) 2001-2022 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#ifdef _GLIBCXX_ASSERTIONS
109 if (!__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 201606L
157// These were used in drafts of SD-6:
158#define __cpp_lib_gcd 201606L
159#define __cpp_lib_lcm 201606L
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 constexpr (is_signed_v<_Ct>)
193 if (__is_constant_evaluated())
194 return __r * __n2; // constant evaluation can detect overflow here.
195
196 bool __overflow = __builtin_mul_overflow(__r, __n2, &__r);
197 __glibcxx_assert(!__overflow);
198 return __r;
199 }
200
201#endif // C++17
202#endif // C++14
203
204#if __cplusplus > 201703L
205
206 // midpoint
207# define __cpp_lib_interpolate 201902L
208
209 template<typename _Tp>
210 constexpr
211 enable_if_t<__and_v<is_arithmetic<_Tp>, is_same<remove_cv_t<_Tp>, _Tp>,
212 __not_<is_same<_Tp, bool>>>,
213 _Tp>
214 midpoint(_Tp __a, _Tp __b) noexcept
215 {
216 if constexpr (is_integral_v<_Tp>)
217 {
218 using _Up = make_unsigned_t<_Tp>;
219
220 int __k = 1;
221 _Up __m = __a;
222 _Up __M = __b;
223 if (__a > __b)
224 {
225 __k = -1;
226 __m = __b;
227 __M = __a;
228 }
229 return __a + __k * _Tp(_Up(__M - __m) / 2);
230 }
231 else // is_floating
232 {
233 constexpr _Tp __lo = numeric_limits<_Tp>::min() * 2;
234 constexpr _Tp __hi = numeric_limits<_Tp>::max() / 2;
235 const _Tp __abs_a = __a < 0 ? -__a : __a;
236 const _Tp __abs_b = __b < 0 ? -__b : __b;
237 if (__abs_a <= __hi && __abs_b <= __hi) [[likely]]
238 return (__a + __b) / 2; // always correctly rounded
239 if (__abs_a < __lo) // not safe to halve __a
240 return __a + __b/2;
241 if (__abs_b < __lo) // not safe to halve __b
242 return __a/2 + __b;
243 return __a/2 + __b/2; // otherwise correctly rounded
244 }
245 }
246
247 template<typename _Tp>
248 constexpr enable_if_t<is_object_v<_Tp>, _Tp*>
249 midpoint(_Tp* __a, _Tp* __b) noexcept
250 {
251 static_assert( sizeof(_Tp) != 0, "type must be complete" );
252 return __a + (__b - __a) / 2;
253 }
254#endif // C++20
255
256#if __cplusplus >= 201703L
257
258#if __cplusplus > 201703L
259#define __cpp_lib_constexpr_numeric 201911L
260#endif
261
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
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
743// Feature test macro for parallel algorithms
744# define __cpp_lib_parallel_algorithm 201603L
745#endif // C++17
746
747#endif /* _GLIBCXX_NUMERIC */
748