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