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 
23 namespace boost { namespace algorithm {
24 
25 template <typename T>
identity_operation(std::multiplies<T>)26 T identity_operation ( std::multiplies<T> ) { return T(1); }
27 
28 template <typename T>
identity_operation(std::plus<T>)29 T 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
40 template <typename T, typename Integer>
41 typename boost::enable_if<boost::is_integral<Integer>, T>::type
power(T x,Integer n)42 power (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
67 template <typename T, typename Integer, typename Operation>
68 typename boost::enable_if<boost::is_integral<Integer>, T>::type
power(T x,Integer n,Operation op)69 power (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