1/* Copyright 2015 The TensorFlow Authors. All Rights Reserved.
2
3Licensed under the Apache License, Version 2.0 (the "License");
4you may not use this file except in compliance with the License.
5You may obtain a copy of the License at
6
7    http://www.apache.org/licenses/LICENSE-2.0
8
9Unless required by applicable law or agreed to in writing, software
10distributed under the License is distributed on an "AS IS" BASIS,
11WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12See the License for the specific language governing permissions and
13limitations 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
27namespace tensorflow {
28namespace {
29
30// Number of arguments for each test of the CeilOrRatio method
31const int kNumTestArguments = 4;
32
33template <typename IntegralType, typename TestDataType>
34void 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
61template <typename UnsignedIntegralType>
62void 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
88template <typename SignedInteger>
89void 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
144template <typename IntegralType>
145static 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
153template <typename IntegralType>
154static IntegralType FloorOfRatioByDivision(IntegralType numerator,
155                                           IntegralType denominator) {
156  return numerator / denominator;
157}
158
159template <typename Integer, bool ComputeCeil>
160static 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
169void 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
182void 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
189TEST(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
203struct GCDTestCase {
204  unsigned int x;
205  unsigned int y;
206  unsigned int gcd;
207};
208
209TEST(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
232template <typename T>
233void 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
241template <typename T>
242void 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
250template <typename T>
251void 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
260TEST(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
286TEST(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