1 // <numeric> -*- C++ -*-
3 // Copyright (C) 2001-2021 Free Software Foundation, Inc.
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)
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.
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.
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/>.
28 * Hewlett-Packard Company
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.
39 * Copyright (c) 1996,1997
40 * Silicon Graphics Computer Systems, Inc.
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.
51 /** @file include/numeric
52 * This is a Standard C++ Library header.
55 #ifndef _GLIBCXX_NUMERIC
56 #define _GLIBCXX_NUMERIC 1
58 #pragma GCC system_header
60 #include <bits/c++config.h>
61 #include <bits/stl_iterator_base_types.h>
62 #include <bits/stl_numeric.h>
64 #ifdef _GLIBCXX_PARALLEL
65 # include <parallel/numeric>
68 #if __cplusplus >= 201402L
69 # include <type_traits>
73 #if __cplusplus >= 201703L
74 # include <bits/stl_function.h>
77 #if __cplusplus > 201703L
82 * @defgroup numerics Numerics
84 * Components for performing numeric operations. Includes support for
85 * complex number types, random number generation, numeric (n-at-a-time)
86 * arrays, generalized numeric algorithms, and mathematical special functions.
89 namespace std _GLIBCXX_VISIBILITY(default)
91 _GLIBCXX_BEGIN_NAMESPACE_VERSION
93 #if __cplusplus >= 201402L
96 // std::abs is not constexpr, doesn't support unsigned integers,
97 // and std::abs(std::numeric_limits<T>::min()) is undefined.
98 template<typename _Up, typename _Tp>
102 static_assert(is_unsigned<_Up>::value, "result type must be unsigned");
103 static_assert(sizeof(_Up) >= sizeof(_Tp),
104 "result type must be at least as wide as the input type");
105 return __val < 0 ? -(_Up)__val : (_Up)__val;
108 template<typename _Up> void __absu(bool) = delete;
110 // GCD implementation, using Stein's algorithm
111 template<typename _Tp>
113 __gcd(_Tp __m, _Tp __n)
115 static_assert(is_unsigned<_Tp>::value, "type must be unsigned");
122 const int __i = std::__countr_zero(__m);
124 const int __j = std::__countr_zero(__n);
126 const int __k = __i < __j ? __i : __j; // min(i, j)
142 __n >>= std::__countr_zero(__n);
146 // LCM implementation
147 template<typename _Tp>
149 __lcm(_Tp __m, _Tp __n)
151 return (__m != 0 && __n != 0)
152 ? (__m / __detail::__gcd(__m, __n)) * __n
155 } // namespace __detail
157 #if __cplusplus >= 201703L
159 #define __cpp_lib_gcd_lcm 201606
160 // These were used in drafts of SD-6:
161 #define __cpp_lib_gcd 201606
162 #define __cpp_lib_lcm 201606
164 /// Greatest common divisor
165 template<typename _Mn, typename _Nn>
166 constexpr common_type_t<_Mn, _Nn>
167 gcd(_Mn __m, _Nn __n) noexcept
169 static_assert(is_integral_v<_Mn>, "std::gcd arguments must be integers");
170 static_assert(is_integral_v<_Nn>, "std::gcd arguments must be integers");
171 static_assert(_Mn(2) != _Mn(1), "std::gcd arguments must not be bool");
172 static_assert(_Nn(2) != _Nn(1), "std::gcd arguments must not be bool");
173 using _Up = make_unsigned_t<common_type_t<_Mn, _Nn>>;
174 return __detail::__gcd(__detail::__absu<_Up>(__m),
175 __detail::__absu<_Up>(__n));
178 /// Least common multiple
179 template<typename _Mn, typename _Nn>
180 constexpr common_type_t<_Mn, _Nn>
181 lcm(_Mn __m, _Nn __n) noexcept
183 static_assert(is_integral_v<_Mn>, "std::lcm arguments must be integers");
184 static_assert(is_integral_v<_Nn>, "std::lcm arguments must be integers");
185 static_assert(_Mn(2) == 2, "std::lcm arguments must not be bool");
186 static_assert(_Nn(2) == 2, "std::lcm arguments must not be bool");
187 using _Up = make_unsigned_t<common_type_t<_Mn, _Nn>>;
188 return __detail::__lcm(__detail::__absu<_Up>(__m),
189 __detail::__absu<_Up>(__n));
195 #if __cplusplus > 201703L
198 # define __cpp_lib_interpolate 201902L
200 template<typename _Tp>
202 enable_if_t<__and_v<is_arithmetic<_Tp>, is_same<remove_cv_t<_Tp>, _Tp>,
203 __not_<is_same<_Tp, bool>>>,
205 midpoint(_Tp __a, _Tp __b) noexcept
207 if constexpr (is_integral_v<_Tp>)
209 using _Up = make_unsigned_t<_Tp>;
220 return __a + __k * _Tp(_Up(__M - __m) / 2);
224 constexpr _Tp __lo = numeric_limits<_Tp>::min() * 2;
225 constexpr _Tp __hi = numeric_limits<_Tp>::max() / 2;
226 const _Tp __abs_a = __a < 0 ? -__a : __a;
227 const _Tp __abs_b = __b < 0 ? -__b : __b;
228 if (__abs_a <= __hi && __abs_b <= __hi) [[likely]]
229 return (__a + __b) / 2; // always correctly rounded
230 if (__abs_a < __lo) // not safe to halve __a
232 if (__abs_b < __lo) // not safe to halve __b
234 return __a/2 + __b/2; // otherwise correctly rounded
238 template<typename _Tp>
239 constexpr enable_if_t<is_object_v<_Tp>, _Tp*>
240 midpoint(_Tp* __a, _Tp* __b) noexcept
242 static_assert( sizeof(_Tp) != 0, "type must be complete" );
243 return __a + (__b - __a) / 2;
247 #if __cplusplus >= 201703L
249 #if __cplusplus > 201703L
250 #define __cpp_lib_constexpr_numeric 201911L
253 /// @addtogroup numeric_ops
257 * @brief Calculate reduction of values in a range.
259 * @param __first Start of range.
260 * @param __last End of range.
261 * @param __init Starting value to add other values to.
262 * @param __binary_op A binary function object.
263 * @return The final sum.
265 * Reduce the values in the range `[first,last)` using a binary operation.
266 * The initial value is `init`. The values are not necessarily processed
269 * This algorithm is similar to `std::accumulate` but is not required to
270 * perform the operations in order from first to last. For operations
271 * that are commutative and associative the result will be the same as
272 * for `std::accumulate`, but for other operations (such as floating point
273 * arithmetic) the result can be different.
275 template<typename _InputIterator, typename _Tp, typename _BinaryOperation>
278 reduce(_InputIterator __first, _InputIterator __last, _Tp __init,
279 _BinaryOperation __binary_op)
281 using value_type = typename iterator_traits<_InputIterator>::value_type;
282 static_assert(is_invocable_r_v<_Tp, _BinaryOperation&, _Tp&, _Tp&>);
283 static_assert(is_convertible_v<value_type, _Tp>);
284 if constexpr (__is_random_access_iter<_InputIterator>::value)
286 while ((__last - __first) >= 4)
288 _Tp __v1 = __binary_op(__first[0], __first[1]);
289 _Tp __v2 = __binary_op(__first[2], __first[3]);
290 _Tp __v3 = __binary_op(__v1, __v2);
291 __init = __binary_op(__init, __v3);
295 for (; __first != __last; ++__first)
296 __init = __binary_op(__init, *__first);
301 * @brief Calculate reduction of values in a range.
303 * @param __first Start of range.
304 * @param __last End of range.
305 * @param __init Starting value to add other values to.
306 * @return The final sum.
308 * Reduce the values in the range `[first,last)` using addition.
309 * Equivalent to calling `std::reduce(first, last, init, std::plus<>())`.
311 template<typename _InputIterator, typename _Tp>
314 reduce(_InputIterator __first, _InputIterator __last, _Tp __init)
315 { return std::reduce(__first, __last, std::move(__init), plus<>()); }
318 * @brief Calculate reduction of values in a range.
320 * @param __first Start of range.
321 * @param __last End of range.
322 * @return The final sum.
324 * Reduce the values in the range `[first,last)` using addition, with
325 * an initial value of `T{}`, where `T` is the iterator's value type.
326 * Equivalent to calling `std::reduce(first, last, T{}, std::plus<>())`.
328 template<typename _InputIterator>
330 inline typename iterator_traits<_InputIterator>::value_type
331 reduce(_InputIterator __first, _InputIterator __last)
333 using value_type = typename iterator_traits<_InputIterator>::value_type;
334 return std::reduce(__first, __last, value_type{}, plus<>());
338 * @brief Combine elements from two ranges and reduce
340 * @param __first1 Start of first range.
341 * @param __last1 End of first range.
342 * @param __first2 Start of second range.
343 * @param __init Starting value to add other values to.
344 * @param __binary_op1 The function used to perform reduction.
345 * @param __binary_op2 The function used to combine values from the ranges.
346 * @return The final sum.
348 * Call `binary_op2(first1[n],first2[n])` for each `n` in `[0,last1-first1)`
349 * and then use `binary_op1` to reduce the values returned by `binary_op2`
350 * to a single value of type `T`.
352 * The range beginning at `first2` must contain at least `last1-first1`
355 template<typename _InputIterator1, typename _InputIterator2, typename _Tp,
356 typename _BinaryOperation1, typename _BinaryOperation2>
359 transform_reduce(_InputIterator1 __first1, _InputIterator1 __last1,
360 _InputIterator2 __first2, _Tp __init,
361 _BinaryOperation1 __binary_op1,
362 _BinaryOperation2 __binary_op2)
364 if constexpr (__and_v<__is_random_access_iter<_InputIterator1>,
365 __is_random_access_iter<_InputIterator2>>)
367 while ((__last1 - __first1) >= 4)
369 _Tp __v1 = __binary_op1(__binary_op2(__first1[0], __first2[0]),
370 __binary_op2(__first1[1], __first2[1]));
371 _Tp __v2 = __binary_op1(__binary_op2(__first1[2], __first2[2]),
372 __binary_op2(__first1[3], __first2[3]));
373 _Tp __v3 = __binary_op1(__v1, __v2);
374 __init = __binary_op1(__init, __v3);
379 for (; __first1 != __last1; ++__first1, (void) ++__first2)
380 __init = __binary_op1(__init, __binary_op2(*__first1, *__first2));
385 * @brief Combine elements from two ranges and reduce
387 * @param __first1 Start of first range.
388 * @param __last1 End of first range.
389 * @param __first2 Start of second range.
390 * @param __init Starting value to add other values to.
391 * @return The final sum.
393 * Call `first1[n]*first2[n]` for each `n` in `[0,last1-first1)` and then
394 * use addition to sum those products to a single value of type `T`.
396 * The range beginning at `first2` must contain at least `last1-first1`
399 template<typename _InputIterator1, typename _InputIterator2, typename _Tp>
402 transform_reduce(_InputIterator1 __first1, _InputIterator1 __last1,
403 _InputIterator2 __first2, _Tp __init)
405 return std::transform_reduce(__first1, __last1, __first2,
407 plus<>(), multiplies<>());
411 * @brief Transform the elements of a range and reduce
413 * @param __first Start of range.
414 * @param __last End of range.
415 * @param __init Starting value to add other values to.
416 * @param __binary_op The function used to perform reduction.
417 * @param __unary_op The function used to transform values from the range.
418 * @return The final sum.
420 * Call `unary_op(first[n])` for each `n` in `[0,last-first)` and then
421 * use `binary_op` to reduce the values returned by `unary_op`
422 * to a single value of type `T`.
424 template<typename _InputIterator, typename _Tp,
425 typename _BinaryOperation, typename _UnaryOperation>
428 transform_reduce(_InputIterator __first, _InputIterator __last, _Tp __init,
429 _BinaryOperation __binary_op, _UnaryOperation __unary_op)
431 if constexpr (__is_random_access_iter<_InputIterator>::value)
433 while ((__last - __first) >= 4)
435 _Tp __v1 = __binary_op(__unary_op(__first[0]),
436 __unary_op(__first[1]));
437 _Tp __v2 = __binary_op(__unary_op(__first[2]),
438 __unary_op(__first[3]));
439 _Tp __v3 = __binary_op(__v1, __v2);
440 __init = __binary_op(__init, __v3);
444 for (; __first != __last; ++__first)
445 __init = __binary_op(__init, __unary_op(*__first));
449 /** @brief Output the cumulative sum of one range to a second range
451 * @param __first Start of input range.
452 * @param __last End of input range.
453 * @param __result Start of output range.
454 * @param __init Initial value.
455 * @param __binary_op Function to perform summation.
456 * @return The end of the output range.
458 * Write the cumulative sum (aka prefix sum, aka scan) of the input range
459 * to the output range. Each element of the output range contains the
460 * running total of all earlier elements (and the initial value),
461 * using `binary_op` for summation.
463 * This function generates an "exclusive" scan, meaning the Nth element
464 * of the output range is the sum of the first N-1 input elements,
465 * so the Nth input element is not included.
467 template<typename _InputIterator, typename _OutputIterator, typename _Tp,
468 typename _BinaryOperation>
471 exclusive_scan(_InputIterator __first, _InputIterator __last,
472 _OutputIterator __result, _Tp __init,
473 _BinaryOperation __binary_op)
475 while (__first != __last)
478 __init = __binary_op(__init, *__first);
480 *__result++ = std::move(__v);
485 /** @brief Output the cumulative sum of one range to a second range
487 * @param __first Start of input range.
488 * @param __last End of input range.
489 * @param __result Start of output range.
490 * @param __init Initial value.
491 * @return The end of the output range.
493 * Write the cumulative sum (aka prefix sum, aka scan) of the input range
494 * to the output range. Each element of the output range contains the
495 * running total of all earlier elements (and the initial value),
496 * using `std::plus<>` for summation.
498 * This function generates an "exclusive" scan, meaning the Nth element
499 * of the output range is the sum of the first N-1 input elements,
500 * so the Nth input element is not included.
502 template<typename _InputIterator, typename _OutputIterator, typename _Tp>
504 inline _OutputIterator
505 exclusive_scan(_InputIterator __first, _InputIterator __last,
506 _OutputIterator __result, _Tp __init)
508 return std::exclusive_scan(__first, __last, __result, std::move(__init),
512 /** @brief Output the cumulative sum of one range to a second range
514 * @param __first Start of input range.
515 * @param __last End of input range.
516 * @param __result Start of output range.
517 * @param __binary_op Function to perform summation.
518 * @param __init Initial value.
519 * @return The end of the output range.
521 * Write the cumulative sum (aka prefix sum, aka scan) of the input range
522 * to the output range. Each element of the output range contains the
523 * running total of all earlier elements (and the initial value),
524 * using `binary_op` for summation.
526 * This function generates an "inclusive" scan, meaning the Nth element
527 * of the output range is the sum of the first N input elements,
528 * so the Nth input element is included.
530 template<typename _InputIterator, typename _OutputIterator,
531 typename _BinaryOperation, typename _Tp>
534 inclusive_scan(_InputIterator __first, _InputIterator __last,
535 _OutputIterator __result, _BinaryOperation __binary_op,
538 for (; __first != __last; ++__first)
539 *__result++ = __init = __binary_op(__init, *__first);
543 /** @brief Output the cumulative sum of one range to a second range
545 * @param __first Start of input range.
546 * @param __last End of input range.
547 * @param __result Start of output range.
548 * @param __binary_op Function to perform summation.
549 * @return The end of the output range.
551 * Write the cumulative sum (aka prefix sum, aka scan) of the input range
552 * to the output range. Each element of the output range contains the
553 * running total of all earlier elements, using `binary_op` for summation.
555 * This function generates an "inclusive" scan, meaning the Nth element
556 * of the output range is the sum of the first N input elements,
557 * so the Nth input element is included.
559 template<typename _InputIterator, typename _OutputIterator,
560 typename _BinaryOperation>
563 inclusive_scan(_InputIterator __first, _InputIterator __last,
564 _OutputIterator __result, _BinaryOperation __binary_op)
566 if (__first != __last)
568 auto __init = *__first;
569 *__result++ = __init;
571 if (__first != __last)
572 __result = std::inclusive_scan(__first, __last, __result,
573 __binary_op, std::move(__init));
578 /** @brief Output the cumulative sum of one range to a second range
580 * @param __first Start of input range.
581 * @param __last End of input range.
582 * @param __result Start of output range.
583 * @return The end of the output range.
585 * Write the cumulative sum (aka prefix sum, aka scan) of the input range
586 * to the output range. Each element of the output range contains the
587 * running total of all earlier elements, using `std::plus<>` for summation.
589 * This function generates an "inclusive" scan, meaning the Nth element
590 * of the output range is the sum of the first N input elements,
591 * so the Nth input element is included.
593 template<typename _InputIterator, typename _OutputIterator>
595 inline _OutputIterator
596 inclusive_scan(_InputIterator __first, _InputIterator __last,
597 _OutputIterator __result)
598 { return std::inclusive_scan(__first, __last, __result, plus<>()); }
600 /** @brief Output the cumulative sum of one range to a second range
602 * @param __first Start of input range.
603 * @param __last End of input range.
604 * @param __result Start of output range.
605 * @param __init Initial value.
606 * @param __binary_op Function to perform summation.
607 * @param __unary_op Function to transform elements of the input range.
608 * @return The end of the output range.
610 * Write the cumulative sum (aka prefix sum, aka scan) of the input range
611 * to the output range. Each element of the output range contains the
612 * running total of all earlier elements (and the initial value),
613 * using `__unary_op` to transform the input elements
614 * and using `__binary_op` for summation.
616 * This function generates an "exclusive" scan, meaning the Nth element
617 * of the output range is the sum of the first N-1 input elements,
618 * so the Nth input element is not included.
620 template<typename _InputIterator, typename _OutputIterator, typename _Tp,
621 typename _BinaryOperation, typename _UnaryOperation>
624 transform_exclusive_scan(_InputIterator __first, _InputIterator __last,
625 _OutputIterator __result, _Tp __init,
626 _BinaryOperation __binary_op,
627 _UnaryOperation __unary_op)
629 while (__first != __last)
632 __init = __binary_op(__init, __unary_op(*__first));
634 *__result++ = std::move(__v);
639 /** @brief Output the cumulative sum of one range to a second range
641 * @param __first Start of input range.
642 * @param __last End of input range.
643 * @param __result Start of output range.
644 * @param __binary_op Function to perform summation.
645 * @param __unary_op Function to transform elements of the input range.
646 * @param __init Initial value.
647 * @return The end of the output range.
649 * Write the cumulative sum (aka prefix sum, aka scan) of the input range
650 * to the output range. Each element of the output range contains the
651 * running total of all earlier elements (and the initial value),
652 * using `__unary_op` to transform the input elements
653 * and using `__binary_op` for summation.
655 * This function generates an "inclusive" scan, meaning the Nth element
656 * of the output range is the sum of the first N input elements,
657 * so the Nth input element is included.
659 template<typename _InputIterator, typename _OutputIterator,
660 typename _BinaryOperation, typename _UnaryOperation, typename _Tp>
663 transform_inclusive_scan(_InputIterator __first, _InputIterator __last,
664 _OutputIterator __result,
665 _BinaryOperation __binary_op,
666 _UnaryOperation __unary_op,
669 for (; __first != __last; ++__first)
670 *__result++ = __init = __binary_op(__init, __unary_op(*__first));
674 /** @brief Output the cumulative sum of one range to a second range
676 * @param __first Start of input range.
677 * @param __last End of input range.
678 * @param __result Start of output range.
679 * @param __binary_op Function to perform summation.
680 * @param __unary_op Function to transform elements of the input range.
681 * @return The end of the output range.
683 * Write the cumulative sum (aka prefix sum, aka scan) of the input range
684 * to the output range. Each element of the output range contains the
685 * running total of all earlier elements,
686 * using `__unary_op` to transform the input elements
687 * and using `__binary_op` for summation.
689 * This function generates an "inclusive" scan, meaning the Nth element
690 * of the output range is the sum of the first N input elements,
691 * so the Nth input element is included.
693 template<typename _InputIterator, typename _OutputIterator,
694 typename _BinaryOperation, typename _UnaryOperation>
697 transform_inclusive_scan(_InputIterator __first, _InputIterator __last,
698 _OutputIterator __result,
699 _BinaryOperation __binary_op,
700 _UnaryOperation __unary_op)
702 if (__first != __last)
704 auto __init = __unary_op(*__first);
705 *__result++ = __init;
707 if (__first != __last)
708 __result = std::transform_inclusive_scan(__first, __last, __result,
709 __binary_op, __unary_op,
715 // @} group numeric_ops
718 _GLIBCXX_END_NAMESPACE_VERSION
721 #if __cplusplus >= 201703L
722 // Parallel STL algorithms
723 # if _PSTL_EXECUTION_POLICIES_DEFINED
724 // If <execution> has already been included, pull in implementations
725 # include <pstl/glue_numeric_impl.h>
727 // Otherwise just pull in forward declarations
728 # include <pstl/glue_numeric_defs.h>
729 # define _PSTL_NUMERIC_FORWARD_DECLARED 1
732 // Feature test macro for parallel algorithms
733 # define __cpp_lib_parallel_algorithm 201603L
736 #endif /* _GLIBCXX_NUMERIC */