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