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