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>
64 #ifdef _GLIBCXX_PARALLEL
65 # include <parallel/numeric>
69 * @defgroup numerics Numerics
71 * Components for performing numeric operations. Includes support for
72 * complex number types, random number generation, numeric (n-at-a-time)
73 * arrays, generalized numeric algorithms, and mathematical special functions.
76 #if __cplusplus >= 201402L
77 #include <type_traits>
79 namespace std _GLIBCXX_VISIBILITY(default)
81 _GLIBCXX_BEGIN_NAMESPACE_VERSION
85 // std::abs is not constexpr and doesn't support unsigned integers.
86 template<typename _Tp>
88 enable_if_t<__and_<is_integral<_Tp>, is_signed<_Tp>>::value, _Tp>
89 __abs_integral(_Tp __val)
90 { return __val < 0 ? -__val : __val; }
92 template<typename _Tp>
94 enable_if_t<__and_<is_integral<_Tp>, is_unsigned<_Tp>>::value, _Tp>
95 __abs_integral(_Tp __val)
98 void __abs_integral(bool) = delete;
100 template<typename _Mn, typename _Nn>
101 constexpr common_type_t<_Mn, _Nn>
102 __gcd(_Mn __m, _Nn __n)
104 return __m == 0 ? __detail::__abs_integral(__n)
105 : __n == 0 ? __detail::__abs_integral(__m)
106 : __detail::__gcd(__n, __m % __n);
109 /// Least common multiple
110 template<typename _Mn, typename _Nn>
111 constexpr common_type_t<_Mn, _Nn>
112 __lcm(_Mn __m, _Nn __n)
114 return (__m != 0 && __n != 0)
115 ? (__detail::__abs_integral(__m) / __detail::__gcd(__m, __n))
116 * __detail::__abs_integral(__n)
119 } // namespace __detail
121 #if __cplusplus >= 201703L
123 #define __cpp_lib_gcd_lcm 201606
124 // These were used in drafts of SD-6:
125 #define __cpp_lib_gcd 201606
126 #define __cpp_lib_lcm 201606
128 /// Greatest common divisor
129 template<typename _Mn, typename _Nn>
130 constexpr common_type_t<_Mn, _Nn>
131 gcd(_Mn __m, _Nn __n)
133 static_assert(is_integral_v<_Mn>, "gcd arguments are integers");
134 static_assert(is_integral_v<_Nn>, "gcd arguments are integers");
135 static_assert(!is_same_v<remove_cv_t<_Mn>, bool>,
136 "gcd arguments are not bools");
137 static_assert(!is_same_v<remove_cv_t<_Nn>, bool>,
138 "gcd arguments are not bools");
139 return __detail::__gcd(__m, __n);
142 /// Least common multiple
143 template<typename _Mn, typename _Nn>
144 constexpr common_type_t<_Mn, _Nn>
145 lcm(_Mn __m, _Nn __n)
147 static_assert(is_integral_v<_Mn>, "lcm arguments are integers");
148 static_assert(is_integral_v<_Nn>, "lcm arguments are integers");
149 static_assert(!is_same_v<remove_cv_t<_Mn>, bool>,
150 "lcm arguments are not bools");
151 static_assert(!is_same_v<remove_cv_t<_Nn>, bool>,
152 "lcm arguments are not bools");
153 return __detail::__lcm(__m, __n);
158 _GLIBCXX_END_NAMESPACE_VERSION
163 #if __cplusplus > 201703L
166 namespace std _GLIBCXX_VISIBILITY(default)
168 _GLIBCXX_BEGIN_NAMESPACE_VERSION
170 # define __cpp_lib_interpolate 201902L
172 template<typename _Tp>
174 enable_if_t<__and_v<is_arithmetic<_Tp>, is_same<remove_cv_t<_Tp>, _Tp>,
175 __not_<is_same<_Tp, bool>>>,
177 midpoint(_Tp __a, _Tp __b) noexcept
179 if constexpr (is_integral_v<_Tp>)
181 using _Up = make_unsigned_t<_Tp>;
192 return __a + __k * _Tp(_Up(__M - __m) / 2);
196 constexpr _Tp __lo = numeric_limits<_Tp>::min() * 2;
197 constexpr _Tp __hi = numeric_limits<_Tp>::max() / 2;
198 const _Tp __abs_a = __a < 0 ? -__a : __a;
199 const _Tp __abs_b = __b < 0 ? -__b : __b;
200 if (__abs_a <= __hi && __abs_b <= __hi) [[likely]]
201 return (__a + __b) / 2; // always correctly rounded
202 if (__abs_a < __lo) // not safe to halve __a
204 if (__abs_b < __lo) // not safe to halve __b
206 return __a/2 + __b/2; // otherwise correctly rounded
210 template<typename _Tp>
212 enable_if_t<__and_v<is_object<_Tp>, bool_constant<sizeof(_Tp) != 0>>, _Tp*>
213 midpoint(_Tp* __a, _Tp* __b) noexcept
215 return __a + (__b - __a) / 2;
217 _GLIBCXX_END_NAMESPACE_VERSION
222 #if __cplusplus > 201402L
223 #include <bits/stl_function.h>
225 namespace std _GLIBCXX_VISIBILITY(default)
227 _GLIBCXX_BEGIN_NAMESPACE_VERSION
229 /// @addtogroup numeric_ops
233 * @brief Calculate reduction of values in a range.
235 * @param __first Start of range.
236 * @param __last End of range.
237 * @param __init Starting value to add other values to.
238 * @param __binary_op A binary function object.
239 * @return The final sum.
241 * Reduce the values in the range `[first,last)` using a binary operation.
242 * The initial value is `init`. The values are not necessarily processed
245 * This algorithm is similar to `std::accumulate` but is not required to
246 * perform the operations in order from first to last. For operations
247 * that are commutative and associative the result will be the same as
248 * for `std::accumulate`, but for other operations (such as floating point
249 * arithmetic) the result can be different.
251 template<typename _InputIterator, typename _Tp, typename _BinaryOperation>
253 reduce(_InputIterator __first, _InputIterator __last, _Tp __init,
254 _BinaryOperation __binary_op)
256 using value_type = typename iterator_traits<_InputIterator>::value_type;
257 static_assert(is_invocable_r_v<_Tp, _BinaryOperation&, _Tp&, _Tp&>);
258 static_assert(is_convertible_v<value_type, _Tp>);
259 if constexpr (__is_random_access_iter<_InputIterator>::value)
261 while ((__last - __first) >= 4)
263 _Tp __v1 = __binary_op(__first[0], __first[1]);
264 _Tp __v2 = __binary_op(__first[2], __first[3]);
265 _Tp __v3 = __binary_op(__v1, __v2);
266 __init = __binary_op(__init, __v3);
270 for (; __first != __last; ++__first)
271 __init = __binary_op(__init, *__first);
276 * @brief Calculate reduction of values in a range.
278 * @param __first Start of range.
279 * @param __last End of range.
280 * @param __init Starting value to add other values to.
281 * @return The final sum.
283 * Reduce the values in the range `[first,last)` using addition.
284 * Equivalent to calling `std::reduce(first, last, init, std::plus<>())`.
286 template<typename _InputIterator, typename _Tp>
288 reduce(_InputIterator __first, _InputIterator __last, _Tp __init)
289 { return std::reduce(__first, __last, std::move(__init), plus<>()); }
292 * @brief Calculate reduction of values in a range.
294 * @param __first Start of range.
295 * @param __last End of range.
296 * @return The final sum.
298 * Reduce the values in the range `[first,last)` using addition, with
299 * an initial value of `T{}`, where `T` is the iterator's value type.
300 * Equivalent to calling `std::reduce(first, last, T{}, std::plus<>())`.
302 template<typename _InputIterator>
303 inline typename iterator_traits<_InputIterator>::value_type
304 reduce(_InputIterator __first, _InputIterator __last)
306 using value_type = typename iterator_traits<_InputIterator>::value_type;
307 return std::reduce(__first, __last, value_type{}, plus<>());
311 * @brief Combine elements from two ranges and reduce
313 * @param __first1 Start of first range.
314 * @param __last1 End of first range.
315 * @param __first2 Start of second range.
316 * @param __init Starting value to add other values to.
317 * @param __binary_op1 The function used to perform reduction.
318 * @param __binary_op2 The function used to combine values from the ranges.
319 * @return The final sum.
321 * Call `binary_op2(first1[n],first2[n])` for each `n` in `[0,last1-first1)`
322 * and then use `binary_op1` to reduce the values returned by `binary_op2`
323 * to a single value of type `T`.
325 * The range beginning at `first2` must contain at least `last1-first1`
328 template<typename _InputIterator1, typename _InputIterator2, typename _Tp,
329 typename _BinaryOperation1, typename _BinaryOperation2>
331 transform_reduce(_InputIterator1 __first1, _InputIterator1 __last1,
332 _InputIterator2 __first2, _Tp __init,
333 _BinaryOperation1 __binary_op1,
334 _BinaryOperation2 __binary_op2)
336 if constexpr (__and_v<__is_random_access_iter<_InputIterator1>,
337 __is_random_access_iter<_InputIterator2>>)
339 while ((__last1 - __first1) >= 4)
341 _Tp __v1 = __binary_op1(__binary_op2(__first1[0], __first2[0]),
342 __binary_op2(__first1[1], __first2[1]));
343 _Tp __v2 = __binary_op1(__binary_op2(__first1[2], __first2[2]),
344 __binary_op2(__first1[3], __first2[3]));
345 _Tp __v3 = __binary_op1(__v1, __v2);
346 __init = __binary_op1(__init, __v3);
351 for (; __first1 != __last1; ++__first1, (void) ++__first2)
352 __init = __binary_op1(__init, __binary_op2(*__first1, *__first2));
357 * @brief Combine elements from two ranges and reduce
359 * @param __first1 Start of first range.
360 * @param __last1 End of first range.
361 * @param __first2 Start of second range.
362 * @param __init Starting value to add other values to.
363 * @return The final sum.
365 * Call `first1[n]*first2[n]` for each `n` in `[0,last1-first1)` and then
366 * use addition to sum those products to a single value of type `T`.
368 * The range beginning at `first2` must contain at least `last1-first1`
371 template<typename _InputIterator1, typename _InputIterator2, typename _Tp>
373 transform_reduce(_InputIterator1 __first1, _InputIterator1 __last1,
374 _InputIterator2 __first2, _Tp __init)
376 return std::transform_reduce(__first1, __last1, __first2,
378 plus<>(), multiplies<>());
382 * @brief Transform the elements of a range and reduce
384 * @param __first Start of range.
385 * @param __last End of range.
386 * @param __init Starting value to add other values to.
387 * @param __binary_op The function used to perform reduction.
388 * @param __unary_op The function used to transform values from the range.
389 * @return The final sum.
391 * Call `unary_op(first[n])` for each `n` in `[0,last-first)` and then
392 * use `binary_op` to reduce the values returned by `unary_op`
393 * to a single value of type `T`.
395 template<typename _InputIterator, typename _Tp,
396 typename _BinaryOperation, typename _UnaryOperation>
398 transform_reduce(_InputIterator __first, _InputIterator __last, _Tp __init,
399 _BinaryOperation __binary_op, _UnaryOperation __unary_op)
401 if constexpr (__is_random_access_iter<_InputIterator>::value)
403 while ((__last - __first) >= 4)
405 _Tp __v1 = __binary_op(__unary_op(__first[0]),
406 __unary_op(__first[1]));
407 _Tp __v2 = __binary_op(__unary_op(__first[2]),
408 __unary_op(__first[3]));
409 _Tp __v3 = __binary_op(__v1, __v2);
410 __init = __binary_op(__init, __v3);
414 for (; __first != __last; ++__first)
415 __init = __binary_op(__init, __unary_op(*__first));
419 /** @brief Output the cumulative sum of one range to a second range
421 * @param __first Start of input range.
422 * @param __last End of input range.
423 * @param __result Start of output range.
424 * @param __init Initial value.
425 * @param __binary_op Function to perform summation.
426 * @return The end of the output range.
428 * Write the cumulative sum (aka prefix sum, aka scan) of the input range
429 * to the output range. Each element of the output range contains the
430 * running total of all earlier elements (and the initial value),
431 * using `binary_op` for summation.
433 * This function generates an "exclusive" scan, meaning the Nth element
434 * of the output range is the sum of the first N-1 input elements,
435 * so the Nth input element is not included.
437 template<typename _InputIterator, typename _OutputIterator, typename _Tp,
438 typename _BinaryOperation>
440 exclusive_scan(_InputIterator __first, _InputIterator __last,
441 _OutputIterator __result, _Tp __init,
442 _BinaryOperation __binary_op)
444 while (__first != __last)
447 __init = __binary_op(__init, *__first);
449 *__result++ = std::move(__v);
454 /** @brief Output the cumulative sum of one range to a second range
456 * @param __first Start of input range.
457 * @param __last End of input range.
458 * @param __result Start of output range.
459 * @param __init Initial value.
460 * @return The end of the output range.
462 * Write the cumulative sum (aka prefix sum, aka scan) of the input range
463 * to the output range. Each element of the output range contains the
464 * running total of all earlier elements (and the initial value),
465 * using `std::plus<>` for summation.
467 * This function generates an "exclusive" scan, meaning the Nth element
468 * of the output range is the sum of the first N-1 input elements,
469 * so the Nth input element is not included.
471 template<typename _InputIterator, typename _OutputIterator, typename _Tp>
472 inline _OutputIterator
473 exclusive_scan(_InputIterator __first, _InputIterator __last,
474 _OutputIterator __result, _Tp __init)
476 return std::exclusive_scan(__first, __last, __result, std::move(__init),
480 /** @brief Output the cumulative sum of one range to a second range
482 * @param __first Start of input range.
483 * @param __last End of input range.
484 * @param __result Start of output range.
485 * @param __binary_op Function to perform summation.
486 * @param __init Initial value.
487 * @return The end of the output range.
489 * Write the cumulative sum (aka prefix sum, aka scan) of the input range
490 * to the output range. Each element of the output range contains the
491 * running total of all earlier elements (and the initial value),
492 * using `binary_op` for summation.
494 * This function generates an "inclusive" scan, meaning the Nth element
495 * of the output range is the sum of the first N input elements,
496 * so the Nth input element is included.
498 template<typename _InputIterator, typename _OutputIterator,
499 typename _BinaryOperation, typename _Tp>
501 inclusive_scan(_InputIterator __first, _InputIterator __last,
502 _OutputIterator __result, _BinaryOperation __binary_op,
505 for (; __first != __last; ++__first)
506 *__result++ = __init = __binary_op(__init, *__first);
510 /** @brief Output the cumulative sum of one range to a second range
512 * @param __first Start of input range.
513 * @param __last End of input range.
514 * @param __result Start of output range.
515 * @param __binary_op Function to perform summation.
516 * @return The end of the output range.
518 * Write the cumulative sum (aka prefix sum, aka scan) of the input range
519 * to the output range. Each element of the output range contains the
520 * running total of all earlier elements, using `binary_op` for summation.
522 * This function generates an "inclusive" scan, meaning the Nth element
523 * of the output range is the sum of the first N input elements,
524 * so the Nth input element is included.
526 template<typename _InputIterator, typename _OutputIterator,
527 typename _BinaryOperation>
529 inclusive_scan(_InputIterator __first, _InputIterator __last,
530 _OutputIterator __result, _BinaryOperation __binary_op)
532 if (__first != __last)
534 auto __init = *__first;
535 *__result++ = __init;
537 if (__first != __last)
538 __result = std::inclusive_scan(__first, __last, __result,
539 __binary_op, std::move(__init));
544 /** @brief Output the cumulative sum of one range to a second range
546 * @param __first Start of input range.
547 * @param __last End of input range.
548 * @param __result Start of output range.
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 `std::plus<>` 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 inline _OutputIterator
561 inclusive_scan(_InputIterator __first, _InputIterator __last,
562 _OutputIterator __result)
563 { return std::inclusive_scan(__first, __last, __result, plus<>()); }
565 /** @brief Output the cumulative sum of one range to a second range
567 * @param __first Start of input range.
568 * @param __last End of input range.
569 * @param __result Start of output range.
570 * @param __init Initial value.
571 * @param __binary_op Function to perform summation.
572 * @param __unary_op Function to transform elements of the input range.
573 * @return The end of the output range.
575 * Write the cumulative sum (aka prefix sum, aka scan) of the input range
576 * to the output range. Each element of the output range contains the
577 * running total of all earlier elements (and the initial value),
578 * using `__unary_op` to transform the input elements
579 * and using `__binary_op` for summation.
581 * This function generates an "exclusive" scan, meaning the Nth element
582 * of the output range is the sum of the first N-1 input elements,
583 * so the Nth input element is not included.
585 template<typename _InputIterator, typename _OutputIterator, typename _Tp,
586 typename _BinaryOperation, typename _UnaryOperation>
588 transform_exclusive_scan(_InputIterator __first, _InputIterator __last,
589 _OutputIterator __result, _Tp __init,
590 _BinaryOperation __binary_op,
591 _UnaryOperation __unary_op)
593 while (__first != __last)
596 __init = __binary_op(__init, __unary_op(*__first));
598 *__result++ = std::move(__v);
603 /** @brief Output the cumulative sum of one range to a second range
605 * @param __first Start of input range.
606 * @param __last End of input range.
607 * @param __result Start of output range.
608 * @param __binary_op Function to perform summation.
609 * @param __unary_op Function to transform elements of the input range.
610 * @param __init Initial value.
611 * @return The end of the output range.
613 * Write the cumulative sum (aka prefix sum, aka scan) of the input range
614 * to the output range. Each element of the output range contains the
615 * running total of all earlier elements (and the initial value),
616 * using `__unary_op` to transform the input elements
617 * and using `__binary_op` for summation.
619 * This function generates an "inclusive" scan, meaning the Nth element
620 * of the output range is the sum of the first N input elements,
621 * so the Nth input element is included.
623 template<typename _InputIterator, typename _OutputIterator,
624 typename _BinaryOperation, typename _UnaryOperation, typename _Tp>
626 transform_inclusive_scan(_InputIterator __first, _InputIterator __last,
627 _OutputIterator __result,
628 _BinaryOperation __binary_op,
629 _UnaryOperation __unary_op,
632 for (; __first != __last; ++__first)
633 *__result++ = __init = __binary_op(__init, __unary_op(*__first));
637 /** @brief Output the cumulative sum of one range to a second range
639 * @param __first Start of input range.
640 * @param __last End of input range.
641 * @param __result Start of output range.
642 * @param __binary_op Function to perform summation.
643 * @param __unary_op Function to transform elements of the input range.
644 * @return The end of the output range.
646 * Write the cumulative sum (aka prefix sum, aka scan) of the input range
647 * to the output range. Each element of the output range contains the
648 * running total of all earlier elements,
649 * using `__unary_op` to transform the input elements
650 * and using `__binary_op` for summation.
652 * This function generates an "inclusive" scan, meaning the Nth element
653 * of the output range is the sum of the first N input elements,
654 * so the Nth input element is included.
656 template<typename _InputIterator, typename _OutputIterator,
657 typename _BinaryOperation, typename _UnaryOperation>
659 transform_inclusive_scan(_InputIterator __first, _InputIterator __last,
660 _OutputIterator __result,
661 _BinaryOperation __binary_op,
662 _UnaryOperation __unary_op)
664 if (__first != __last)
666 auto __init = __unary_op(*__first);
667 *__result++ = __init;
669 if (__first != __last)
670 __result = std::transform_inclusive_scan(__first, __last, __result,
671 __binary_op, __unary_op,
677 // @} group numeric_ops
679 _GLIBCXX_END_NAMESPACE_VERSION
682 // Parallel STL algorithms
683 # if _PSTL_EXECUTION_POLICIES_DEFINED
684 // If <execution> has already been included, pull in implementations
685 # include <pstl/glue_numeric_impl.h>
687 // Otherwise just pull in forward declarations
688 # include <pstl/glue_numeric_defs.h>
689 # define _PSTL_NUMERIC_FORWARD_DECLARED 1
692 // Feature test macro for parallel algorithms
693 # define __cpp_lib_parallel_algorithm 201603L
696 #endif /* _GLIBCXX_NUMERIC */