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