1/*
2 Copyright (c) Marshall Clow 2014.
3
4 Distributed under the Boost Software License, Version 1.0. (See accompanying
5 file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
6
7 Revision history:
8 2 Dec 2014 mtc First version; power
9
10*/
11
12/// \file algorithm.hpp
13/// \brief Misc Algorithms
14/// \author Marshall Clow
15///
16
17#ifndef BOOST_ALGORITHM_HPP
18#define BOOST_ALGORITHM_HPP
19
20#include <boost/utility/enable_if.hpp> // for boost::disable_if
21#include <boost/type_traits/is_integral.hpp>
22
23namespace boost { namespace algorithm {
24
25template <typename T>
26T identity_operation ( std::multiplies<T> ) { return T(1); }
27
28template <typename T>
29T identity_operation ( std::plus<T> ) { return T(0); }
30
31
32/// \fn power ( T x, Integer n )
33/// \return the value "x" raised to the power "n"
34///
35/// \param x The value to be exponentiated
36/// \param n The exponent (must be >= 0)
37///
38// \remark Taken from Knuth, The Art of Computer Programming, Volume 2:
39// Seminumerical Algorithms, Section 4.6.3
40template <typename T, typename Integer>
41typename boost::enable_if<boost::is_integral<Integer>, T>::type
42power (T x, Integer n) {
43 T y = 1; // Should be "T y{1};"
44 if (n == 0) return y;
45 while (true) {
46 if (n % 2 == 1) {
47 y = x * y;
48 if (n == 1)
49 return y;
50 }
51 n = n / 2;
52 x = x * x;
53 }
54 return y;
55 }
56
57/// \fn power ( T x, Integer n, Operation op )
58/// \return the value "x" raised to the power "n"
59/// using the operaton "op".
60///
61/// \param x The value to be exponentiated
62/// \param n The exponent (must be >= 0)
63/// \param op The operation used
64///
65// \remark Taken from Knuth, The Art of Computer Programming, Volume 2:
66// Seminumerical Algorithms, Section 4.6.3
67template <typename T, typename Integer, typename Operation>
68typename boost::enable_if<boost::is_integral<Integer>, T>::type
69power (T x, Integer n, Operation op) {
70 T y = identity_operation(op);
71 if (n == 0) return y;
72 while (true) {
73 if (n % 2 == 1) {
74 y = op(x, y);
75 if (n == 1)
76 return y;
77 }
78 n = n / 2;
79 x = op(x, x);
80 }
81 return y;
82 }
83
84}}
85
86#endif // BOOST_ALGORITHM_HPP
87