1 // <numeric> -*- C++ -*-
3 // Copyright (C) 2001-2020 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>
63 #include <ext/numeric_traits.h>
65 #ifdef _GLIBCXX_PARALLEL
66 # include <parallel/numeric>
70 * @defgroup numerics Numerics
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.
77 #if __cplusplus >= 201402L
78 #include <type_traits>
80 namespace std _GLIBCXX_VISIBILITY(default)
82 _GLIBCXX_BEGIN_NAMESPACE_VERSION
86 // std::abs is not constexpr, doesn't support unsigned integers,
87 // and std::abs(std::numeric_limits<T>::min()) is undefined.
88 template<typename _Up, typename _Tp>
92 static_assert(is_unsigned<_Up>::value, "result type must be unsigned");
93 static_assert(sizeof(_Up) >= sizeof(_Tp),
94 "result type must be at least as wide as the input type");
95 return __val < 0 ? -(_Up)__val : (_Up)__val;
98 template<typename _Up> void __absu(bool) = delete;
100 // GCD implementation
101 template<typename _Tp>
103 __gcd(_Tp __m, _Tp __n)
105 static_assert(is_unsigned<_Tp>::value, "type must be unsigned");
106 return __m == 0 ? __n
108 : __detail::__gcd(__n, _Tp(__m % __n));
111 // LCM implementation
112 template<typename _Tp>
114 __lcm(_Tp __m, _Tp __n)
116 return (__m != 0 && __n != 0)
117 ? (__m / __detail::__gcd(__m, __n)) * __n
120 } // namespace __detail
122 #if __cplusplus >= 201703L
124 #define __cpp_lib_gcd_lcm 201606
125 // These were used in drafts of SD-6:
126 #define __cpp_lib_gcd 201606
127 #define __cpp_lib_lcm 201606
129 /// Greatest common divisor
130 template<typename _Mn, typename _Nn>
131 constexpr common_type_t<_Mn, _Nn>
132 gcd(_Mn __m, _Nn __n) noexcept
134 static_assert(is_integral_v<_Mn>, "std::gcd arguments must be integers");
135 static_assert(is_integral_v<_Nn>, "std::gcd arguments must be integers");
136 static_assert(_Mn(2) != _Mn(1), "std::gcd arguments must not be bool");
137 static_assert(_Nn(2) != _Nn(1), "std::gcd arguments must not be bool");
138 using _Up = make_unsigned_t<common_type_t<_Mn, _Nn>>;
139 return __detail::__gcd(__detail::__absu<_Up>(__m),
140 __detail::__absu<_Up>(__n));
143 /// Least common multiple
144 template<typename _Mn, typename _Nn>
145 constexpr common_type_t<_Mn, _Nn>
146 lcm(_Mn __m, _Nn __n) noexcept
148 static_assert(is_integral_v<_Mn>, "std::lcm arguments must be integers");
149 static_assert(is_integral_v<_Nn>, "std::lcm arguments must be integers");
150 static_assert(_Mn(2) == 2, "std::lcm arguments must not be bool");
151 static_assert(_Nn(2) == 2, "std::lcm arguments must not be bool");
152 using _Up = make_unsigned_t<common_type_t<_Mn, _Nn>>;
153 return __detail::__lcm(__detail::__absu<_Up>(__m),
154 __detail::__absu<_Up>(__n));
159 _GLIBCXX_END_NAMESPACE_VERSION
164 #if __cplusplus > 201703L
167 namespace std _GLIBCXX_VISIBILITY(default)
169 _GLIBCXX_BEGIN_NAMESPACE_VERSION
171 # define __cpp_lib_interpolate 201902L
173 template<typename _Tp>
175 enable_if_t<__and_v<is_arithmetic<_Tp>, is_same<remove_cv_t<_Tp>, _Tp>,
176 __not_<is_same<_Tp, bool>>>,
178 midpoint(_Tp __a, _Tp __b) noexcept
180 if constexpr (is_integral_v<_Tp>)
182 using _Up = make_unsigned_t<_Tp>;
193 return __a + __k * _Tp(_Up(__M - __m) / 2);
197 constexpr _Tp __lo = numeric_limits<_Tp>::min() * 2;
198 constexpr _Tp __hi = numeric_limits<_Tp>::max() / 2;
199 const _Tp __abs_a = __a < 0 ? -__a : __a;
200 const _Tp __abs_b = __b < 0 ? -__b : __b;
201 if (__abs_a <= __hi && __abs_b <= __hi) [[likely]]
202 return (__a + __b) / 2; // always correctly rounded
203 if (__abs_a < __lo) // not safe to halve __a
205 if (__abs_b < __lo) // not safe to halve __b
207 return __a/2 + __b/2; // otherwise correctly rounded
211 template<typename _Tp>
212 constexpr enable_if_t<is_object_v<_Tp>, _Tp*>
213 midpoint(_Tp* __a, _Tp* __b) noexcept
215 static_assert( sizeof(_Tp) != 0, "type must be complete" );
216 return __a + (__b - __a) / 2;
218 _GLIBCXX_END_NAMESPACE_VERSION
223 #if __cplusplus > 201402L
224 #include <bits/stl_function.h>
226 namespace std _GLIBCXX_VISIBILITY(default)
228 _GLIBCXX_BEGIN_NAMESPACE_VERSION
230 #if __cplusplus > 201703L
231 #define __cpp_lib_constexpr_numeric 201911L
234 /// @addtogroup numeric_ops
238 * @brief Calculate reduction of values in a range.
240 * @param __first Start of range.
241 * @param __last End of range.
242 * @param __init Starting value to add other values to.
243 * @param __binary_op A binary function object.
244 * @return The final sum.
246 * Reduce the values in the range `[first,last)` using a binary operation.
247 * The initial value is `init`. The values are not necessarily processed
250 * This algorithm is similar to `std::accumulate` but is not required to
251 * perform the operations in order from first to last. For operations
252 * that are commutative and associative the result will be the same as
253 * for `std::accumulate`, but for other operations (such as floating point
254 * arithmetic) the result can be different.
256 template<typename _InputIterator, typename _Tp, typename _BinaryOperation>
259 reduce(_InputIterator __first, _InputIterator __last, _Tp __init,
260 _BinaryOperation __binary_op)
262 using value_type = typename iterator_traits<_InputIterator>::value_type;
263 static_assert(is_invocable_r_v<_Tp, _BinaryOperation&, _Tp&, _Tp&>);
264 static_assert(is_convertible_v<value_type, _Tp>);
265 if constexpr (__is_random_access_iter<_InputIterator>::value)
267 while ((__last - __first) >= 4)
269 _Tp __v1 = __binary_op(__first[0], __first[1]);
270 _Tp __v2 = __binary_op(__first[2], __first[3]);
271 _Tp __v3 = __binary_op(__v1, __v2);
272 __init = __binary_op(__init, __v3);
276 for (; __first != __last; ++__first)
277 __init = __binary_op(__init, *__first);
282 * @brief Calculate reduction of values in a range.
284 * @param __first Start of range.
285 * @param __last End of range.
286 * @param __init Starting value to add other values to.
287 * @return The final sum.
289 * Reduce the values in the range `[first,last)` using addition.
290 * Equivalent to calling `std::reduce(first, last, init, std::plus<>())`.
292 template<typename _InputIterator, typename _Tp>
295 reduce(_InputIterator __first, _InputIterator __last, _Tp __init)
296 { return std::reduce(__first, __last, std::move(__init), plus<>()); }
299 * @brief Calculate reduction of values in a range.
301 * @param __first Start of range.
302 * @param __last End of range.
303 * @return The final sum.
305 * Reduce the values in the range `[first,last)` using addition, with
306 * an initial value of `T{}`, where `T` is the iterator's value type.
307 * Equivalent to calling `std::reduce(first, last, T{}, std::plus<>())`.
309 template<typename _InputIterator>
311 inline typename iterator_traits<_InputIterator>::value_type
312 reduce(_InputIterator __first, _InputIterator __last)
314 using value_type = typename iterator_traits<_InputIterator>::value_type;
315 return std::reduce(__first, __last, value_type{}, plus<>());
319 * @brief Combine elements from two ranges and reduce
321 * @param __first1 Start of first range.
322 * @param __last1 End of first range.
323 * @param __first2 Start of second range.
324 * @param __init Starting value to add other values to.
325 * @param __binary_op1 The function used to perform reduction.
326 * @param __binary_op2 The function used to combine values from the ranges.
327 * @return The final sum.
329 * Call `binary_op2(first1[n],first2[n])` for each `n` in `[0,last1-first1)`
330 * and then use `binary_op1` to reduce the values returned by `binary_op2`
331 * to a single value of type `T`.
333 * The range beginning at `first2` must contain at least `last1-first1`
336 template<typename _InputIterator1, typename _InputIterator2, typename _Tp,
337 typename _BinaryOperation1, typename _BinaryOperation2>
340 transform_reduce(_InputIterator1 __first1, _InputIterator1 __last1,
341 _InputIterator2 __first2, _Tp __init,
342 _BinaryOperation1 __binary_op1,
343 _BinaryOperation2 __binary_op2)
345 if constexpr (__and_v<__is_random_access_iter<_InputIterator1>,
346 __is_random_access_iter<_InputIterator2>>)
348 while ((__last1 - __first1) >= 4)
350 _Tp __v1 = __binary_op1(__binary_op2(__first1[0], __first2[0]),
351 __binary_op2(__first1[1], __first2[1]));
352 _Tp __v2 = __binary_op1(__binary_op2(__first1[2], __first2[2]),
353 __binary_op2(__first1[3], __first2[3]));
354 _Tp __v3 = __binary_op1(__v1, __v2);
355 __init = __binary_op1(__init, __v3);
360 for (; __first1 != __last1; ++__first1, (void) ++__first2)
361 __init = __binary_op1(__init, __binary_op2(*__first1, *__first2));
366 * @brief Combine elements from two ranges and reduce
368 * @param __first1 Start of first range.
369 * @param __last1 End of first range.
370 * @param __first2 Start of second range.
371 * @param __init Starting value to add other values to.
372 * @return The final sum.
374 * Call `first1[n]*first2[n]` for each `n` in `[0,last1-first1)` and then
375 * use addition to sum those products to a single value of type `T`.
377 * The range beginning at `first2` must contain at least `last1-first1`
380 template<typename _InputIterator1, typename _InputIterator2, typename _Tp>
383 transform_reduce(_InputIterator1 __first1, _InputIterator1 __last1,
384 _InputIterator2 __first2, _Tp __init)
386 return std::transform_reduce(__first1, __last1, __first2,
388 plus<>(), multiplies<>());
392 * @brief Transform the elements of a range and reduce
394 * @param __first Start of range.
395 * @param __last End of range.
396 * @param __init Starting value to add other values to.
397 * @param __binary_op The function used to perform reduction.
398 * @param __unary_op The function used to transform values from the range.
399 * @return The final sum.
401 * Call `unary_op(first[n])` for each `n` in `[0,last-first)` and then
402 * use `binary_op` to reduce the values returned by `unary_op`
403 * to a single value of type `T`.
405 template<typename _InputIterator, typename _Tp,
406 typename _BinaryOperation, typename _UnaryOperation>
409 transform_reduce(_InputIterator __first, _InputIterator __last, _Tp __init,
410 _BinaryOperation __binary_op, _UnaryOperation __unary_op)
412 if constexpr (__is_random_access_iter<_InputIterator>::value)
414 while ((__last - __first) >= 4)
416 _Tp __v1 = __binary_op(__unary_op(__first[0]),
417 __unary_op(__first[1]));
418 _Tp __v2 = __binary_op(__unary_op(__first[2]),
419 __unary_op(__first[3]));
420 _Tp __v3 = __binary_op(__v1, __v2);
421 __init = __binary_op(__init, __v3);
425 for (; __first != __last; ++__first)
426 __init = __binary_op(__init, __unary_op(*__first));
430 /** @brief Output the cumulative sum of one range to a second range
432 * @param __first Start of input range.
433 * @param __last End of input range.
434 * @param __result Start of output range.
435 * @param __init Initial value.
436 * @param __binary_op Function to perform summation.
437 * @return The end of the output range.
439 * Write the cumulative sum (aka prefix sum, aka scan) of the input range
440 * to the output range. Each element of the output range contains the
441 * running total of all earlier elements (and the initial value),
442 * using `binary_op` for summation.
444 * This function generates an "exclusive" scan, meaning the Nth element
445 * of the output range is the sum of the first N-1 input elements,
446 * so the Nth input element is not included.
448 template<typename _InputIterator, typename _OutputIterator, typename _Tp,
449 typename _BinaryOperation>
452 exclusive_scan(_InputIterator __first, _InputIterator __last,
453 _OutputIterator __result, _Tp __init,
454 _BinaryOperation __binary_op)
456 while (__first != __last)
459 __init = __binary_op(__init, *__first);
461 *__result++ = std::move(__v);
466 /** @brief Output the cumulative sum of one range to a second range
468 * @param __first Start of input range.
469 * @param __last End of input range.
470 * @param __result Start of output range.
471 * @param __init Initial value.
472 * @return The end of the output range.
474 * Write the cumulative sum (aka prefix sum, aka scan) of the input range
475 * to the output range. Each element of the output range contains the
476 * running total of all earlier elements (and the initial value),
477 * using `std::plus<>` for summation.
479 * This function generates an "exclusive" scan, meaning the Nth element
480 * of the output range is the sum of the first N-1 input elements,
481 * so the Nth input element is not included.
483 template<typename _InputIterator, typename _OutputIterator, typename _Tp>
485 inline _OutputIterator
486 exclusive_scan(_InputIterator __first, _InputIterator __last,
487 _OutputIterator __result, _Tp __init)
489 return std::exclusive_scan(__first, __last, __result, std::move(__init),
493 /** @brief Output the cumulative sum of one range to a second range
495 * @param __first Start of input range.
496 * @param __last End of input range.
497 * @param __result Start of output range.
498 * @param __binary_op Function to perform summation.
499 * @param __init Initial value.
500 * @return The end of the output range.
502 * Write the cumulative sum (aka prefix sum, aka scan) of the input range
503 * to the output range. Each element of the output range contains the
504 * running total of all earlier elements (and the initial value),
505 * using `binary_op` for summation.
507 * This function generates an "inclusive" scan, meaning the Nth element
508 * of the output range is the sum of the first N input elements,
509 * so the Nth input element is included.
511 template<typename _InputIterator, typename _OutputIterator,
512 typename _BinaryOperation, typename _Tp>
515 inclusive_scan(_InputIterator __first, _InputIterator __last,
516 _OutputIterator __result, _BinaryOperation __binary_op,
519 for (; __first != __last; ++__first)
520 *__result++ = __init = __binary_op(__init, *__first);
524 /** @brief Output the cumulative sum of one range to a second range
526 * @param __first Start of input range.
527 * @param __last End of input range.
528 * @param __result Start of output range.
529 * @param __binary_op Function to perform summation.
530 * @return The end of the output range.
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, using `binary_op` for summation.
536 * This function generates an "inclusive" scan, meaning the Nth element
537 * of the output range is the sum of the first N input elements,
538 * so the Nth input element is included.
540 template<typename _InputIterator, typename _OutputIterator,
541 typename _BinaryOperation>
544 inclusive_scan(_InputIterator __first, _InputIterator __last,
545 _OutputIterator __result, _BinaryOperation __binary_op)
547 if (__first != __last)
549 auto __init = *__first;
550 *__result++ = __init;
552 if (__first != __last)
553 __result = std::inclusive_scan(__first, __last, __result,
554 __binary_op, std::move(__init));
559 /** @brief Output the cumulative sum of one range to a second range
561 * @param __first Start of input range.
562 * @param __last End of input range.
563 * @param __result Start of output range.
564 * @return The end of the output range.
566 * Write the cumulative sum (aka prefix sum, aka scan) of the input range
567 * to the output range. Each element of the output range contains the
568 * running total of all earlier elements, using `std::plus<>` for summation.
570 * This function generates an "inclusive" scan, meaning the Nth element
571 * of the output range is the sum of the first N input elements,
572 * so the Nth input element is included.
574 template<typename _InputIterator, typename _OutputIterator>
576 inline _OutputIterator
577 inclusive_scan(_InputIterator __first, _InputIterator __last,
578 _OutputIterator __result)
579 { return std::inclusive_scan(__first, __last, __result, plus<>()); }
581 /** @brief Output the cumulative sum of one range to a second range
583 * @param __first Start of input range.
584 * @param __last End of input range.
585 * @param __result Start of output range.
586 * @param __init Initial value.
587 * @param __binary_op Function to perform summation.
588 * @param __unary_op Function to transform elements of the input range.
589 * @return The end of the output range.
591 * Write the cumulative sum (aka prefix sum, aka scan) of the input range
592 * to the output range. Each element of the output range contains the
593 * running total of all earlier elements (and the initial value),
594 * using `__unary_op` to transform the input elements
595 * and using `__binary_op` for summation.
597 * This function generates an "exclusive" scan, meaning the Nth element
598 * of the output range is the sum of the first N-1 input elements,
599 * so the Nth input element is not included.
601 template<typename _InputIterator, typename _OutputIterator, typename _Tp,
602 typename _BinaryOperation, typename _UnaryOperation>
605 transform_exclusive_scan(_InputIterator __first, _InputIterator __last,
606 _OutputIterator __result, _Tp __init,
607 _BinaryOperation __binary_op,
608 _UnaryOperation __unary_op)
610 while (__first != __last)
613 __init = __binary_op(__init, __unary_op(*__first));
615 *__result++ = std::move(__v);
620 /** @brief Output the cumulative sum of one range to a second range
622 * @param __first Start of input range.
623 * @param __last End of input range.
624 * @param __result Start of output range.
625 * @param __binary_op Function to perform summation.
626 * @param __unary_op Function to transform elements of the input range.
627 * @param __init Initial value.
628 * @return The end of the output range.
630 * Write the cumulative sum (aka prefix sum, aka scan) of the input range
631 * to the output range. Each element of the output range contains the
632 * running total of all earlier elements (and the initial value),
633 * using `__unary_op` to transform the input elements
634 * and using `__binary_op` for summation.
636 * This function generates an "inclusive" scan, meaning the Nth element
637 * of the output range is the sum of the first N input elements,
638 * so the Nth input element is included.
640 template<typename _InputIterator, typename _OutputIterator,
641 typename _BinaryOperation, typename _UnaryOperation, typename _Tp>
644 transform_inclusive_scan(_InputIterator __first, _InputIterator __last,
645 _OutputIterator __result,
646 _BinaryOperation __binary_op,
647 _UnaryOperation __unary_op,
650 for (; __first != __last; ++__first)
651 *__result++ = __init = __binary_op(__init, __unary_op(*__first));
655 /** @brief Output the cumulative sum of one range to a second range
657 * @param __first Start of input range.
658 * @param __last End of input range.
659 * @param __result Start of output range.
660 * @param __binary_op Function to perform summation.
661 * @param __unary_op Function to transform elements of the input range.
662 * @return The end of the output range.
664 * Write the cumulative sum (aka prefix sum, aka scan) of the input range
665 * to the output range. Each element of the output range contains the
666 * running total of all earlier elements,
667 * using `__unary_op` to transform the input elements
668 * and using `__binary_op` for summation.
670 * This function generates an "inclusive" scan, meaning the Nth element
671 * of the output range is the sum of the first N input elements,
672 * so the Nth input element is included.
674 template<typename _InputIterator, typename _OutputIterator,
675 typename _BinaryOperation, typename _UnaryOperation>
678 transform_inclusive_scan(_InputIterator __first, _InputIterator __last,
679 _OutputIterator __result,
680 _BinaryOperation __binary_op,
681 _UnaryOperation __unary_op)
683 if (__first != __last)
685 auto __init = __unary_op(*__first);
686 *__result++ = __init;
688 if (__first != __last)
689 __result = std::transform_inclusive_scan(__first, __last, __result,
690 __binary_op, __unary_op,
696 // @} group numeric_ops
698 _GLIBCXX_END_NAMESPACE_VERSION
701 // Parallel STL algorithms
702 # if _PSTL_EXECUTION_POLICIES_DEFINED
703 // If <execution> has already been included, pull in implementations
704 # include <pstl/glue_numeric_impl.h>
706 // Otherwise just pull in forward declarations
707 # include <pstl/glue_numeric_defs.h>
708 # define _PSTL_NUMERIC_FORWARD_DECLARED 1
711 // Feature test macro for parallel algorithms
712 # define __cpp_lib_parallel_algorithm 201603L
715 #endif /* _GLIBCXX_NUMERIC */