libstdc++-v3/include/std/numeric

   1 // <numeric> -*- C++ -*-
   2
   3 // Copyright (C) 2001-2024 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 #define __glibcxx_want_constexpr_numeric
  83 #define __glibcxx_want_gcd
  84 #define __glibcxx_want_gcd_lcm
  85 #define __glibcxx_want_interpolate
  86 #define __glibcxx_want_lcm
  87 #define __glibcxx_want_parallel_algorithm
  88 #define __glibcxx_want_ranges_iota
  89 #define __glibcxx_want_saturation_arithmetic
  90 #include <bits/version.h>
  91
  92 #ifdef __glibcxx_saturation_arithmetic // C++ >= 26
  93 # include <bits/sat_arith.h>
  94 #endif
  95
  96 /**
  97  * @defgroup numerics Numerics
  98  *
  99  * Components for performing numeric operations. Includes support for
 100  * complex number types, random number generation, numeric (n-at-a-time)
 101  * arrays, generalized numeric algorithms, and mathematical special functions.
 102  */
 103
 104 namespace std _GLIBCXX_VISIBILITY(default)
 105 {
 106 _GLIBCXX_BEGIN_NAMESPACE_VERSION
 107
 108 #if __cplusplus >= 201402L
 109 namespace __detail
 110 {
 111   // Like std::abs, but supports unsigned types and returns the specified type,
 112   // so |std::numeric_limits<_Tp>::min()| is OK if representable in _Res.
 113   template<typename _Res, typename _Tp>
 114     constexpr _Res
 115     __abs_r(_Tp __val)
 116     {
 117       static_assert(sizeof(_Res) >= sizeof(_Tp),
 118           "result type must be at least as wide as the input type");
 119
 120       if (__val >= 0)
 121         return __val;
 122 #ifdef _GLIBCXX_ASSERTIONS
 123       if (!__is_constant_evaluated()) // overflow already detected in constexpr
 124         __glibcxx_assert(__val != __gnu_cxx::__int_traits<_Res>::__min);
 125 #endif
 126       return -static_cast<_Res>(__val);
 127     }
 128
 129   template<typename> void __abs_r(bool) = delete;
 130
 131   // GCD implementation, using Stein's algorithm
 132   template<typename _Tp>
 133     constexpr _Tp
 134     __gcd(_Tp __m, _Tp __n)
 135     {
 136       static_assert(is_unsigned<_Tp>::value, "type must be unsigned");
 137
 138       if (__m == 0)
 139         return __n;
 140       if (__n == 0)
 141         return __m;
 142
 143       const int __i = std::__countr_zero(__m);
 144       __m >>= __i;
 145       const int __j = std::__countr_zero(__n);
 146       __n >>= __j;
 147       const int __k = __i < __j ? __i : __j; // min(i, j)
 148
 149       while (true)
 150         {
 151           if (__m > __n)
 152             {
 153               _Tp __tmp = __m;
 154               __m = __n;
 155               __n = __tmp;
 156             }
 157
 158           __n -= __m;
 159
 160           if (__n == 0)
 161             return __m << __k;
 162
 163           __n >>= std::__countr_zero(__n);
 164         }
 165     }
 166 } // namespace __detail
 167 #endif // C++14
 168
 169 #ifdef __cpp_lib_gcd_lcm // C++ >= 17
 170   /// Greatest common divisor
 171   template<typename _Mn, typename _Nn>
 172     constexpr common_type_t<_Mn, _Nn>
 173     gcd(_Mn __m, _Nn __n) noexcept
 174     {
 175       static_assert(is_integral_v<_Mn> && is_integral_v<_Nn>,
 176                     "std::gcd arguments must be integers");
 177       static_assert(_Mn(2) == 2 && _Nn(2) == 2,
 178                     "std::gcd arguments must not be bool");
 179       using _Ct = common_type_t<_Mn, _Nn>;
 180       const _Ct __m2 = __detail::__abs_r<_Ct>(__m);
 181       const _Ct __n2 = __detail::__abs_r<_Ct>(__n);
 182       return __detail::__gcd<make_unsigned_t<_Ct>>(__m2, __n2);
 183     }
 184
 185   /// Least common multiple
 186   template<typename _Mn, typename _Nn>
 187     constexpr common_type_t<_Mn, _Nn>
 188     lcm(_Mn __m, _Nn __n) noexcept
 189     {
 190       static_assert(is_integral_v<_Mn> && is_integral_v<_Nn>,
 191                     "std::lcm arguments must be integers");
 192       static_assert(_Mn(2) == 2 && _Nn(2) == 2,
 193                     "std::lcm arguments must not be bool");
 194       using _Ct = common_type_t<_Mn, _Nn>;
 195       const _Ct __m2 = __detail::__abs_r<_Ct>(__m);
 196       const _Ct __n2 = __detail::__abs_r<_Ct>(__n);
 197       if (__m2 == 0 || __n2 == 0)
 198         return 0;
 199       _Ct __r = __m2 / __detail::__gcd<make_unsigned_t<_Ct>>(__m2, __n2);
 200
 201       if constexpr (is_signed_v<_Ct>)
 202         if (__is_constant_evaluated())
 203           return __r * __n2; // constant evaluation can detect overflow here.
 204
 205       bool __overflow = __builtin_mul_overflow(__r, __n2, &__r);
 206       __glibcxx_assert(!__overflow);
 207       return __r;
 208     }
 209
 210 #endif // __cpp_lib_gcd_lcm
 211
 212   // midpoint
 213 #ifdef __cpp_lib_interpolate // C++ >= 20
 214   template<typename _Tp>
 215     constexpr
 216     enable_if_t<__and_v<is_arithmetic<_Tp>, is_same<remove_cv_t<_Tp>, _Tp>,
 217                         __not_<is_same<_Tp, bool>>>,
 218                 _Tp>
 219     midpoint(_Tp __a, _Tp __b) noexcept
 220     {
 221       if constexpr (is_integral_v<_Tp>)
 222         {
 223           using _Up = make_unsigned_t<_Tp>;
 224
 225           int __k = 1;
 226           _Up __m = __a;
 227           _Up __M = __b;
 228           if (__a > __b)
 229             {
 230               __k = -1;
 231               __m = __b;
 232               __M = __a;
 233             }
 234           return __a + __k * _Tp(_Up(__M - __m) / 2);
 235         }
 236       else // is_floating
 237         {
 238           constexpr _Tp __lo = numeric_limits<_Tp>::min() * 2;
 239           constexpr _Tp __hi = numeric_limits<_Tp>::max() / 2;
 240           const _Tp __abs_a = __a < 0 ? -__a : __a;
 241           const _Tp __abs_b = __b < 0 ? -__b : __b;
 242           if (__abs_a <= __hi && __abs_b <= __hi) [[likely]]
 243             return (__a + __b) / 2; // always correctly rounded
 244           if (__abs_a < __lo) // not safe to halve __a
 245             return __a + __b/2;
 246           if (__abs_b < __lo) // not safe to halve __b
 247             return __a/2 + __b;
 248           return __a/2 + __b/2;     // otherwise correctly rounded
 249         }
 250     }
 251
 252   template<typename _Tp>
 253     constexpr enable_if_t<is_object_v<_Tp>, _Tp*>
 254     midpoint(_Tp* __a, _Tp* __b) noexcept
 255     {
 256       static_assert( sizeof(_Tp) != 0, "type must be complete" );
 257       return __a  + (__b - __a) / 2;
 258     }
 259 #endif // __cpp_lib_interpolate
 260
 261 #if __cplusplus >= 201703L
 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 && _GLIBCXX_HOSTED
 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 #endif // C++17
 743
 744 #endif /* _GLIBCXX_NUMERIC */