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