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