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
104namespace std _GLIBCXX_VISIBILITY(default)
105{
106_GLIBCXX_BEGIN_NAMESPACE_VERSION
107
108#if __cplusplus >= 201402L
109namespace __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 */
745