1 /* Copyright 2015 The TensorFlow Authors. All Rights Reserved.
2
3 Licensed under the Apache License, Version 2.0 (the "License");
4 you may not use this file except in compliance with the License.
5 You may obtain a copy of the License at
6
7 http://www.apache.org/licenses/LICENSE-2.0
8
9 Unless required by applicable law or agreed to in writing, software
10 distributed under the License is distributed on an "AS IS" BASIS,
11 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12 See the License for the specific language governing permissions and
13 limitations under the License.
14 ==============================================================================*/
15
16 #include "tensorflow/core/lib/math/math_util.h"
17
18 #include <cmath>
19 #include <limits>
20 #include <vector>
21
22 #include "tensorflow/core/platform/logging.h"
23 #include "tensorflow/core/platform/test.h"
24 #include "tensorflow/core/platform/test_benchmark.h"
25 #include "tensorflow/core/platform/types.h"
26
27 namespace tensorflow {
28 namespace {
29
30 // Number of arguments for each test of the CeilOrRatio method
31 const int kNumTestArguments = 4;
32
33 template <typename IntegralType, typename TestDataType>
TestCeilOfRatio(const TestDataType test_data[][kNumTestArguments],int num_tests)34 void TestCeilOfRatio(const TestDataType test_data[][kNumTestArguments],
35 int num_tests) {
36 for (int i = 0; i < num_tests; ++i) {
37 const IntegralType numerator = test_data[i][0];
38 const IntegralType denominator = test_data[i][1];
39 const IntegralType expected_floor = test_data[i][2];
40 const IntegralType expected_ceil = test_data[i][3];
41 // Make sure the two ways to compute the floor return the same thing.
42 IntegralType floor_1 = MathUtil::FloorOfRatio(numerator, denominator);
43 IntegralType floor_2 = MathUtil::CeilOrFloorOfRatio<IntegralType, false>(
44 numerator, denominator);
45 EXPECT_EQ(floor_1, floor_2);
46 EXPECT_EQ(expected_floor, floor_1)
47 << "FloorOfRatio fails with numerator = " << numerator
48 << ", denominator = " << denominator << " "
49 << (8 * sizeof(IntegralType)) << " bits";
50 IntegralType ceil_1 = MathUtil::CeilOfRatio(numerator, denominator);
51 IntegralType ceil_2 = MathUtil::CeilOrFloorOfRatio<IntegralType, true>(
52 numerator, denominator);
53 EXPECT_EQ(ceil_1, ceil_2);
54 EXPECT_EQ(expected_ceil, ceil_1)
55 << "CeilOfRatio fails with numerator = " << numerator
56 << ", denominator = " << denominator << " "
57 << (8 * sizeof(IntegralType)) << " bits";
58 }
59 }
60
61 template <typename UnsignedIntegralType>
TestCeilOfRatioUnsigned(uint64 kMax)62 void TestCeilOfRatioUnsigned(uint64 kMax) {
63 const int kNumTests = 12;
64 const uint64 kTestData[kNumTests][kNumTestArguments] = {
65 // Numerator | Denominator | Expected floor of ratio | Expected ceil of
66 // ratio |
67 // When numerator = 0, the result is always zero
68 {0, 1, 0, 0},
69 {0, 2, 0, 0},
70 {0, kMax, 0, 0},
71 // Try some non-extreme cases
72 {1, 1, 1, 1},
73 {5, 2, 2, 3},
74 // Try with huge positive numerator
75 {kMax, 1, kMax, kMax},
76 {kMax, 2, kMax / 2, kMax / 2 + ((kMax % 2 != 0) ? 1 : 0)},
77 {kMax, 3, kMax / 3, kMax / 3 + ((kMax % 3 != 0) ? 1 : 0)},
78 // Try with a huge positive denominator
79 {1, kMax, 0, 1},
80 {2, kMax, 0, 1},
81 {3, kMax, 0, 1},
82 // Try with a huge numerator and a huge denominator
83 {kMax, kMax, 1, 1},
84 };
85 TestCeilOfRatio<UnsignedIntegralType, uint64>(kTestData, kNumTests);
86 }
87
88 template <typename SignedInteger>
TestCeilOfRatioSigned(int64 kMin,int64 kMax)89 void TestCeilOfRatioSigned(int64 kMin, int64 kMax) {
90 const int kNumTests = 30;
91 const int64 kTestData[kNumTests][kNumTestArguments] = {
92 // Numerator | Denominator | Expected floor of ratio | Expected ceil of
93 // ratio |
94 // When numerator = 0, the result is always zero
95 {0, 1, 0, 0},
96 {0, -1, 0, 0},
97 {0, 2, 0, 0},
98 {0, kMin, 0, 0},
99 {0, kMax, 0, 0},
100 // Try all four combinations of 1 and -1
101 {1, 1, 1, 1},
102 {-1, 1, -1, -1},
103 {1, -1, -1, -1},
104 {-1, -1, 1, 1},
105 // Try all four combinations of +/-5 divided by +/- 2
106 {5, 2, 2, 3},
107 {-5, 2, -3, -2},
108 {5, -2, -3, -2},
109 {-5, -2, 2, 3},
110 // Try with huge positive numerator
111 {kMax, 1, kMax, kMax},
112 {kMax, -1, -kMax, -kMax},
113 {kMax, 2, kMax / 2, kMax / 2 + ((kMax % 2 != 0) ? 1 : 0)},
114 {kMax, 3, kMax / 3, kMax / 3 + ((kMax % 3 != 0) ? 1 : 0)},
115 // Try with huge negative numerator
116 {kMin, 1, kMin, kMin},
117 {kMin, 2, kMin / 2 - ((kMin % 2 != 0) ? 1 : 0), kMin / 2},
118 {kMin, 3, kMin / 3 - ((kMin % 3 != 0) ? 1 : 0), kMin / 3},
119 // Try with a huge positive denominator
120 {1, kMax, 0, 1},
121 {2, kMax, 0, 1},
122 {3, kMax, 0, 1},
123 // Try with a huge negative denominator
124 {1, kMin, -1, 0},
125 {2, kMin, -1, 0},
126 {3, kMin, -1, 0},
127 // Try with a huge numerator and a huge denominator
128 {kMin, kMin, 1, 1},
129 {kMin, kMax, -2, -1},
130 {kMax, kMin, -1, 0},
131 {kMax, kMax, 1, 1},
132 };
133 TestCeilOfRatio<SignedInteger, int64>(kTestData, kNumTests);
134 }
135
136 // ------------------------------------------------------------------------ //
137 // Benchmarking CeilOrFloorOfRatio
138 //
139 // We compare with other implementations that are unsafe in general.
140 // ------------------------------------------------------------------------ //
141
142 // An implementation of CeilOfRatio that is correct for small enough values,
143 // and provided that the numerator and denominator are both positive
144 template <typename IntegralType>
CeilOfRatioDenomMinusOne(IntegralType numerator,IntegralType denominator)145 static IntegralType CeilOfRatioDenomMinusOne(IntegralType numerator,
146 IntegralType denominator) {
147 const IntegralType kOne(1);
148 return (numerator + denominator - kOne) / denominator;
149 }
150
151 // An implementation of FloorOfRatio that is correct when the denominator is
152 // positive and the numerator non-negative
153 template <typename IntegralType>
FloorOfRatioByDivision(IntegralType numerator,IntegralType denominator)154 static IntegralType FloorOfRatioByDivision(IntegralType numerator,
155 IntegralType denominator) {
156 return numerator / denominator;
157 }
158
159 template <typename Integer, bool ComputeCeil>
CeilOrFloorOfRatioArithmetic(Integer numerator,Integer denominator)160 static Integer CeilOrFloorOfRatioArithmetic(Integer numerator,
161 Integer denominator) {
162 if (ComputeCeil) {
163 return CeilOfRatioDenomMinusOne(numerator, denominator);
164 } else {
165 return FloorOfRatioByDivision(numerator, denominator);
166 }
167 }
168
TestThatCeilOfRatioDenomMinusOneIsIncorrect(int64 numerator,int64 denominator,int64 expected_error)169 void TestThatCeilOfRatioDenomMinusOneIsIncorrect(int64 numerator,
170 int64 denominator,
171 int64 expected_error) {
172 const int64 correct_result = MathUtil::CeilOfRatio(numerator, denominator);
173 const int64 result_by_denom_minus_one =
174 CeilOfRatioDenomMinusOne(numerator, denominator);
175 EXPECT_EQ(result_by_denom_minus_one + expected_error, correct_result)
176 << "numerator = " << numerator << " denominator = " << denominator
177 << " expected error = " << expected_error
178 << " Actual difference: " << (correct_result - result_by_denom_minus_one);
179 }
180
181 // Here we demonstrate why not to use CeilOfRatioDenomMinusOne
TestThatCeilOfRatioDenomMinusOneIsIncorrect()182 void TestThatCeilOfRatioDenomMinusOneIsIncorrect() {
183 // It does not work with negative values
184 TestThatCeilOfRatioDenomMinusOneIsIncorrect(-1LL, -2LL, -1LL);
185
186 // This would also fail if given kint64max because of signed integer overflow.
187 }
188
TEST(MathUtil,CeilOfRatio)189 TEST(MathUtil, CeilOfRatio) {
190 TestCeilOfRatioUnsigned<uint8>(kuint8max);
191 TestCeilOfRatioUnsigned<uint16>(kuint16max);
192 TestCeilOfRatioUnsigned<uint32>(kuint32max);
193 TestCeilOfRatioUnsigned<uint64>(kuint64max);
194 TestCeilOfRatioSigned<int8>(kint8min, kint8max);
195 TestCeilOfRatioSigned<int16>(kint16min, kint16max);
196 TestCeilOfRatioSigned<int32>(kint32min, kint32max);
197 TestCeilOfRatioSigned<int64>(kint64min, kint64max);
198 #if 0
199 TestThatCeilOfRatioDenomMinusOneIsIncorrect();
200 #endif
201 }
202
203 struct GCDTestCase {
204 unsigned int x;
205 unsigned int y;
206 unsigned int gcd;
207 };
208
TEST(MathUtil,GCD)209 TEST(MathUtil, GCD) {
210 std::vector<GCDTestCase> testcases({
211 {10, 20, 10}, //
212 {27, 8, 1}, //
213 {4, 3, 1}, //
214 {6, 8, 2}, //
215 {5, 0, 5}, //
216 {5, 5, 5}, //
217 {0, 0, 0} //
218 });
219
220 for (const auto& tc : testcases) {
221 EXPECT_EQ(tc.gcd, MathUtil::GCD<uint32>(tc.x, tc.y));
222 EXPECT_EQ(tc.gcd, MathUtil::GCD<uint32>(tc.y, tc.x));
223 EXPECT_EQ(tc.gcd, MathUtil::GCD<uint64>(tc.x, tc.y));
224 EXPECT_EQ(tc.gcd, MathUtil::GCD<uint64>(tc.y, tc.x));
225 }
226
227 const uint64 biggish_prime = 1666666667;
228 EXPECT_EQ(biggish_prime,
229 MathUtil::GCD<uint64>(biggish_prime * 3, biggish_prime * 4));
230 }
231
232 template <typename T>
TestOneIPowN()233 void TestOneIPowN() {
234 const T one{1};
235 for (int i = 0; i < 1024; ++i) {
236 // Computations are exact.
237 EXPECT_EQ(MathUtil::IPow(one, i), one);
238 }
239 }
240
241 template <typename T>
TestTwoIPowN()242 void TestTwoIPowN() {
243 int limit = std::is_integral<T>::value ? std::numeric_limits<T>::digits : 63;
244 for (int i = 0; i < limit; ++i) {
245 // Computations are exact.
246 EXPECT_EQ(MathUtil::IPow(T{2}, i), static_cast<T>(1ull << i));
247 }
248 }
249
250 template <typename T>
TestFloatIPow(const int max_exponent,const T start,const T end,const T step)251 void TestFloatIPow(const int max_exponent, const T start, const T end,
252 const T step) {
253 for (T f = start; f < end; f += step) {
254 for (int i = 0; i < max_exponent; ++i) {
255 EXPECT_FLOAT_EQ(MathUtil::IPow(f, i), pow(f, i));
256 }
257 }
258 }
259
TEST(MathUtil,IPow)260 TEST(MathUtil, IPow) {
261 TestOneIPowN<double>();
262 TestOneIPowN<float>();
263 TestOneIPowN<int>();
264 TestOneIPowN<int64>();
265 TestTwoIPowN<double>();
266 TestTwoIPowN<float>();
267 TestTwoIPowN<int>();
268 TestTwoIPowN<int64>();
269
270 EXPECT_EQ(MathUtil::IPow(3, 0), 1);
271 EXPECT_EQ(MathUtil::IPow(3, 1), 3);
272 EXPECT_EQ(MathUtil::IPow(3, 2), 9);
273 EXPECT_EQ(MathUtil::IPow(3, 3), 27);
274 EXPECT_EQ(MathUtil::IPow(3, 4), 81);
275 EXPECT_EQ(MathUtil::IPow(3, 5), 243);
276
277 TestFloatIPow<float>(13, -16.0f, 16.0f, 1.0f / 8);
278 TestFloatIPow<double>(13, -16.0, 16.0, 1.0 / 8);
279
280 TestFloatIPow<float>(13, -1.0f / (1 << 12), -1.0f / (1 << 12),
281 1.0f / (1 << 16));
282 TestFloatIPow<double>(13, -1.0 / (1 << 12), -1.0 / (1 << 12),
283 1.0 / (1 << 16));
284 }
285
TEST(MathUtil,IPowEdgeCases)286 TEST(MathUtil, IPowEdgeCases) {
287 constexpr const double kInf = std::numeric_limits<double>::infinity();
288
289 EXPECT_EQ(MathUtil::IPow(-12345.0, 79), -kInf);
290 EXPECT_EQ(MathUtil::IPow(-12345.0, 80), +kInf);
291
292 // The semantics of the edge cases that follow are defined in the standard:
293 // http://en.cppreference.com/w/cpp/numeric/math/pow for a summary.
294
295 // 1 - These edge cases apply.
296 // pow(+0, exp), where exp is a positive odd integer, returns +0
297 EXPECT_EQ(MathUtil::IPow(+0.0, 3), +0.0);
298 // pow(-0, exp), where exp is a positive odd integer, returns -0
299 EXPECT_EQ(MathUtil::IPow(-0.0, 3), -0.0);
300 // pow(±0, exp), where exp is positive non-integer or a positive even integer,
301 // returns +0
302 EXPECT_EQ(MathUtil::IPow(+0.0, 42), +0.0);
303 EXPECT_EQ(MathUtil::IPow(-0.0, 42), +0.0);
304 // pow(base, ±0) returns 1 for any base, even when base is NaN
305 EXPECT_EQ(MathUtil::IPow(-kInf, 0.0), 1.0);
306 EXPECT_EQ(MathUtil::IPow(-2.0, 0.0), 1.0);
307 EXPECT_EQ(MathUtil::IPow(-1.0, 0.0), 1.0);
308 EXPECT_EQ(MathUtil::IPow(-0.0, 0.0), 1.0);
309 EXPECT_EQ(MathUtil::IPow(+0.0, 0.0), 1.0);
310 EXPECT_EQ(MathUtil::IPow(+1.0, 0.0), 1.0);
311 EXPECT_EQ(MathUtil::IPow(+2.0, 0.0), 1.0);
312 EXPECT_EQ(MathUtil::IPow(+kInf, 0.0), 1.0);
313 EXPECT_EQ(MathUtil::IPow(std::numeric_limits<double>::quiet_NaN(), 0.0), 1.0);
314 // pow(-∞, exp) returns -∞ if exp is a positive odd integer
315 EXPECT_EQ(MathUtil::IPow(-kInf, 43), -kInf);
316 // pow(-∞, exp) returns +∞ if exp is a positive non-integer or even integer
317 EXPECT_EQ(MathUtil::IPow(-kInf, 42), +kInf);
318 // pow(+∞, exp) returns +∞ for any positive exp
319 EXPECT_EQ(MathUtil::IPow(+kInf, 42), +kInf);
320 EXPECT_EQ(MathUtil::IPow(+kInf, 43), +kInf);
321
322 // 2 - These do not apply due to the restricted exp range.
323 // pow(+0, exp), where exp is a negative odd integer, returns +∞ and raises
324 // FE_DIVBYZERO pow(-0, exp), where exp is a negative odd integer, returns -∞
325 // and raises FE_DIVBYZERO pow(±0, exp), where exp is negative, finite, and is
326 // an even integer or a non-integer, returns +∞ and raises FE_DIVBYZERO
327 // pow(-1, ±∞) returns 1
328 // pow(+1, exp) returns 1 for any exp, even when exp is NaN
329 // pow(±0, -∞) returns +∞ and may raise FE_DIVBYZERO
330 // pow(base, exp) returns NaN and raises FE_INVALID if base is finite and
331 // negative and exp is finite and non-integer. pow(base, -∞) returns +∞ for
332 // any |base|<1 pow(base, -∞) returns +0 for any |base|>1 pow(base, +∞)
333 // returns +0 for any |base|<1 pow(base, +∞) returns +∞ for any |base|>1
334 // pow(-∞, exp) returns -0 if exp is a negative odd integer
335 // pow(-∞, exp) returns +0 if exp is a negative non-integer or even integer
336 // pow(+∞, exp) returns +0 for any negative exp
337 }
338
339 } // namespace
340 } // namespace tensorflow
341