1// Copyright 2001 and onwards Google Inc.
2//
3// This class is intended to contain a collection of useful (static)
4// mathematical functions, properly coded (by consulting numerical
5// recipes or another authoritative source first).
6
7#ifndef UTIL_MATH_MATHUTIL_H__
8#define UTIL_MATH_MATHUTIL_H__
9
10#include <math.h>
11#include <algorithm>
12using std::min;
13using std::max;
14using std::swap;
15using std::reverse;
16
17#include <vector>
18using std::vector;
19
20#include "base/basictypes.h"
21#include "base/logging.h"
22#include "util/math/mathlimits.h"
23
24// Returns the sign of x:
25// -1 if x < 0,
26// +1 if x > 0,
27// 0 if x = 0.
28// Consider instead using MathUtil::Sign below for readability
29// and floating-point correctness.
30template <class T>
31inline T sgn(const T x) {
32 return (x == 0 ? 0 : (x < 0 ? -1 : 1));
33}
34
35// ========================================================================= //
36
37class MathUtil {
38 public:
39
40 // Return type of RealRootsForQuadratic (below). The enum values are
41 // chosen to be sensible if converted to bool or int, and should not be
42 // changed lightly.
43 enum QuadraticRootType {kNoRealRoots = 0, kAmbiguous = 1, kTwoRealRoots = 2};
44
45 // Returns the QuadraticRootType of the equation a * x^2 + b * x + c = 0.
46 // Normal cases are kNoRealRoots, in which case *r1 and *r2 are not
47 // changed; and kTwoRealRoots, in which case the root(s) are placed in
48 // *r1 and *r2, order not specified. The kAmbiguous return value
49 // indicates that the disciminant is zero, within floating-point error
50 // (i.e. that changing an input by epsilon<double> would change the sign
51 // of the discriminant). The resulting roots are equal, as if the
52 // discriminant were exactly zero.
53 //
54 // A special case occurs when a==0; see DegenerateQuadraticRoots().
55 // See also QuadraticIsAmbiguous() and RealQuadraticRoots().
56 static inline QuadraticRootType RealRootsForQuadratic(long double a,
57 long double b,
58 long double c,
59 long double *r1,
60 long double *r2) {
61 // Deal with degenerate cases where leading coefficients vanish.
62 if (a == 0.0) {
63 return DegenerateQuadraticRoots(b, c, r1, r2);
64 }
65
66 // General case: the quadratic formula, rearranged for greater numerical
67 // stability.
68
69 // If the discriminant is zero to numerical precision, regardless of
70 // sign, treat it as zero and return kAmbiguous. We use the double
71 // rather than long double value for epsilon because in practice inputs
72 // are generally calculated in double precision.
73 const long double discriminant = QuadraticDiscriminant(a, b, c);
74 if (QuadraticIsAmbiguous(a, b, c, discriminant,
75 MathLimits<double>::kEpsilon)) {
76 *r2 = *r1 = -b / 2 / a; // The quadratic is (2*a*x + b)^2 = 0.
77 return kAmbiguous;
78 }
79
80 if (discriminant < 0) {
81 // The discriminant is definitely negative so there are no real roots.
82 return kNoRealRoots;
83 }
84
85 RealQuadraticRoots(a, b, c, discriminant, r1, r2);
86 return kTwoRealRoots;
87 }
88
89 // Returns the discriminant of the quadratic equation a * x^2 + b * x + c = 0.
90 static inline long double QuadraticDiscriminant(long double a,
91 long double b,
92 long double c) {
93 return b * b - 4 * a * c;
94 }
95
96 // Returns true if the discriminant is zero within floating-point error,
97 // in the sense that changing one of the coefficients by epsilon (e.g. a
98 // -> a + a*epsilon) could change the sign of the discriminant. [When the
99 // discriminant is exactly 0 the quadratic is (2*a*x + b)^2 = 0 and the
100 // root is - b / (2*a).]
101 static inline bool QuadraticIsAmbiguous(long double a,
102 long double b,
103 long double c,
104 long double discriminant,
105 long double epsilon) {
106 // Discriminants below kTolerance in absolute value are considered zero
107 // because changing the final bit of one of the inputs can change the
108 // sign of the discriminant.
109 const double kTolerance = epsilon * max(fabs(2 * b * b), fabs(4 * a * c));
110 return (fabs(discriminant) <= kTolerance);
111 }
112
113 // Returns in *r1 and *r2 the roots of a "normal" quadratic equation
114 // whose discriminant (b*b - 4*a*c) is known and positive. Preconditions
115 // (will DCHECK and return false if not satisfied): a != 0, discriminant > 0.
116 static inline bool RealQuadraticRoots(long double a,
117 long double b,
118 long double c,
119 long double discriminant,
120 long double *r1,
121 long double *r2) {
122 if (discriminant <= 0 || a == 0) {
123 // A case that should have been excluded by the caller.
124 DCHECK(false);
125 return false;
126 }
127
128 // The discriminant is positive so there are two distinct roots.
129 // According to Numerical Recipes (p. 184), it would be a mistake to
130 // solve for the roots using
131 //
132 // r1 = 2c / (-b + sqrt(b^2 - 4ac)),
133 // r2 = 2c / (-b - sqrt(b^2 - 4ac)).
134 //
135 // If a*c is small, then one of the roots above will involve the
136 // subtraction of b from a very nearly equal quantity (the discriminant),
137 // producing a very inaccurate root. Avoid the risk of cancellation with
138 // the following rearrangement. (Note we don't use sgn(b) because we
139 // need sgn(0) = +1 or -1.)
140 long double const q = -0.5 *
141 (b + ((b >= 0) ? sqrt(discriminant) : -sqrt(discriminant)));
142 *r1 = q / a; // If a is very small this produces +/- HUGE_VAL.
143 *r2 = c / q; // q cannot be too small.
144 return true;
145 }
146
147 // Returns the root of the degenerate quadratic equation b * x + c = 0,
148 // following the interface of RealRootsForQuadratic. To be consistent
149 // with that function as a->0, the degenerate quadratic is considered to
150 // have two real roots, one of which is +/- HUGE_VAL and one of which is
151 // -c / b. If both a and b are 0, so the equation is c = 0, the response
152 // is kNoRealRoots if c != 0 or kAmbiguous if c == 0 (since the
153 // discriminant is zero).
154 static QuadraticRootType DegenerateQuadraticRoots(long double b,
155 long double c,
156 long double *r1,
157 long double *r2);
158
159 // Solves for the real roots of x^3+ax^2+bx+c=0, returns true iff
160 // all three are real, in which case the roots are stored (in any
161 // order) in r1, r2, r3; otherwise, exactly one real root exists and
162 // it is stored in r1.
163 static bool RealRootsForCubic(long double a,
164 long double b,
165 long double c,
166 long double *r1,
167 long double *r2,
168 long double *r3);
169
170
171 // ----------------------------------------------------------------------
172 // Sigmoid
173 // A sigmoid function is a differentiably s curve that ranges between
174 // 0 and 1:
175 // f(x) = 1/(1+e^(-lambda x))
176 // --------------------------------------------------------------------
177 static double Sigmoid(double x, double lambda) {
178 return 1/(1+exp(-lambda*x));
179 }
180
181 // ----------------------------------------------------------------------
182 // InverseSigmoid
183 // Inverts Sigmoid such that InverseSigmoid(Sigmoid(x)) == x for all x
184 // Note that all inputs must be in (-1, 1)
185 // ----------------------------------------------------------------------
186 static double InverseSigmoid(double const x, double const lambda) {
187 return -log((1.0 / x - 1)) / lambda;
188 }
189
190 // ----------------------------------------------------------------------
191 // Sigmoid2
192 // A nicer way of specifying a sigmoid. A sigmoid is a smooth s curve
193 // that ranges from 0 to 1. We describe a sigmoid with three values:
194 //
195 // start: the x value at which f(x) = tolerance
196 // finish: the x value at which f(x) = 1-tolerance
197 //
198 // So if we was a smoothly transitioning function from, say, x=1 to
199 // x=10 with the property that anything outside the domain [1, 10]
200 // will still be within 10% of either f(1) or f(10) then we set: start
201 // = 1 finish = 10 tolerance = 0.1
202 // --------------------------------------------------------------------
203 static double Sigmoid2(double x, double start_x,
204 double finish_x, double tolerance) {
205 DCHECK_GT(tolerance, 0);
206 DCHECK_LT(tolerance, 1);
207 DCHECK_NE(finish_x - start_x, 0);
208 double lambda = log((1-tolerance)/tolerance)*2/(finish_x - start_x);
209 return Sigmoid(x - 0.5 * (start_x + finish_x), lambda);
210 }
211
212 // Returns the greatest common divisor of two unsigned integers x and y
213 static unsigned int GCD(unsigned int x, unsigned int y) {
214 while (y != 0) {
215 unsigned int r = x % y;
216 x = y;
217 y = r;
218 }
219 return x;
220 }
221
222 // Returns the greatest common divisor of two unsigned integers x and y,
223 // and assigns a, and b such that a*x + b*y = gcd(x, y).
224 static unsigned int ExtendedGCD(unsigned int x, unsigned int y,
225 int* a, int* b);
226
227 // Returns the least common multiple of two unsigned integers. Returns zero
228 // if either is zero.
229 static unsigned int LeastCommonMultiple(unsigned int a, unsigned int b) {
230 if (a > b) {
231 return (a / MathUtil::GCD(a, b)) * b;
232 } else if (a < b) {
233 return (b / MathUtil::GCD(b, a)) * a;
234 } else {
235 return a;
236 }
237 }
238
239 // Converts a non-zero double value representing an odds into its
240 // probability value.
241 static double OddsToProbability(double odds) {
242 DCHECK_GE(odds, 0.0);
243 return odds / (1.0 + odds);
244 }
245
246 // Converts a probability with range [0-1.0) into its odds value.
247 static double ProbabilityToOdds(double prob) {
248 DCHECK_GE(prob, 0.0);
249 DCHECK_LT(prob, 1.0);
250 return prob / (1.0 - prob);
251 }
252
253 // --------------------------------------------------------------------
254 // ShardsToRead
255 // Resharding helper. Suppose we have N input shards and M output
256 // shards sharded by modulo of the same hash function. If we want
257 // to write a subset of the output shards, which input shards should
258 // we read?
259 //
260 // Inputs:
261 // shards_to_write gives the desired subset of the M output shards.
262 // shards_to_read gives the number N of the input shards.
263 // Outputs:
264 // shards_to_read gives the subset of the N input shards to read.
265 // --------------------------------------------------------------------
266 static void ShardsToRead(const vector<bool>& shards_to_write,
267 vector<bool>* shards_to_read);
268
269 // --------------------------------------------------------------------
270 // Round, IntRound
271 // These functions round a floating-point number to an integer. They
272 // work for positive or negative numbers.
273 //
274 // Values that are halfway between two integers may be rounded up or
275 // down, for example IntRound(0.5) == 0 and IntRound(1.5) == 2. This
276 // allows these functions to be implemented efficiently on Intel
277 // processors (see the template specializations at the bottom of this
278 // file). You should not use these functions if you care about which
279 // way such half-integers are rounded.
280 //
281 // Example usage:
282 // double y, z;
283 // int x = IntRound(y + 3.7);
284 // int64 b = Round<int64>(0.3 * z);
285 //
286 // Note that the floating-point template parameter is typically inferred
287 // from the argument type, i.e. there is no need to specify it explicitly.
288 // --------------------------------------------------------------------
289 template <class IntOut, class FloatIn>
290 static IntOut Round(FloatIn x) {
291 COMPILE_ASSERT(!MathLimits<FloatIn>::kIsInteger, FloatIn_is_integer);
292 COMPILE_ASSERT(MathLimits<IntOut>::kIsInteger, IntOut_is_not_integer);
293
294 // We don't use sgn(x) below because there is no need to distinguish the
295 // (x == 0) case. Also note that there are specialized faster versions
296 // of this function for Intel processors at the bottom of this file.
297 return static_cast<IntOut>(x < 0 ? (x - 0.5) : (x + 0.5));
298 }
299
300 // Example usage: IntRound(3.6) (no need for IntRound<double>(3.6)).
301 template <class FloatIn>
302 static int IntRound(FloatIn x) { return Round<int>(x); }
303
304 // --------------------------------------------------------------------
305 // FastIntRound, FastInt64Round
306 // Fast routines for converting floating-point numbers to integers.
307 //
308 // These routines are approximately 6 times faster than the default
309 // implementation of IntRound() on Intel processors (12 times faster on
310 // the Pentium 3). They are also more than 5 times faster than simply
311 // casting a "double" to an "int" using static_cast<int>. This is
312 // because casts are defined to truncate towards zero, which on Intel
313 // processors requires changing the rounding mode and flushing the
314 // floating-point pipeline (unless programs are compiled specifically
315 // for the Pentium 4, which has a new instruction to avoid this).
316 //
317 // Numbers that are halfway between two integers may be rounded up or
318 // down. This is because the conversion is done using the default
319 // rounding mode, which rounds towards the closest even number in case
320 // of ties. So for example, FastIntRound(0.5) == 0, but
321 // FastIntRound(1.5) == 2. These functions should only be used with
322 // applications that don't care about which way such half-integers are
323 // rounded.
324 //
325 // There are template specializations of Round() which call these
326 // functions (for "int" and "int64" only), but it's safer to call them
327 // directly.
328 //
329 // This functions are equivalent to lrint() and llrint() as defined in
330 // the ISO C99 standard. Unfortunately this standard does not seem to
331 // widely adopted yet and these functions are not available by default.
332 // --------------------------------------------------------------------
333
334 static int32 FastIntRound(double x) {
335 // This function is not templatized because gcc doesn't seem to be able
336 // to deal with inline assembly code in templatized functions, and there
337 // is no advantage to passing an argument type of "float" on Intel
338 // architectures anyway.
339
340#if defined __GNUC__ && (defined __i386__ || defined __SSE2__)
341#if defined __SSE2__
342 // SSE2.
343 int32 result;
344 __asm__ __volatile__
345 ("cvtsd2si %1, %0"
346 : "=r" (result) // Output operand is a register
347 : "x" (x)); // Input operand is an xmm register
348 return result;
349#elif defined __i386__
350 // FPU stack. Adapted from /usr/include/bits/mathinline.h.
351 int32 result;
352 __asm__ __volatile__
353 ("fistpl %0"
354 : "=m" (result) // Output operand is a memory location
355 : "t" (x) // Input operand is top of FP stack
356 : "st"); // Clobbers (pops) top of FP stack
357 return result;
358#endif // if defined __x86_64__ || ...
359#else
360 return Round<int32, double>(x);
361#endif // if defined __GNUC__ && ...
362 }
363
364 static int64 FastInt64Round(double x) {
365#if defined __GNUC__ && (defined __i386__ || defined __x86_64__)
366#if defined __x86_64__
367 // SSE2.
368 int64 result;
369 __asm__ __volatile__
370 ("cvtsd2si %1, %0"
371 : "=r" (result) // Output operand is a register
372 : "x" (x)); // Input operand is an xmm register
373 return result;
374#elif defined __i386__
375 // There is no CVTSD2SI in i386 to produce a 64 bit int, even with SSE2.
376 // FPU stack. Adapted from /usr/include/bits/mathinline.h.
377 int64 result;
378 __asm__ __volatile__
379 ("fistpll %0"
380 : "=m" (result) // Output operand is a memory location
381 : "t" (x) // Input operand is top of FP stack
382 : "st"); // Clobbers (pops) top of FP stack
383 return result;
384#endif // if defined __i386__
385#else
386 return Round<int64, double>(x);
387#endif // if defined __GNUC__ && ...
388 }
389
390 // Return Not a Number.
391 // Consider instead using MathLimits<double>::kNaN directly.
392 static double NaN() { return MathLimits<double>::kNaN; }
393
394 // the sine cardinal function
395 static double Sinc(double x) {
396 if (fabs(x) < 1E-8) return 1.0;
397 return sin(x) / x;
398 }
399
400 // Returns an approximation An for the n-th element of the harmonic
401 // serices Hn = 1 + ... + 1/n. Sets error e such that |An-Hn| < e.
402 static double Harmonic(int64 n, double *e);
403
404 // Returns Stirling's Approximation for log(n!) which has an error
405 // of at worst 1/(1260*n^5).
406 static double Stirling(double n);
407
408 // Returns the log of the binomial coefficient C(n, k), known in the
409 // vernacular as "N choose K". Why log? Because the integer number
410 // for non-trivial N and K would overflow.
411 // Note that if k > 15, this uses Stirling's approximation of log(n!).
412 // The relative error is about 1/(1260*k^5) (which is 7.6e-10 when k=16).
413 static double LogCombinations(int n, int k);
414
415 // Rounds "f" to the nearest float which has its last "bits" bits of
416 // the mantissa set to zero. This rounding will introduce a
417 // fractional error of at most 2^(bits - 24). Useful for values
418 // stored in compressed files, when super-accurate numbers aren't
419 // needed and the random-looking low-order bits foil compressors.
420 // This routine should be really fast when inlined with "bits" set
421 // to a constant.
422 // Precondition: 1 <= bits <= 23, f != NaN
423 static float RoundOffBits(const float f, const int bits) {
424 const int32 f_rep = bit_cast<int32>(f);
425
426 // Set low-order "bits" bits to zero.
427 int32 g_rep = f_rep & ~((1 << bits) - 1);
428
429 // Round mantissa up if we need to. Note that we do round-to-even,
430 // a.k.a. round-up-if-odd.
431 const int32 lowbits = f_rep & ((1 << bits) - 1);
432 if (lowbits > (1 << (bits - 1)) ||
433 (lowbits == (1 << (bits - 1)) && (f_rep & (1 << bits)))) {
434 g_rep += (1 << bits);
435 // NOTE: overflow does a really nice thing here - if all the
436 // rest of the mantissa bits are 1, the carry carries over into
437 // the exponent and increments it by 1, which is exactly what we
438 // want. It even gets to +/-INF properly.
439 }
440 return bit_cast<float>(g_rep);
441 }
442 // Same, but for doubles. 1 <= bits <= 52, error at most 2^(bits - 53).
443 static double RoundOffBits(const double f, const int bits) {
444 const int64 f_rep = bit_cast<int64>(f);
445 int64 g_rep = f_rep & ~((1LL << bits) - 1);
446 const int64 lowbits = f_rep & ((1LL << bits) - 1);
447 if (lowbits > (1LL << (bits - 1)) ||
448 (lowbits == (1LL << (bits - 1)) && (f_rep & (1LL << bits)))) {
449 g_rep += (1LL << bits);
450 }
451 return bit_cast<double>(g_rep);
452 }
453
454 // Largest of two values.
455 // Works correctly for special floating point values.
456 // Note: 0.0 and -0.0 are not differentiated by Max (Max(0.0, -0.0) is -0.0),
457 // which should be OK because, although they (can) have different
458 // bit representation, they are observably the same when examined
459 // with arithmetic and (in)equality operators.
460 template<typename T>
461 static T Max(const T x, const T y) {
462 return MathLimits<T>::IsNaN(x) || x > y ? x : y;
463 }
464
465 // Smallest of two values
466 // Works correctly for special floating point values.
467 // Note: 0.0 and -0.0 are not differentiated by Min (Min(-0.0, 0.0) is 0.0),
468 // which should be OK: see the comment for Max above.
469 template<typename T>
470 static T Min(const T x, const T y) {
471 return MathLimits<T>::IsNaN(x) || x < y ? x : y;
472 }
473
474 // Absolute value of x
475 // Works correctly for unsigned types and
476 // for special floating point values.
477 // Note: 0.0 and -0.0 are not differentiated by Abs (Abs(0.0) is -0.0),
478 // which should be OK: see the comment for Max above.
479 template<typename T>
480 static T Abs(const T x) {
481 return x > 0 ? x : -x;
482 }
483
484 // The sign of x
485 // (works for unsigned types and special floating point values as well):
486 // -1 if x < 0,
487 // +1 if x > 0,
488 // 0 if x = 0.
489 // nan if x is nan.
490 template<typename T>
491 static T Sign(const T x) {
492 return MathLimits<T>::IsNaN(x) ? x : (x == 0 ? 0 : (x > 0 ? 1 : -1));
493 }
494
495 // Returns the square of x
496 template <typename T>
497 static T Square(const T x) {
498 return x * x;
499 }
500
501 // Absolute value of the difference between two numbers.
502 // Works correctly for signed types and special floating point values.
503 template<typename T>
504 static typename MathLimits<T>::UnsignedType AbsDiff(const T x, const T y) {
505 return x > y ? x - y : y - x;
506 }
507
508 // CAVEAT: Floating point computation instability for x86 CPUs
509 // can frequently stem from the difference of when floating point values
510 // are transferred from registers to memory and back,
511 // which can depend simply on the level of optimization.
512 // The reason is that the registers use a higher-precision representation.
513 // Hence, instead of relying on approximate floating point equality below
514 // it might be better to reorganize the code with volatile modifiers
515 // for the floating point variables so as to control when
516 // the loss of precision occurs.
517
518 // If two (usually floating point) numbers are within a certain
519 // absolute margin of error.
520 // NOTE: this "misbehaves" is one is trying to capture provisons for errors
521 // that are relative, i.e. larger if the numbers involved are larger.
522 // Consider using WithinFraction or WithinFractionOrMargin below.
523 //
524 // This and other Within* NearBy* functions below
525 // work correctly for signed types and special floating point values.
526 template<typename T>
527 static bool WithinMargin(const T x, const T y, const T margin) {
528 DCHECK_GE(margin, 0);
529 // this is a little faster than x <= y + margin && x >= y - margin
530 return AbsDiff(x, y) <= margin;
531 }
532
533 // If two (usually floating point) numbers are within a certain
534 // fraction of their magnitude.
535 // CAVEAT: Although this works well if computation errors are relative
536 // both for large magnitude numbers > 1 and for small magnitude numbers < 1,
537 // zero is never within a fraction of any
538 // non-zero finite number (fraction is required to be < 1).
539 template<typename T>
540 static bool WithinFraction(const T x, const T y, const T fraction);
541
542 // If two (usually floating point) numbers are within a certain
543 // fraction of their magnitude or within a certain absolute margin of error.
544 // This is the same as the following but faster:
545 // WithinFraction(x, y, fraction) || WithinMargin(x, y, margin)
546 // E.g. WithinFraction(0.0, 1e-10, 1e-5) is false but
547 // WithinFractionOrMargin(0.0, 1e-10, 1e-5, 1e-5) is true.
548 template<typename T>
549 static bool WithinFractionOrMargin(const T x, const T y,
550 const T fraction, const T margin);
551
552 // NearBy* functions below are geared as replacements for CHECK_EQ()
553 // over floating-point numbers.
554
555 // Same as WithinMargin(x, y, MathLimits<T>::kStdError)
556 // Works as == for integer types.
557 template<typename T>
558 static bool NearByMargin(const T x, const T y) {
559 return AbsDiff(x, y) <= MathLimits<T>::kStdError;
560 }
561
562 // Same as WithinFraction(x, y, MathLimits<T>::kStdError)
563 // Works as == for integer types.
564 template<typename T>
565 static bool NearByFraction(const T x, const T y) {
566 return WithinFraction(x, y, MathLimits<T>::kStdError);
567 }
568
569 // Same as WithinFractionOrMargin(x, y, MathLimits<T>::kStdError,
570 // MathLimits<T>::kStdError)
571 // Works as == for integer types.
572 template<typename T>
573 static bool NearByFractionOrMargin(const T x, const T y) {
574 return WithinFractionOrMargin(x, y, MathLimits<T>::kStdError,
575 MathLimits<T>::kStdError);
576 }
577
578 // Tests whether two values are close enough to each other to be considered
579 // equal. This function is intended to be used mainly as a replacement for
580 // equality tests of floating point values in CHECK()s, and as a replacement
581 // for equality comparison against golden files. It is the same as == for
582 // integer types. The purpose of AlmostEquals() is to avoid false positive
583 // error reports in unit tests and regression tests due to minute differences
584 // in floating point arithmetic (for example, due to a different compiler).
585 //
586 // We cannot use simple equality to compare floating point values
587 // because floating point expressions often accumulate inaccuracies, and
588 // new compilers may introduce further variations in the values.
589 //
590 // Two values x and y are considered "almost equals" if:
591 // (a) Both values are very close to zero: x and y are in the range
592 // [-standard_error, standard_error]
593 // Normal calculations producing these values are likely to be dealing
594 // with meaningless residual values.
595 // -or-
596 // (b) The difference between the values is small:
597 // abs(x - y) <= standard_error
598 // -or-
599 // (c) The values are finite and the relative difference between them is
600 // small:
601 // abs (x - y) <= standard_error * max(abs(x), abs(y))
602 // -or-
603 // (d) The values are both positive infinity or both negative infinity.
604 //
605 // Cases (b) and (c) are the same as MathUtils::NearByFractionOrMargin(x, y),
606 // for finite values of x and y.
607 //
608 // standard_error is the corresponding MathLimits<T>::kStdError constant.
609 // It is equivalent to 5 bits of mantissa error. See
610 // google3/util/math/mathlimits.cc.
611 //
612 // Caveat:
613 // AlmostEquals() is not appropriate for checking long sequences of
614 // operations where errors may cascade (like extended sums or matrix
615 // computations), or where significant cancellation may occur
616 // (e.g., the expression (x+y)-x, where x is much larger than y).
617 // Both cases may produce errors in excess of standard_error.
618 // In any case, you should not test the results of calculations which have
619 // not been vetted for possible cancellation errors and the like.
620 template<typename T>
621 static bool AlmostEquals(const T x, const T y) {
622 if (x == y) // Covers +inf and -inf, and is a shortcut for finite values.
623 return true;
624 if (!MathLimits<T>::IsFinite(x) || !MathLimits<T>::IsFinite(y))
625 return false;
626
627 if (MathUtil::Abs<T>(x) <= MathLimits<T>::kStdError &&
628 MathUtil::Abs<T>(y) <= MathLimits<T>::kStdError)
629 return true;
630
631 return MathUtil::NearByFractionOrMargin<T>(x, y);
632 }
633
634 // Returns the clamped value to be between low and high inclusively.
635 template<typename T>
636 static const T& Clamp(const T& low, const T& high, const T& value) {
637 return std::max(low, std::min(value, high));
638 }
639
640 // Clamps value to be between min and max inclusively.
641 template<typename T>
642 static void ClampValue(const T& low, const T& high, T* value) {
643 *value = Clamp(low, high, *value);
644 }
645};
646
647// ========================================================================= //
648
649#if (defined __i386__ || defined __x86_64__) && defined __GNUC__
650
651// We define template specializations of Round() to get the more efficient
652// Intel versions when possible. Note that gcc does not currently support
653// partial specialization of templatized functions.
654
655template<>
656inline int32 MathUtil::Round<int32, double>(double x) {
657 return FastIntRound(x);
658}
659
660template<>
661inline int32 MathUtil::Round<int32, float>(float x) {
662 return FastIntRound(x);
663}
664
665template<>
666inline int64 MathUtil::Round<int64, double>(double x) {
667 return FastInt64Round(x);
668}
669
670template<>
671inline int64 MathUtil::Round<int64, float>(float x) {
672 return FastInt64Round(x);
673}
674
675#endif
676
677template<typename T>
678bool MathUtil::WithinFraction(const T x, const T y, const T fraction) {
679 // not just "0 <= fraction" to fool the compiler for unsigned types
680 DCHECK((0 < fraction || 0 == fraction) && fraction < 1);
681
682 // Template specialization will convert the if() condition to a constant,
683 // which will cause the compiler to generate code for either the "if" part
684 // or the "then" part. In this way we avoid a compiler warning
685 // about a potential integer overflow in crosstool v12 (gcc 4.3.1).
686 if (MathLimits<T>::kIsInteger) {
687 return x == y;
688 } else {
689 // IsFinite checks are to make kPosInf and kNegInf not within fraction
690 return (MathLimits<T>::IsFinite(x) || MathLimits<T>::IsFinite(y)) &&
691 (AbsDiff(x, y) <= fraction * Max(Abs(x), Abs(y)));
692 }
693}
694
695template<typename T>
696bool MathUtil::WithinFractionOrMargin(const T x, const T y,
697 const T fraction, const T margin) {
698 // not just "0 <= fraction" to fool the compiler for unsigned types
699 DCHECK((0 < fraction || 0 == fraction) && fraction < 1 && margin >= 0);
700
701 // Template specialization will convert the if() condition to a constant,
702 // which will cause the compiler to generate code for either the "if" part
703 // or the "then" part. In this way we avoid a compiler warning
704 // about a potential integer overflow in crosstool v12 (gcc 4.3.1).
705 if (MathLimits<T>::kIsInteger) {
706 return x == y;
707 } else {
708 // IsFinite checks are to make kPosInf and kNegInf not within fraction
709 return (MathLimits<T>::IsFinite(x) || MathLimits<T>::IsFinite(y)) &&
710 (AbsDiff(x, y) <= Max(margin, fraction * Max(Abs(x), Abs(y))));
711 }
712}
713
714#endif // UTIL_MATH_MATHUTIL_H__
715