1// Copyright 2008 Google Inc. All Rights Reserved.
2
3#include "util/math/mathutil.h"
4#include <vector>
5using std::vector;
6
7#include "base/integral_types.h"
8#include "base/logging.h"
9
10MathUtil::QuadraticRootType MathUtil::DegenerateQuadraticRoots(
11 long double b,
12 long double c,
13 long double *r1,
14 long double *r2) {
15 // This degenerate quadratic is really a linear equation b * x = -c.
16 if (b == 0.0) {
17 // The equation is constant, c == 0.
18 if (c == 0.0) {
19 // Quadratic equation is 0==0; treat as ambiguous, as if a==epsilon.
20 *r1 = *r2 = 0.0;
21 return kAmbiguous;
22 }
23 return kNoRealRoots;
24 }
25 // The linear equation has a single root at x = -c / b, not a double
26 // one. Respond as if a==epsilon: The other root is at "infinity",
27 // which we signal with HUGE_VAL so that the behavior stays consistent
28 // as a->0.
29 *r1 = -c / b;
30 *r2 = HUGE_VAL;
31 return kTwoRealRoots;
32}
33
34bool MathUtil::RealRootsForCubic(long double const a,
35 long double const b,
36 long double const c,
37 long double *const r1,
38 long double *const r2,
39 long double *const r3) {
40 // According to Numerical Recipes (pp. 184-5), what
41 // follows is an arrangement of computations to
42 // compute the roots of a cubic that minimizes
43 // roundoff error (as pointed out by A.J. Glassman).
44
45 long double const a_squared = a*a, a_third = a/3.0, b_tripled = 3.0*b;
46 long double const Q = (a_squared - b_tripled) / 9.0;
47 long double const R = (2.0*a_squared*a - 3.0*a*b_tripled + 27.0*c) / 54.0;
48
49 long double const R_squared = R*R;
50 long double const Q_cubed = Q*Q*Q;
51 long double const root_Q = sqrt(Q);
52
53 if (R_squared < Q_cubed) {
54 long double const two_pi_third = 2.0 * M_PI / 3.0;
55 long double const theta_third = acos(R / sqrt(Q_cubed)) / 3.0;
56 long double const minus_two_root_Q = -2.0 * root_Q;
57
58 *r1 = minus_two_root_Q * cos(theta_third) - a_third;
59 *r2 = minus_two_root_Q * cos(theta_third + two_pi_third) - a_third;
60 *r3 = minus_two_root_Q * cos(theta_third - two_pi_third) - a_third;
61
62 return true;
63 }
64
65 long double const A =
66 -sgn(R) * pow(fabs(R) + sqrt(R_squared - Q_cubed), 1.0/3.0L);
67
68 if (A != 0.0) { // in which case, B from NR is zero
69 *r1 = A + Q / A - a_third;
70 return false;
71 }
72
73 *r1 = *r2 = *r3 = -a_third;
74 return true;
75}
76
77// Returns the greatest common divisor of two unsigned integers x and y,
78// and assigns a, and b such that a*x + b*y = gcd(x, y).
79unsigned int MathUtil::ExtendedGCD(unsigned int x, unsigned int y,
80 int* a, int* b) {
81 *a = 1;
82 *b = 0;
83 int c = 0;
84 int d = 1;
85 // before and after each loop:
86 // current_x == a * original_x + b * original_y
87 // current_y == c * original_x + d * original_y
88 while (y != 0) {
89 // div() takes int parameters; there is no version that takes unsigned int
90 div_t r = div(static_cast<int>(x), static_cast<int>(y));
91 x = y;
92 y = r.rem;
93
94 int tmp = c;
95 c = *a - r.quot * c;
96 *a = tmp;
97
98 tmp = d;
99 d = *b - r.quot * d;
100 *b = tmp;
101 }
102 return x;
103}
104
105
106void MathUtil::ShardsToRead(const vector<bool>& shards_to_write,
107 vector<bool>* shards_to_read) {
108 const int N = shards_to_read->size();
109 const int M = shards_to_write.size();
110 CHECK(N > 0 || M == 0) << ": have shards to write but not to read";
111
112 // Input shard n of N can contribute to output shard m of M if there
113 // exists a record with sharding hash x s.t. n = x % N and m = x % M.
114 // Equivalently, there must exist s and t s.t. x = tN + n = sM + m,
115 // i.e., tN - sM = m - n. Since G = gcd(N, M) evenly divides tN - sM,
116 // G must also evenly divide m - n. Proof in the other direction is
117 // left as an exercise.
118 // Given output shard m, we should, therefore, read input shards n
119 // that satisfy (n - m) = kG, i.e., n = m + kG. Let 0 <= n < N.
120 // Then, 0 <= m + kG < N and, finally, -m / G <= k < (N - m) / G.
121
122 const int G = GCD(N, M);
123 shards_to_read->assign(N, false);
124 for (int m = 0; m < M; m++) {
125 if (!shards_to_write[m]) continue;
126 const int k_min = -m / G;
127 const int k_max = k_min + N / G;
128 for (int k = k_min; k < k_max; k++) {
129 (*shards_to_read)[m + k * G] = true;
130 }
131 }
132}
133
134double MathUtil::Harmonic(int64 const n, double *const e) {
135 CHECK_GT(n, 0);
136
137 // Hn ~ ln(n) + 0.5772156649 +
138 // + 1/(2n) - 1/(12n^2) + 1/(120n^4) - error,
139 // with 0 < error < 1/(256*n^4).
140
141 double const
142 d = static_cast<double>(n),
143 d2 = d * d,
144 d4 = d2 * d2;
145
146 return (log(d) + 0.5772156649) // ln + Gamma constant
147 + 1 / (2 * d) - 1 / (12 * d2) + 1 / (120 * d4)
148 - (*e = 1 / (256 * d4));
149}
150
151// The formula is extracted from the following page
152// http://en.wikipedia.org/w/index.php?title=Stirling%27s_approximation
153double MathUtil::Stirling(double n) {
154 static const double kLog2Pi = log(2 * M_PI);
155 const double logN = log(n);
156 return (n * logN
157 - n
158 + 0.5 * (kLog2Pi + logN) // 0.5 * log(2 * M_PI * n)
159 + 1 / (12 * n)
160 - 1 / (360 * n * n * n));
161}
162
163double MathUtil::LogCombinations(int n, int k) {
164 CHECK_GE(n, k);
165 CHECK_GT(n, 0);
166 CHECK_GE(k, 0);
167
168 // use symmetry to pick the shorter calculation
169 if (k > n / 2) {
170 k = n - k;
171 }
172
173 // If we have more than 30 logarithms to calculate, we'll use
174 // Stirling's approximation for log(n!).
175 if (k > 15) {
176 return Stirling(n) - Stirling(k) - Stirling(n - k);
177 } else {
178 double result = 0;
179 for (int i = 1; i <= k; i++) {
180 result += log(n - k + i) - log(i);
181 }
182 return result;
183 }
184}
185