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