|
1 | 1 | /**
|
2 | 2 | * @file
|
3 |
| - * @brief Faster computation for \f$a^b\f$ |
4 |
| - * |
| 3 | + * @brief Exponentiating by squaring is a general method for fast computation of large positive integer powers of a number. |
| 4 | + * (https://en.wikipedia.org/wiki/Exponentiation_by_squaring) |
| 5 | + *@details |
5 | 6 | * Program that computes \f$a^b\f$ in \f$O(logN)\f$ time.
|
6 | 7 | * It is based on formula that:
|
7 | 8 | * 1. if \f$b\f$ is even:
|
|
10 | 11 | * \cdot a^\frac{b-1}{2} \cdot a = {a^\frac{b-1}{2}}^2 \cdot a\f$
|
11 | 12 | *
|
12 | 13 | * We can compute \f$a^b\f$ recursively using above algorithm.
|
| 14 | + * @author |
| 15 | + * @see |
13 | 16 | */
|
14 | 17 |
|
15 |
| -#include <cassert> |
16 |
| -#include <cmath> |
17 |
| -#include <cstdint> |
18 |
| -#include <cstdlib> |
19 |
| -#include <ctime> |
20 |
| -#include <iostream> |
| 18 | +#include <cassert> /// for assert |
| 19 | +#include <cmath> /// for std::pow |
| 20 | +#include <cstdint> /// for int64_t |
| 21 | +#include <cstdlib> /// for std::rand |
| 22 | +#include <ctime> /// for std::time |
| 23 | +#include <iostream> /// for IO operations |
| 24 | + |
21 | 25 |
|
22 | 26 | /**
|
23 |
| - * algorithm implementation for \f$a^b\f$ |
| 27 | + * @namespace math |
| 28 | + * @brief algorithm implementation for \f$a^b\f$ |
24 | 29 | */
|
25 |
| -template <typename T> |
26 |
| -double fast_power_recursive(T a, T b) { |
27 |
| - // negative power. a^b = 1 / (a^-b) |
28 |
| - if (b < 0) |
29 |
| - return 1.0 / fast_power_recursive(a, -b); |
30 |
| - |
31 |
| - if (b == 0) |
32 |
| - return 1; |
33 |
| - T bottom = fast_power_recursive(a, b >> 1); |
34 |
| - // Since it is integer division b/2 = (b-1)/2 where b is odd. |
35 |
| - // Therefore, case2 is easily solved by integer division. |
36 |
| - |
37 |
| - double result; |
38 |
| - if ((b & 1) == 0) // case1 |
39 |
| - result = bottom * bottom; |
40 |
| - else // case2 |
41 |
| - result = bottom * bottom * a; |
42 |
| - return result; |
43 |
| -} |
| 30 | + |
| 31 | +namespace math { |
| 32 | + |
| 33 | +/** |
| 34 | + * @brief Functions for fast computation of large positive integer powers of a number. |
| 35 | + * @param a The base |
| 36 | + * @param b The exponent |
| 37 | + * @returns The result of \f$a^b\f$ |
| 38 | + */ |
| 39 | + |
| 40 | + template <typename T> |
| 41 | + double fast_power_recursive(T a, T b) { |
| 42 | + /*When the base number is 0 and the exponent is non-positive, it is defined as meaningless |
| 43 | + */ |
| 44 | + if(a==0 && b<=0){ |
| 45 | + return NAN; |
| 46 | + } |
| 47 | + |
| 48 | + // negative power. a^b = 1 / (a^-b) |
| 49 | + if (b < 0) |
| 50 | + return 1.0 / fast_power_recursive(a, -b); |
| 51 | + |
| 52 | + if (b == 0) |
| 53 | + return 1; |
| 54 | + T bottom = fast_power_recursive(a, b >> 1); |
| 55 | + // Since it is integer division b/2 = (b-1)/2 where b is odd. |
| 56 | + // Therefore, case2 is easily solved by integer division. |
| 57 | + |
| 58 | + double result; |
| 59 | + if ((b & 1) == 0) // case1 |
| 60 | + result = bottom * bottom; |
| 61 | + else // case2 |
| 62 | + result = bottom * bottom * a; |
| 63 | + return result; |
| 64 | + } |
44 | 65 |
|
45 | 66 | /**
|
46 | 67 | Same algorithm with little different formula.
|
47 | 68 | It still calculates in \f$O(\log N)\f$
|
48 | 69 | */
|
49 |
| -template <typename T> |
50 |
| -double fast_power_linear(T a, T b) { |
51 |
| - // negative power. a^b = 1 / (a^-b) |
52 |
| - if (b < 0) |
53 |
| - return 1.0 / fast_power_linear(a, -b); |
54 |
| - |
55 |
| - double result = 1; |
56 |
| - while (b) { |
57 |
| - if (b & 1) |
58 |
| - result = result * a; |
59 |
| - a = a * a; |
60 |
| - b = b >> 1; |
| 70 | + template <typename T> |
| 71 | + double fast_power_linear(T a, T b) { |
| 72 | + /*When the base number is 0 and the exponent is non-positive, it is defined as meaningless |
| 73 | + */ |
| 74 | + if(a==0 && b<=0){ |
| 75 | + return NAN; |
| 76 | + } |
| 77 | + |
| 78 | + // negative power. a^b = 1 / (a^-b) |
| 79 | + if (b < 0) |
| 80 | + return 1.0 / fast_power_linear(a, -b); |
| 81 | + |
| 82 | + double result = 1; |
| 83 | + while (b) { |
| 84 | + if (b & 1) |
| 85 | + result = result * a; |
| 86 | + a = a * a; |
| 87 | + b = b >> 1; |
| 88 | + } |
| 89 | + return result; |
61 | 90 | }
|
62 |
| - return result; |
63 |
| -} |
| 91 | + |
| 92 | +}// namespace math |
64 | 93 |
|
65 | 94 | /**
|
66 |
| - * Main function |
| 95 | + * @brief Self-test implementations |
| 96 | + * @returns void |
67 | 97 | */
|
68 |
| -int main() { |
| 98 | +static void test() { |
| 99 | + /* The following program will generate and test 1000 pairs of random base and exponential combinations |
| 100 | + (ranging from -10 to 9), simulating power operations. The results of verifying fast_power_recursive(a, b) |
| 101 | + and fast_power_linear(a, b) are identical to those of the standard library functions std::pow(a, b) |
| 102 | + */ |
69 | 103 | std::srand(std::time(nullptr));
|
70 | 104 | std::ios_base::sync_with_stdio(false);
|
| 105 | + /*When the exponent is negative, it is often unreliable to use the == operator directly. |
| 106 | + When comparing comparison results, we use a small threshold (epsilon) to determine whether they are "close enough." |
| 107 | + */ |
| 108 | + const double epsilon = 1e-8; |
71 | 109 |
|
72 |
| - std::cout << "Testing..." << std::endl; |
73 |
| - for (int i = 0; i < 20; i++) { |
| 110 | + for (int i = 0; i < 1000; i++) { |
74 | 111 | int a = std::rand() % 20 - 10;
|
75 | 112 | int b = std::rand() % 20 - 10;
|
76 |
| - std::cout << std::endl << "Calculating " << a << "^" << b << std::endl; |
77 |
| - assert(fast_power_recursive(a, b) == std::pow(a, b)); |
78 |
| - assert(fast_power_linear(a, b) == std::pow(a, b)); |
| 113 | + /*When the base number is 0 and the exponent is non-positive, it is defined as meaningless |
| 114 | + */ |
| 115 | + if(a==0&&b<=0){ |
| 116 | + continue; |
| 117 | + } |
| 118 | + double result_recursive = math::fast_power_recursive(a, b); |
| 119 | + double result_linear = math::fast_power_linear(a, b); |
| 120 | + double result_pow = std::pow(a, b); |
79 | 121 |
|
80 |
| - std::cout << "------ " << a << "^" << b << " = " |
81 |
| - << fast_power_recursive(a, b) << std::endl; |
| 122 | + assert(std::fabs(result_recursive - result_pow) < epsilon); |
| 123 | + assert(std::fabs(result_linear - result_pow) < epsilon); |
82 | 124 | }
|
83 | 125 |
|
84 |
| - int64_t a, b; |
85 |
| - std::cin >> a >> b; |
86 |
| - |
87 |
| - std::cout << a << "^" << b << " = " << fast_power_recursive(a, b) |
88 |
| - << std::endl; |
| 126 | + std::cout << "All tests have successfully passed!\n"; |
| 127 | +} |
89 | 128 |
|
90 |
| - std::cout << a << "^" << b << " = " << fast_power_linear(a, b) << std::endl; |
| 129 | +/** |
| 130 | + * @brief Main function |
| 131 | + * @param argc commandline argument count (ignored) |
| 132 | + * @param argv commandline array of arguments (ignored) |
| 133 | + * @returns 0 on exit |
| 134 | + */ |
| 135 | +int main() { |
91 | 136 |
|
| 137 | + test(); // run self-test implementations |
| 138 | + // std::cout << math::fast_power_recursive(-10, -10) << "\n"<<std::pow(-10, -10)<<std::endl; |
92 | 139 | return 0;
|
93 | 140 | }
|
0 commit comments