1 // Copyright 2016 Google Inc. 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 // Author: ericv@google.com (Eric Veach)
17
18 #include "s2/s2predicates.h"
19 #include "s2/s2predicates_internal.h"
20
21 #include <algorithm>
22 #include <cfloat>
23 #include <cmath>
24 #include "s2/base/commandlineflags.h"
25 #include "s2/base/stringprintf.h"
26 #include <gtest/gtest.h>
27 #include "s2/third_party/absl/base/casts.h"
28 #include "s2/s1angle.h"
29 #include "s2/s1chord_angle.h"
30 #include "s2/s2edge_distances.h"
31 #include "s2/s2pointutil.h"
32 #include "s2/s2testing.h"
33 #include "s2/util/math/exactfloat/exactfloat.h"
34 #include "s2/util/math/vector.h"
35
36 DEFINE_int32(consistency_iters, 5000,
37 "Number of iterations for precision consistency tests");
38
39 using std::min;
40 using std::numeric_limits;
41 using std::pow;
42 using std::vector;
43
44 namespace s2pred {
45
TEST(epsilon_for_digits,recursion)46 TEST(epsilon_for_digits, recursion) {
47 EXPECT_EQ(1.0, epsilon_for_digits(0));
48 EXPECT_EQ(std::ldexp(1.0, -24), epsilon_for_digits(24));
49 EXPECT_EQ(std::ldexp(1.0, -53), epsilon_for_digits(53));
50 EXPECT_EQ(std::ldexp(1.0, -64), epsilon_for_digits(64));
51 EXPECT_EQ(std::ldexp(1.0, -106), epsilon_for_digits(106));
52 EXPECT_EQ(std::ldexp(1.0, -113), epsilon_for_digits(113));
53 }
54
TEST(rounding_epsilon,vs_numeric_limits)55 TEST(rounding_epsilon, vs_numeric_limits) {
56 // Check that rounding_epsilon<T>() returns the expected value for "float"
57 // and "double". We explicitly do not test "long double" since if this type
58 // is implemented using double-double arithmetic then the numeric_limits
59 // epsilon() value is completely unrelated to the maximum rounding error.
60 EXPECT_EQ(0.5 * numeric_limits<float>::epsilon(),
61 rounding_epsilon<float>());
62 EXPECT_EQ(0.5 * numeric_limits<double>::epsilon(),
63 rounding_epsilon<double>());
64 }
65
TEST(Sign,CollinearPoints)66 TEST(Sign, CollinearPoints) {
67 // The following points happen to be *exactly collinear* along a line that it
68 // approximate tangent to the surface of the unit sphere. In fact, C is the
69 // exact midpoint of the line segment AB. All of these points are close
70 // enough to unit length to satisfy S2::IsUnitLength().
71 S2Point a(0.72571927877036835, 0.46058825605889098, 0.51106749730504852);
72 S2Point b(0.7257192746638208, 0.46058826573818168, 0.51106749441312738);
73 S2Point c(0.72571927671709457, 0.46058826089853633, 0.51106749585908795);
74 EXPECT_EQ(c - a, b - c);
75 EXPECT_NE(0, Sign(a, b, c));
76 EXPECT_EQ(Sign(a, b, c), Sign(b, c, a));
77 EXPECT_EQ(Sign(a, b, c), -Sign(c, b, a));
78
79 // The points "x1" and "x2" are exactly proportional, i.e. they both lie
80 // on a common line through the origin. Both points are considered to be
81 // normalized, and in fact they both satisfy (x == x.Normalize()).
82 // Therefore the triangle (x1, x2, -x1) consists of three distinct points
83 // that all lie on a common line through the origin.
84 S2Point x1(0.99999999999999989, 1.4901161193847655e-08, 0);
85 S2Point x2(1, 1.4901161193847656e-08, 0);
86 EXPECT_EQ(x1, x1.Normalize());
87 EXPECT_EQ(x2, x2.Normalize());
88 EXPECT_NE(0, Sign(x1, x2, -x1));
89 EXPECT_EQ(Sign(x1, x2, -x1), Sign(x2, -x1, x1));
90 EXPECT_EQ(Sign(x1, x2, -x1), -Sign(-x1, x2, x1));
91
92 // Here are two more points that are distinct, exactly proportional, and
93 // that satisfy (x == x.Normalize()).
94 S2Point x3 = S2Point(1, 1, 1).Normalize();
95 S2Point x4 = 0.99999999999999989 * x3;
96 EXPECT_EQ(x3, x3.Normalize());
97 EXPECT_EQ(x4, x4.Normalize());
98 EXPECT_NE(x3, x4);
99 EXPECT_NE(0, Sign(x3, x4, -x3));
100
101 // The following two points demonstrate that Normalize() is not idempotent,
102 // i.e. y0.Normalize() != y0.Normalize().Normalize(). Both points satisfy
103 // S2::IsNormalized(), though, and the two points are exactly proportional.
104 S2Point y0 = S2Point(1, 1, 0);
105 S2Point y1 = y0.Normalize();
106 S2Point y2 = y1.Normalize();
107 EXPECT_NE(y1, y2);
108 EXPECT_EQ(y2, y2.Normalize());
109 EXPECT_NE(0, Sign(y1, y2, -y1));
110 EXPECT_EQ(Sign(y1, y2, -y1), Sign(y2, -y1, y1));
111 EXPECT_EQ(Sign(y1, y2, -y1), -Sign(-y1, y2, y1));
112 }
113
114 // This test repeatedly constructs some number of points that are on or nearly
115 // on a given great circle. Then it chooses one of these points as the
116 // "origin" and sorts the other points in CCW order around it. Of course,
117 // since the origin is on the same great circle as the points being sorted,
118 // nearly all of these tests are degenerate. It then does various consistency
119 // checks to verify that the points are indeed sorted in CCW order.
120 //
121 // It is easier to think about what this test is doing if you imagine that the
122 // points are in general position rather than on a great circle.
123 class SignTest : public testing::Test {
124 protected:
125 // The following method is used to sort a collection of points in CCW order
126 // around a given origin. It returns true if A comes before B in the CCW
127 // ordering (starting at an arbitrary fixed direction).
128 class LessCCW {
129 public:
LessCCW(const S2Point & origin,const S2Point & start)130 LessCCW(const S2Point& origin, const S2Point& start)
131 : origin_(origin), start_(start) {
132 }
operator ()(const S2Point & a,const S2Point & b) const133 bool operator()(const S2Point& a, const S2Point& b) const {
134 // OrderedCCW() acts like "<=", so we need to invert the comparison.
135 return !s2pred::OrderedCCW(start_, b, a, origin_);
136 }
137 private:
138 const S2Point origin_;
139 const S2Point start_;
140 };
141
142 // Given a set of points with no duplicates, first remove "origin" from
143 // "points" (if it exists) and then sort the remaining points in CCW order
144 // around "origin" putting the result in "sorted".
SortCCW(const vector<S2Point> & points,const S2Point & origin,vector<S2Point> * sorted)145 static void SortCCW(const vector<S2Point>& points, const S2Point& origin,
146 vector<S2Point>* sorted) {
147 // Make a copy of the points with "origin" removed.
148 sorted->clear();
149 std::remove_copy(points.begin(), points.end(), back_inserter(*sorted),
150 origin);
151
152 // Sort the points CCW around the origin starting at (*sorted)[0].
153 LessCCW less(origin, (*sorted)[0]);
154 std::sort(sorted->begin(), sorted->end(), less);
155 }
156
157 // Given a set of points sorted circularly CCW around "origin", and the
158 // index "start" of a point A, count the number of CCW triangles OAB over
159 // all sorted points B not equal to A. Also check that the results of the
160 // CCW tests are consistent with the hypothesis that the points are sorted.
CountCCW(const vector<S2Point> & sorted,const S2Point & origin,int start)161 static int CountCCW(const vector<S2Point>& sorted, const S2Point& origin,
162 int start) {
163 int num_ccw = 0;
164 int last_sign = 1;
165 const int n = sorted.size();
166 for (int j = 1; j < n; ++j) {
167 int sign = Sign(origin, sorted[start], sorted[(start + j) % n]);
168 EXPECT_NE(0, sign);
169 if (sign > 0) ++num_ccw;
170
171 // Since the points are sorted around the origin, we expect to see a
172 // (possibly empty) sequence of CCW triangles followed by a (possibly
173 // empty) sequence of CW triangles.
174 EXPECT_FALSE(sign > 0 && last_sign < 0);
175 last_sign = sign;
176 }
177 return num_ccw;
178 }
179
180 // Test exhaustively whether the points in "sorted" are sorted circularly
181 // CCW around "origin".
TestCCW(const vector<S2Point> & sorted,const S2Point & origin)182 static void TestCCW(const vector<S2Point>& sorted,
183 const S2Point& origin) {
184 const int n = sorted.size();
185 int total_num_ccw = 0;
186 int last_num_ccw = CountCCW(sorted, origin, n - 1);
187 for (int start = 0; start < n; ++start) {
188 int num_ccw = CountCCW(sorted, origin, start);
189 // Each iteration we increase the start index by 1, therefore the number
190 // of CCW triangles should decrease by at most 1.
191 EXPECT_GE(num_ccw, last_num_ccw - 1);
192 total_num_ccw += num_ccw;
193 last_num_ccw = num_ccw;
194 }
195 // We have tested all triangles of the form OAB. Exactly half of these
196 // should be CCW.
197 EXPECT_EQ(n * (n-1) / 2, total_num_ccw);
198 }
199
AddNormalized(const S2Point & a,vector<S2Point> * points)200 static void AddNormalized(const S2Point& a, vector<S2Point>* points) {
201 points->push_back(a.Normalize());
202 }
203
204 // Add two points A1 and A2 that are slightly offset from A along the
205 // tangent toward B, and such that A, A1, and A2 are exactly collinear
206 // (i.e. even with infinite-precision arithmetic).
AddTangentPoints(const S2Point & a,const S2Point & b,vector<S2Point> * points)207 static void AddTangentPoints(const S2Point& a, const S2Point& b,
208 vector<S2Point>* points) {
209 Vector3_d dir = S2::RobustCrossProd(a, b).CrossProd(a).Normalize();
210 if (dir == S2Point(0, 0, 0)) return;
211 for (;;) {
212 S2Point delta = 1e-15 * S2Testing::rnd.RandDouble() * dir;
213 if ((a + delta) != a && (a + delta) - a == a - (a - delta) &&
214 S2::IsUnitLength(a + delta) && S2::IsUnitLength(a - delta)) {
215 points->push_back(a + delta);
216 points->push_back(a - delta);
217 return;
218 }
219 }
220 }
221
222 // Add zero or more (but usually one) point that is likely to trigger
223 // Sign() degeneracies among the given points.
AddDegeneracy(vector<S2Point> * points)224 static void AddDegeneracy(vector<S2Point>* points) {
225 S2Testing::Random* rnd = &S2Testing::rnd;
226 S2Point a = (*points)[rnd->Uniform(points->size())];
227 S2Point b = (*points)[rnd->Uniform(points->size())];
228 int coord = rnd->Uniform(3);
229 switch (rnd->Uniform(8)) {
230 case 0:
231 // Add a random point (not uniformly distributed) along the great
232 // circle AB.
233 AddNormalized(rnd->UniformDouble(-1, 1) * a +
234 rnd->UniformDouble(-1, 1) * b, points);
235 break;
236 case 1:
237 // Perturb one coordinate by the minimum amount possible.
238 a[coord] = nextafter(a[coord], rnd->OneIn(2) ? 2 : -2);
239 AddNormalized(a, points);
240 break;
241 case 2:
242 // Perturb one coordinate by up to 1e-15.
243 a[coord] += 1e-15 * rnd->UniformDouble(-1, 1);
244 AddNormalized(a, points);
245 break;
246 case 3:
247 // Scale a point just enough so that it is different while still being
248 // considered normalized.
249 a *= rnd->OneIn(2) ? (1 + 2e-16) : (1 - 1e-16);
250 if (S2::IsUnitLength(a)) points->push_back(a);
251 break;
252 case 4: {
253 // Add the intersection point of AB with X=0, Y=0, or Z=0.
254 S2Point dir(0, 0, 0);
255 dir[coord] = rnd->OneIn(2) ? 1 : -1;
256 Vector3_d norm = S2::RobustCrossProd(a, b).Normalize();
257 if (norm.Norm2() > 0) {
258 AddNormalized(S2::RobustCrossProd(dir, norm), points);
259 }
260 break;
261 }
262 case 5:
263 // Add two closely spaced points along the tangent at A to the great
264 // circle through AB.
265 AddTangentPoints(a, b, points);
266 break;
267 case 6:
268 // Add two closely spaced points along the tangent at A to the great
269 // circle through A and the X-axis.
270 AddTangentPoints(a, S2Point(1, 0, 0), points);
271 break;
272 case 7:
273 // Add the negative of a point.
274 points->push_back(-a);
275 break;
276 }
277 }
278
279 // Sort the points around the given origin, and then do some consistency
280 // checks to verify that they are actually sorted.
SortAndTest(const vector<S2Point> & points,const S2Point & origin)281 static void SortAndTest(const vector<S2Point>& points,
282 const S2Point& origin) {
283 vector<S2Point> sorted;
284 SortCCW(points, origin, &sorted);
285 TestCCW(sorted, origin);
286 }
287
288 // Construct approximately "n" points near the great circle through A and B,
289 // then sort them and test whether they are sorted.
TestGreatCircle(S2Point a,S2Point b,int n)290 static void TestGreatCircle(S2Point a, S2Point b, int n) {
291 a = a.Normalize();
292 b = b.Normalize();
293 vector<S2Point> points;
294 points.push_back(a);
295 points.push_back(b);
296 while (points.size() < n) {
297 AddDegeneracy(&points);
298 }
299 // Remove any (0, 0, 0) points that were accidentically created, then sort
300 // the points and remove duplicates.
301 points.erase(std::remove(points.begin(), points.end(), S2Point(0, 0, 0)),
302 points.end());
303 std::sort(points.begin(), points.end());
304 points.erase(std::unique(points.begin(), points.end()), points.end());
305 EXPECT_GE(points.size(), n / 2);
306
307 SortAndTest(points, a);
308 SortAndTest(points, b);
309 for (const S2Point& origin : points) {
310 SortAndTest(points, origin);
311 }
312 }
313 };
314
TEST_F(SignTest,StressTest)315 TEST_F(SignTest, StressTest) {
316 // The run time of this test is *cubic* in the parameter below.
317 static const int kNumPointsPerCircle = 20;
318
319 // This test is randomized, so it is beneficial to run it several times.
320 for (int iter = 0; iter < 3; ++iter) {
321 // The most difficult great circles are the ones in the X-Y, Y-Z, and X-Z
322 // planes, for two reasons. First, when one or more coordinates are close
323 // to zero then the perturbations can be much smaller, since floating
324 // point numbers are spaced much more closely together near zero. (This
325 // tests the handling of things like underflow.) The second reason is
326 // that most of the cases of SymbolicallyPerturbedSign() can only be
327 // reached when one or more input point coordinates are zero.
328 TestGreatCircle(S2Point(1, 0, 0), S2Point(0, 1, 0), kNumPointsPerCircle);
329 TestGreatCircle(S2Point(1, 0, 0), S2Point(0, 0, 1), kNumPointsPerCircle);
330 TestGreatCircle(S2Point(0, -1, 0), S2Point(0, 0, 1), kNumPointsPerCircle);
331
332 // This tests a great circle where at least some points have X, Y, and Z
333 // coordinates with exactly the same mantissa. One useful property of
334 // such points is that when they are scaled (e.g. multiplying by 1+eps),
335 // all such points are exactly collinear with the origin.
336 TestGreatCircle(S2Point(1 << 25, 1, -8), S2Point(-4, -(1 << 20), 1),
337 kNumPointsPerCircle);
338 }
339 }
340
341 class StableSignTest : public testing::Test {
342 protected:
343 // Estimate the probability that S2::StableSign() will not be able to compute
344 // the determinant sign of a triangle A, B, C consisting of three points
345 // that are as collinear as possible and spaced the given distance apart.
GetFailureRate(double km)346 double GetFailureRate(double km) {
347 const int kIters = 1000;
348 int failure_count = 0;
349 double m = tan(S2Testing::KmToAngle(km).radians());
350 for (int iter = 0; iter < kIters; ++iter) {
351 S2Point a, x, y;
352 S2Testing::GetRandomFrame(&a, &x, &y);
353 S2Point b = (a - m * x).Normalize();
354 S2Point c = (a + m * x).Normalize();
355 int sign = s2pred::StableSign(a, b, c);
356 if (sign != 0) {
357 EXPECT_EQ(s2pred::ExactSign(a, b, c, true), sign);
358 } else {
359 ++failure_count;
360 }
361 }
362 double rate = static_cast<double>(failure_count) / kIters;
363 S2_LOG(INFO) << "StableSign failure rate for " << km << " km = " << rate;
364 return rate;
365 }
366 };
367
TEST_F(StableSignTest,FailureRate)368 TEST_F(StableSignTest, FailureRate) {
369 // Verify that StableSign() is able to handle most cases where the three
370 // points are as collinear as possible. (For reference, TriageSign() fails
371 // virtually 100% of the time on this test.)
372 //
373 // Note that the failure rate *decreases* as the points get closer together,
374 // and the decrease is approximately linear. For example, the failure rate
375 // is 0.4% for collinear points spaced 1km apart, but only 0.0004% for
376 // collinear points spaced 1 meter apart.
377
378 EXPECT_LT(GetFailureRate(1.0), 0.01); // 1km spacing: < 1% (actual 0.4%)
379 EXPECT_LT(GetFailureRate(10.0), 0.1); // 10km spacing: < 10% (actual 4%)
380 }
381
382 // Given 3 points A, B, C that are exactly coplanar with the origin and where
383 // A < B < C in lexicographic order, verify that ABC is counterclockwise (if
384 // expected == 1) or clockwise (if expected == -1) using ExpensiveSign().
385 //
386 // This method is intended specifically for checking the cases where
387 // symbolic perturbations are needed to break ties.
CheckSymbolicSign(int expected,const S2Point & a,const S2Point & b,const S2Point & c)388 static void CheckSymbolicSign(int expected, const S2Point& a,
389 const S2Point& b, const S2Point& c) {
390 S2_CHECK_LT(a, b);
391 S2_CHECK_LT(b, c);
392 S2_CHECK_EQ(0, a.DotProd(b.CrossProd(c)));
393
394 // Use ASSERT rather than EXPECT to suppress spurious error messages.
395 ASSERT_EQ(expected, ExpensiveSign(a, b, c));
396 ASSERT_EQ(expected, ExpensiveSign(b, c, a));
397 ASSERT_EQ(expected, ExpensiveSign(c, a, b));
398 ASSERT_EQ(-expected, ExpensiveSign(c, b, a));
399 ASSERT_EQ(-expected, ExpensiveSign(b, a, c));
400 ASSERT_EQ(-expected, ExpensiveSign(a, c, b));
401 }
402
TEST(Sign,SymbolicPerturbationCodeCoverage)403 TEST(Sign, SymbolicPerturbationCodeCoverage) {
404 // The purpose of this test is simply to get code coverage of
405 // SymbolicallyPerturbedSign(). Let M_1, M_2, ... be the sequence of
406 // submatrices whose determinant sign is tested by that function. Then the
407 // i-th test below is a 3x3 matrix M (with rows A, B, C) such that:
408 //
409 // det(M) = 0
410 // det(M_j) = 0 for j < i
411 // det(M_i) != 0
412 // A < B < C in lexicographic order.
413 //
414 // I checked that reversing the sign of any of the "return" statements in
415 // SymbolicallyPerturbedSign() will cause this test to fail.
416
417 // det(M_1) = b0*c1 - b1*c0
418 CheckSymbolicSign(1,
419 S2Point(-3, -1, 0), S2Point(-2, 1, 0), S2Point(1, -2, 0));
420
421 // det(M_2) = b2*c0 - b0*c2
422 CheckSymbolicSign(1,
423 S2Point(-6, 3, 3), S2Point(-4, 2, -1), S2Point(-2, 1, 4));
424
425 // det(M_3) = b1*c2 - b2*c1
426 CheckSymbolicSign(1, S2Point(0, -1, -1), S2Point(0, 1, -2), S2Point(0, 2, 1));
427 // From this point onward, B or C must be zero, or B is proportional to C.
428
429 // det(M_4) = c0*a1 - c1*a0
430 CheckSymbolicSign(1, S2Point(-1, 2, 7), S2Point(2, 1, -4), S2Point(4, 2, -8));
431
432 // det(M_5) = c0
433 CheckSymbolicSign(1,
434 S2Point(-4, -2, 7), S2Point(2, 1, -4), S2Point(4, 2, -8));
435
436 // det(M_6) = -c1
437 CheckSymbolicSign(1, S2Point(0, -5, 7), S2Point(0, -4, 8), S2Point(0, -2, 4));
438
439 // det(M_7) = c2*a0 - c0*a2
440 CheckSymbolicSign(1,
441 S2Point(-5, -2, 7), S2Point(0, 0, -2), S2Point(0, 0, -1));
442
443 // det(M_8) = c2
444 CheckSymbolicSign(1, S2Point(0, -2, 7), S2Point(0, 0, 1), S2Point(0, 0, 2));
445 // From this point onward, C must be zero.
446
447 // det(M_9) = a0*b1 - a1*b0
448 CheckSymbolicSign(1, S2Point(-3, 1, 7), S2Point(-1, -4, 1), S2Point(0, 0, 0));
449
450 // det(M_10) = -b0
451 CheckSymbolicSign(1,
452 S2Point(-6, -4, 7), S2Point(-3, -2, 1), S2Point(0, 0, 0));
453
454 // det(M_11) = b1
455 CheckSymbolicSign(-1, S2Point(0, -4, 7), S2Point(0, -2, 1), S2Point(0, 0, 0));
456
457 // det(M_12) = a0
458 CheckSymbolicSign(-1,
459 S2Point(-1, -4, 5), S2Point(0, 0, -3), S2Point(0, 0, 0));
460
461 // det(M_13) = 1
462 CheckSymbolicSign(1, S2Point(0, -4, 5), S2Point(0, 0, -5), S2Point(0, 0, 0));
463 }
464
465 enum Precision { DOUBLE, LONG_DOUBLE, EXACT, SYMBOLIC, NUM_PRECISIONS };
466
467 static const char* kPrecisionNames[] = {
468 "double", "long double", "exact", "symbolic"
469 };
470
471 // A helper class that keeps track of how often each precision was used and
472 // generates a string for logging purposes.
473 class PrecisionStats {
474 public:
475 PrecisionStats();
Tally(Precision precision)476 void Tally(Precision precision) { ++counts_[precision]; }
477 string ToString();
478
479 private:
480 int counts_[NUM_PRECISIONS];
481 };
482
PrecisionStats()483 PrecisionStats::PrecisionStats() {
484 for (int& count : counts_) count = 0;
485 }
486
ToString()487 string PrecisionStats::ToString() {
488 string result;
489 int total = 0;
490 for (int i = 0; i < NUM_PRECISIONS; ++i) {
491 StringAppendF(&result, "%s=%6d, ", kPrecisionNames[i], counts_[i]);
492 total += counts_[i];
493 }
494 StringAppendF(&result, "total=%6d", total);
495 return result;
496 }
497
498 // Chooses a random S2Point that is often near the intersection of one of the
499 // coodinates planes or coordinate axes with the unit sphere. (It is possible
500 // to represent very small perturbations near such points.)
ChoosePoint()501 static S2Point ChoosePoint() {
502 S2Point x = S2Testing::RandomPoint();
503 for (int i = 0; i < 3; ++i) {
504 if (S2Testing::rnd.OneIn(3)) {
505 x[i] *= pow(1e-50, S2Testing::rnd.RandDouble());
506 }
507 }
508 return x.Normalize();
509 }
510
511 // The following helper classes allow us to test the various distance
512 // calculation methods using a common test framework.
513 class Sin2Distances {
514 public:
515 template <class T>
Triage(const Vector3<T> & x,const Vector3<T> & a,const Vector3<T> & b)516 static int Triage(const Vector3<T>& x,
517 const Vector3<T>& a, const Vector3<T>& b) {
518 return TriageCompareSin2Distances(x, a, b);
519 }
520 };
521
522 class CosDistances {
523 public:
524 template <class T>
Triage(const Vector3<T> & x,const Vector3<T> & a,const Vector3<T> & b)525 static int Triage(const Vector3<T>& x,
526 const Vector3<T>& a, const Vector3<T>& b) {
527 return TriageCompareCosDistances(x, a, b);
528 }
529 };
530
531 // Compares distances greater than 90 degrees using sin^2(distance).
532 class MinusSin2Distances {
533 public:
534 template <class T>
Triage(const Vector3<T> & x,const Vector3<T> & a,const Vector3<T> & b)535 static int Triage(const Vector3<T>& x,
536 const Vector3<T>& a, const Vector3<T>& b) {
537 return -TriageCompareSin2Distances(-x, a, b);
538 }
539 };
540
541 // Verifies that CompareDistances(x, a, b) == expected_sign, and furthermore
542 // checks that the minimum required precision is "expected_prec" when the
543 // distance calculation method defined by CompareDistancesWrapper is used.
544 template <class CompareDistancesWrapper>
TestCompareDistances(S2Point x,S2Point a,S2Point b,int expected_sign,Precision expected_prec)545 void TestCompareDistances(S2Point x, S2Point a, S2Point b,
546 int expected_sign, Precision expected_prec) {
547 // Don't normalize the arguments unless necessary (to allow testing points
548 // that differ only in magnitude).
549 if (!S2::IsUnitLength(x)) x = x.Normalize();
550 if (!S2::IsUnitLength(a)) a = a.Normalize();
551 if (!S2::IsUnitLength(b)) b = b.Normalize();
552
553 int dbl_sign = CompareDistancesWrapper::Triage(x, a, b);
554 int ld_sign = CompareDistancesWrapper::Triage(ToLD(x), ToLD(a), ToLD(b));
555 int exact_sign = ExactCompareDistances(ToExact(x), ToExact(a), ToExact(b));
556 int actual_sign = (exact_sign != 0 ? exact_sign :
557 SymbolicCompareDistances(x, a, b));
558
559 // Check that the signs are correct (if non-zero), and also that if dbl_sign
560 // is non-zero then so is ld_sign, etc.
561 EXPECT_EQ(expected_sign, actual_sign);
562 if (exact_sign != 0) EXPECT_EQ(exact_sign, actual_sign);
563 if (ld_sign != 0) EXPECT_EQ(exact_sign, ld_sign);
564 if (dbl_sign != 0) EXPECT_EQ(ld_sign, dbl_sign);
565
566 Precision actual_prec = (dbl_sign ? DOUBLE :
567 ld_sign ? LONG_DOUBLE :
568 exact_sign ? EXACT : SYMBOLIC);
569 EXPECT_EQ(expected_prec, actual_prec);
570
571 // Make sure that the top-level function returns the expected result.
572 EXPECT_EQ(expected_sign, CompareDistances(x, a, b));
573
574 // Check that reversing the arguments negates the result.
575 EXPECT_EQ(-expected_sign, CompareDistances(x, b, a));
576 }
577
TEST(CompareDistances,Coverage)578 TEST(CompareDistances, Coverage) {
579 // This test attempts to exercise all the code paths in all precisions.
580
581 // Test TriageCompareSin2Distances.
582 TestCompareDistances<Sin2Distances>(
583 S2Point(1, 1, 1), S2Point(1, 1 - 1e-15, 1), S2Point(1, 1, 1 + 2e-15),
584 -1, DOUBLE);
585 TestCompareDistances<Sin2Distances>(
586 S2Point(1, 1, 0), S2Point(1, 1 - 1e-15, 1e-21), S2Point(1, 1 - 1e-15, 0),
587 1, DOUBLE);
588 TestCompareDistances<Sin2Distances>(
589 S2Point(2, 0, 0), S2Point(2, -1, 0), S2Point(2, 1, 1e-8),
590 -1, LONG_DOUBLE);
591 TestCompareDistances<Sin2Distances>(
592 S2Point(2, 0, 0), S2Point(2, -1, 0), S2Point(2, 1, 1e-100),
593 -1, EXACT);
594 TestCompareDistances<Sin2Distances>(
595 S2Point(1, 0, 0), S2Point(1, -1, 0), S2Point(1, 1, 0),
596 1, SYMBOLIC);
597 TestCompareDistances<Sin2Distances>(
598 S2Point(1, 0, 0), S2Point(1, 0, 0), S2Point(1, 0, 0),
599 0, SYMBOLIC);
600
601 // Test TriageCompareCosDistances.
602 TestCompareDistances<CosDistances>(
603 S2Point(1, 1, 1), S2Point(1, -1, 0), S2Point(-1, 1, 3e-15),
604 1, DOUBLE);
605 TestCompareDistances<CosDistances>(
606 S2Point(1, 0, 0), S2Point(1, 1e-30, 0), S2Point(-1, 1e-40, 0),
607 -1, DOUBLE);
608 TestCompareDistances<CosDistances>(
609 S2Point(1, 1, 1), S2Point(1, -1, 0), S2Point(-1, 1, 3e-18),
610 1, LONG_DOUBLE);
611 TestCompareDistances<CosDistances>(
612 S2Point(1, 1, 1), S2Point(1, -1, 0), S2Point(-1, 1, 1e-100),
613 1, EXACT);
614 TestCompareDistances<CosDistances>(
615 S2Point(1, 1, 1), S2Point(1, -1, 0), S2Point(-1, 1, 0),
616 -1, SYMBOLIC);
617 TestCompareDistances<CosDistances>(
618 S2Point(1, 1, 1), S2Point(1, -1, 0), S2Point(1, -1, 0),
619 0, SYMBOLIC);
620
621 // Test TriageCompareSin2Distances using distances greater than 90 degrees.
622 TestCompareDistances<MinusSin2Distances>(
623 S2Point(1, 1, 0), S2Point(-1, -1 + 1e-15, 0), S2Point(-1, -1, 0),
624 -1, DOUBLE);
625 TestCompareDistances<MinusSin2Distances>(
626 S2Point(-1, -1, 0), S2Point(1, 1 - 1e-15, 0),
627 S2Point(1, 1 - 1e-15, 1e-21), 1, DOUBLE);
628 TestCompareDistances<MinusSin2Distances>(
629 S2Point(-1, -1, 0), S2Point(2, 1, 0), S2Point(2, 1, 1e-8),
630 1, LONG_DOUBLE);
631 TestCompareDistances<MinusSin2Distances>(
632 S2Point(-1, -1, 0), S2Point(2, 1, 0), S2Point(2, 1, 1e-30),
633 1, EXACT);
634 TestCompareDistances<MinusSin2Distances>(
635 S2Point(-1, -1, 0), S2Point(2, 1, 0), S2Point(1, 2, 0),
636 -1, SYMBOLIC);
637 }
638
639 // Checks that the result at one level of precision is consistent with the
640 // result at the next higher level of precision. Returns the minimum
641 // precision that yielded a non-zero result.
642 template <class CompareDistancesWrapper>
TestCompareDistancesConsistency(const S2Point & x,const S2Point & a,const S2Point & b)643 Precision TestCompareDistancesConsistency(const S2Point& x, const S2Point& a,
644 const S2Point& b) {
645 int dbl_sign = CompareDistancesWrapper::Triage(x, a, b);
646 int ld_sign = CompareDistancesWrapper::Triage(ToLD(x), ToLD(a), ToLD(b));
647 int exact_sign = ExactCompareDistances(ToExact(x), ToExact(a), ToExact(b));
648 if (dbl_sign != 0) EXPECT_EQ(ld_sign, dbl_sign);
649 if (ld_sign != 0) EXPECT_EQ(exact_sign, ld_sign);
650 if (exact_sign != 0) {
651 EXPECT_EQ(exact_sign, CompareDistances(x, a, b));
652 return (ld_sign == 0) ? EXACT : (dbl_sign == 0) ? LONG_DOUBLE : DOUBLE;
653 } else {
654 // Unlike the other methods, SymbolicCompareDistances has the
655 // precondition that the exact sign must be zero.
656 int symbolic_sign = SymbolicCompareDistances(x, a, b);
657 EXPECT_EQ(symbolic_sign, CompareDistances(x, a, b));
658 return SYMBOLIC;
659 }
660 }
661
TEST(CompareDistances,Consistency)662 TEST(CompareDistances, Consistency) {
663 // This test chooses random point pairs that are nearly equidistant from a
664 // target point, and then checks that the answer given by a method at one
665 // level of precision is consistent with the answer given at the next higher
666 // level of precision.
667 //
668 // The way the .cc file is structured, we can only do comparisons using a
669 // specific precision if we also choose the specific distance calculation
670 // method. The code below checks that the Cos, Sin2, and MinusSin2 methods
671 // are consistent across their entire valid range of inputs, and also
672 // simulates the logic in CompareDistance that chooses which method to use
673 // in order to gather statistics about how often each precision is needed.
674 // (These statistics are only useful for coverage purposes, not benchmarks,
675 // since the input points are chosen to be pathological worst cases.)
676 TestCompareDistancesConsistency<CosDistances>(
677 S2Point(1, 0, 0), S2Point(0, -1, 0), S2Point(0, 1, 0));
678 auto& rnd = S2Testing::rnd;
679 PrecisionStats sin2_stats, cos_stats, minus_sin2_stats;
680 for (int iter = 0; iter < FLAGS_consistency_iters; ++iter) {
681 rnd.Reset(iter + 1); // Easier to reproduce a specific case.
682 S2Point x = ChoosePoint();
683 S2Point dir = ChoosePoint();
684 S1Angle r = S1Angle::Radians(M_PI_2 * pow(1e-30, rnd.RandDouble()));
685 if (rnd.OneIn(2)) r = S1Angle::Radians(M_PI_2) - r;
686 if (rnd.OneIn(2)) r = S1Angle::Radians(M_PI_2) + r;
687 S2Point a = S2::InterpolateAtDistance(r, x, dir);
688 S2Point b = S2::InterpolateAtDistance(r, x, -dir);
689 Precision prec = TestCompareDistancesConsistency<CosDistances>(x, a, b);
690 if (r.degrees() >= 45 && r.degrees() <= 135) cos_stats.Tally(prec);
691 // The Sin2 method is only valid if both distances are less than 90
692 // degrees, and similarly for the MinusSin2 method. (In the actual
693 // implementation these methods are only used if both distances are less
694 // than 45 degrees or greater than 135 degrees respectively.)
695 if (r.radians() < M_PI_2 - 1e-14) {
696 prec = TestCompareDistancesConsistency<Sin2Distances>(x, a, b);
697 if (r.degrees() < 45) {
698 // Don't skew the statistics by recording degenerate inputs.
699 if (a == b) {
700 EXPECT_EQ(SYMBOLIC, prec);
701 } else {
702 sin2_stats.Tally(prec);
703 }
704 }
705 } else if (r.radians() > M_PI_2 + 1e-14) {
706 prec = TestCompareDistancesConsistency<MinusSin2Distances>(x, a, b);
707 if (r.degrees() > 135) minus_sin2_stats.Tally(prec);
708 }
709 }
710 S2_LOG(ERROR) << "\nsin2: " << sin2_stats.ToString()
711 << "\ncos: " << cos_stats.ToString()
712 << "\n-sin2: " << minus_sin2_stats.ToString();
713 }
714
715 // Helper classes for testing the various distance calculation methods.
716 class Sin2Distance {
717 public:
718 template <class T>
Triage(const Vector3<T> & x,const Vector3<T> & y,S1ChordAngle r)719 static int Triage(const Vector3<T>& x, const Vector3<T>& y, S1ChordAngle r) {
720 return TriageCompareSin2Distance(x, y, absl::implicit_cast<T>(r.length2()));
721 }
722 };
723
724 class CosDistance {
725 public:
726 template <class T>
Triage(const Vector3<T> & x,const Vector3<T> & y,S1ChordAngle r)727 static int Triage(const Vector3<T>& x, const Vector3<T>& y, S1ChordAngle r) {
728 return TriageCompareCosDistance(x, y, absl::implicit_cast<T>(r.length2()));
729 }
730 };
731
732 // Verifies that CompareDistance(x, y, r) == expected_sign, and furthermore
733 // checks that the minimum required precision is "expected_prec" when the
734 // distance calculation method defined by CompareDistanceWrapper is used.
735 template <class CompareDistanceWrapper>
TestCompareDistance(S2Point x,S2Point y,S1ChordAngle r,int expected_sign,Precision expected_prec)736 void TestCompareDistance(S2Point x, S2Point y, S1ChordAngle r,
737 int expected_sign, Precision expected_prec) {
738 // Don't normalize the arguments unless necessary (to allow testing points
739 // that differ only in magnitude).
740 if (!S2::IsUnitLength(x)) x = x.Normalize();
741 if (!S2::IsUnitLength(y)) y = y.Normalize();
742
743 int dbl_sign = CompareDistanceWrapper::Triage(x, y, r);
744 int ld_sign = CompareDistanceWrapper::Triage(ToLD(x), ToLD(y), r);
745 int exact_sign = ExactCompareDistance(ToExact(x), ToExact(y), r.length2());
746
747 // Check that the signs are correct (if non-zero), and also that if dbl_sign
748 // is non-zero then so is ld_sign, etc.
749 EXPECT_EQ(expected_sign, exact_sign);
750 if (ld_sign != 0) EXPECT_EQ(exact_sign, ld_sign);
751 if (dbl_sign != 0) EXPECT_EQ(ld_sign, dbl_sign);
752
753 Precision actual_prec = (dbl_sign ? DOUBLE : ld_sign ? LONG_DOUBLE : EXACT);
754 EXPECT_EQ(expected_prec, actual_prec);
755
756 // Make sure that the top-level function returns the expected result.
757 EXPECT_EQ(expected_sign, CompareDistance(x, y, r));
758
759 // Mathematically, if d(X, Y) < r then d(-X, Y) > (Pi - r). Unfortunately
760 // there can be rounding errors when computing the supplementary distance,
761 // so to ensure the two distances are exactly supplementary we need to do
762 // the following.
763 S1ChordAngle r_supp = S1ChordAngle::Straight() - r;
764 r = S1ChordAngle::Straight() - r_supp;
765 EXPECT_EQ(-CompareDistance(x, y, r), CompareDistance(-x, y, r_supp));
766 }
767
TEST(CompareDistance,Coverage)768 TEST(CompareDistance, Coverage) {
769 // Test TriageCompareSin2Distance.
770 TestCompareDistance<Sin2Distance>(
771 S2Point(1, 1, 1), S2Point(1, 1 - 1e-15, 1),
772 S1ChordAngle::Radians(1e-15), -1, DOUBLE);
773 TestCompareDistance<Sin2Distance>(
774 S2Point(1, 0, 0), S2Point(1, 1, 0),
775 S1ChordAngle::Radians(M_PI_4), -1, LONG_DOUBLE);
776 TestCompareDistance<Sin2Distance>(
777 S2Point(1, 1e-40, 0), S2Point(1 + DBL_EPSILON, 1e-40, 0),
778 S1ChordAngle::Radians(0.9 * DBL_EPSILON * 1e-40), 1, EXACT);
779 TestCompareDistance<Sin2Distance>(
780 S2Point(1, 1e-40, 0), S2Point(1 + DBL_EPSILON, 1e-40, 0),
781 S1ChordAngle::Radians(1.1 * DBL_EPSILON * 1e-40), -1, EXACT);
782 TestCompareDistance<Sin2Distance>(
783 S2Point(1, 0, 0), S2Point(1 + DBL_EPSILON, 0, 0),
784 S1ChordAngle::Zero(), 0, EXACT);
785
786 // Test TriageCompareCosDistance.
787 TestCompareDistance<CosDistance>(
788 S2Point(1, 0, 0), S2Point(1, 1e-8, 0),
789 S1ChordAngle::Radians(1e-7), -1, DOUBLE);
790 TestCompareDistance<CosDistance>(
791 S2Point(1, 0, 0), S2Point(-1, 1e-8, 0),
792 S1ChordAngle::Radians(M_PI - 1e-7), 1, DOUBLE);
793 TestCompareDistance<CosDistance>(
794 S2Point(1, 1, 0), S2Point(1, -1 - 2 * DBL_EPSILON, 0),
795 S1ChordAngle::Right(), 1, DOUBLE);
796 TestCompareDistance<CosDistance>(
797 S2Point(1, 1, 0), S2Point(1, -1 - DBL_EPSILON, 0),
798 S1ChordAngle::Right(), 1, LONG_DOUBLE);
799 TestCompareDistance<CosDistance>(
800 S2Point(1, 1, 0), S2Point(1, -1, 1e-30),
801 S1ChordAngle::Right(), 0, EXACT);
802 // The angle between these two points is exactly 60 degrees.
803 TestCompareDistance<CosDistance>(
804 S2Point(1, 1, 0), S2Point(0, 1, 1),
805 S1ChordAngle::FromLength2(1), 0, EXACT);
806 }
807
808 // Checks that the result at one level of precision is consistent with the
809 // result at the next higher level of precision. Returns the minimum
810 // precision that yielded a non-zero result.
811 template <class CompareDistanceWrapper>
TestCompareDistanceConsistency(const S2Point & x,const S2Point & y,S1ChordAngle r)812 Precision TestCompareDistanceConsistency(const S2Point& x, const S2Point& y,
813 S1ChordAngle r) {
814 int dbl_sign = CompareDistanceWrapper::Triage(x, y, r);
815 int ld_sign = CompareDistanceWrapper::Triage(ToLD(x), ToLD(y), r);
816 int exact_sign = ExactCompareDistance(ToExact(x), ToExact(y), r.length2());
817 EXPECT_EQ(exact_sign, CompareDistance(x, y, r));
818 if (dbl_sign != 0) EXPECT_EQ(ld_sign, dbl_sign);
819 if (ld_sign != 0) EXPECT_EQ(exact_sign, ld_sign);
820 return (ld_sign == 0) ? EXACT : (dbl_sign == 0) ? LONG_DOUBLE : DOUBLE;
821 }
822
TEST(CompareDistance,Consistency)823 TEST(CompareDistance, Consistency) {
824 // This test chooses random inputs such that the distance between points X
825 // and Y is very close to the threshold distance "r". It then checks that
826 // the answer given by a method at one level of precision is consistent with
827 // the answer given at the next higher level of precision. See also the
828 // comments in the CompareDistances consistency test.
829 auto& rnd = S2Testing::rnd;
830 PrecisionStats sin2_stats, cos_stats;
831 for (int iter = 0; iter < FLAGS_consistency_iters; ++iter) {
832 rnd.Reset(iter + 1); // Easier to reproduce a specific case.
833 S2Point x = ChoosePoint();
834 S2Point dir = ChoosePoint();
835 S1Angle r = S1Angle::Radians(M_PI_2 * pow(1e-30, rnd.RandDouble()));
836 if (rnd.OneIn(2)) r = S1Angle::Radians(M_PI_2) - r;
837 if (rnd.OneIn(5)) r = S1Angle::Radians(M_PI_2) + r;
838 S2Point y = S2::InterpolateAtDistance(r, x, dir);
839 Precision prec = TestCompareDistanceConsistency<CosDistance>(
840 x, y, S1ChordAngle(r));
841 if (r.degrees() >= 45) cos_stats.Tally(prec);
842 if (r.radians() < M_PI_2 - 1e-14) {
843 prec = TestCompareDistanceConsistency<Sin2Distance>(
844 x, y, S1ChordAngle(r));
845 if (r.degrees() < 45) sin2_stats.Tally(prec);
846 }
847 }
848 S2_LOG(ERROR) << "\nsin2: " << sin2_stats.ToString()
849 << "\ncos: " << cos_stats.ToString();
850 }
851
852 // Verifies that CompareEdgeDistance(x, a0, a1, r) == expected_sign, and
853 // furthermore checks that the minimum required precision is "expected_prec".
TestCompareEdgeDistance(S2Point x,S2Point a0,S2Point a1,S1ChordAngle r,int expected_sign,Precision expected_prec)854 void TestCompareEdgeDistance(S2Point x, S2Point a0, S2Point a1, S1ChordAngle r,
855 int expected_sign, Precision expected_prec) {
856 // Don't normalize the arguments unless necessary (to allow testing points
857 // that differ only in magnitude).
858 if (!S2::IsUnitLength(x)) x = x.Normalize();
859 if (!S2::IsUnitLength(a0)) a0 = a0.Normalize();
860 if (!S2::IsUnitLength(a1)) a1 = a1.Normalize();
861
862 int dbl_sign = TriageCompareEdgeDistance(x, a0, a1, r.length2());
863 int ld_sign = TriageCompareEdgeDistance(ToLD(x), ToLD(a0), ToLD(a1),
864 ToLD(r.length2()));
865 int exact_sign = ExactCompareEdgeDistance(x, a0, a1, r);
866
867 // Check that the signs are correct (if non-zero), and also that if dbl_sign
868 // is non-zero then so is ld_sign, etc.
869 EXPECT_EQ(expected_sign, exact_sign);
870 if (ld_sign != 0) EXPECT_EQ(exact_sign, ld_sign);
871 if (dbl_sign != 0) EXPECT_EQ(ld_sign, dbl_sign);
872
873 Precision actual_prec = (dbl_sign ? DOUBLE : ld_sign ? LONG_DOUBLE : EXACT);
874 EXPECT_EQ(expected_prec, actual_prec);
875
876 // Make sure that the top-level function returns the expected result.
877 EXPECT_EQ(expected_sign, CompareEdgeDistance(x, a0, a1, r));
878 }
879
TEST(CompareEdgeDistance,Coverage)880 TEST(CompareEdgeDistance, Coverage) {
881 // Test TriageCompareLineSin2Distance.
882 TestCompareEdgeDistance(
883 S2Point(1, 1e-10, 1e-15), S2Point(1, 0, 0), S2Point(0, 1, 0),
884 S1ChordAngle::Radians(1e-15 + DBL_EPSILON), -1, DOUBLE);
885 TestCompareEdgeDistance(
886 S2Point(1, 1, 1e-15), S2Point(1, 0, 0), S2Point(0, 1, 0),
887 S1ChordAngle::Radians(1e-15 + DBL_EPSILON), -1, LONG_DOUBLE);
888 TestCompareEdgeDistance(
889 S2Point(1, 1, 1e-40), S2Point(1, 0, 0), S2Point(0, 1, 0),
890 S1ChordAngle::Radians(1e-40), -1, EXACT);
891 TestCompareEdgeDistance(
892 S2Point(1, 1, 0), S2Point(1, 0, 0), S2Point(0, 1, 0),
893 S1ChordAngle::Zero(), 0, EXACT);
894
895 // Test TriageCompareLineCos2Distance.
896 TestCompareEdgeDistance(
897 S2Point(1e-15, 0, 1), S2Point(1, 0, 0), S2Point(0, 1, 0),
898 S1ChordAngle::Radians(M_PI_2 - 1e-15 - 3 * DBL_EPSILON),
899 1, DOUBLE);
900 TestCompareEdgeDistance(
901 S2Point(1e-15, 0, 1), S2Point(1, 0, 0), S2Point(0, 1, 0),
902 S1ChordAngle::Radians(M_PI_2 - 1e-15 - DBL_EPSILON),
903 1, LONG_DOUBLE);
904 TestCompareEdgeDistance(
905 S2Point(1e-40, 0, 1), S2Point(1, 0, 0), S2Point(0, 1, 0),
906 S1ChordAngle::Right(), -1, EXACT);
907 TestCompareEdgeDistance(
908 S2Point(0, 0, 1), S2Point(1, 0, 0), S2Point(0, 1, 0),
909 S1ChordAngle::Right(), 0, EXACT);
910
911 // Test cases where the closest point is an edge endpoint.
912 TestCompareEdgeDistance(
913 S2Point(1e-15, -1, 0), S2Point(1, 0, 0), S2Point(1, 1, 0),
914 S1ChordAngle::Right(), -1, DOUBLE);
915 TestCompareEdgeDistance(
916 S2Point(1e-18, -1, 0), S2Point(1, 0, 0), S2Point(1, 1, 0),
917 S1ChordAngle::Right(), -1, LONG_DOUBLE);
918 TestCompareEdgeDistance(
919 S2Point(1e-100, -1, 0), S2Point(1, 0, 0), S2Point(1, 1, 0),
920 S1ChordAngle::Right(), -1, EXACT);
921 TestCompareEdgeDistance(
922 S2Point(0, -1, 0), S2Point(1, 0, 0), S2Point(1, 1, 0),
923 S1ChordAngle::Right(), 0, EXACT);
924 }
925
926 // Checks that the result at one level of precision is consistent with the
927 // result at the next higher level of precision. Returns the minimum
928 // precision that yielded a non-zero result.
TestCompareEdgeDistanceConsistency(const S2Point & x,const S2Point & a0,const S2Point & a1,S1ChordAngle r)929 Precision TestCompareEdgeDistanceConsistency(
930 const S2Point& x, const S2Point& a0, const S2Point& a1, S1ChordAngle r) {
931 int dbl_sign = TriageCompareEdgeDistance(x, a0, a1, r.length2());
932 int ld_sign = TriageCompareEdgeDistance(ToLD(x), ToLD(a0), ToLD(a1),
933 ToLD(r.length2()));
934 int exact_sign = ExactCompareEdgeDistance(x, a0, a1, r);
935 EXPECT_EQ(exact_sign, CompareEdgeDistance(x, a0, a1, r));
936 if (dbl_sign != 0) EXPECT_EQ(ld_sign, dbl_sign);
937 if (ld_sign != 0) EXPECT_EQ(exact_sign, ld_sign);
938 return (ld_sign == 0) ? EXACT : (dbl_sign == 0) ? LONG_DOUBLE : DOUBLE;
939 }
940
TEST(CompareEdgeDistance,Consistency)941 TEST(CompareEdgeDistance, Consistency) {
942 // This test chooses random inputs such that the distance between "x" and
943 // the line (a0, a1) is very close to the threshold distance "r". It then
944 // checks that the answer given by a method at one level of precision is
945 // consistent with the answer given at the next higher level of precision.
946 // See also the comments in the CompareDistances consistency test.
947 auto& rnd = S2Testing::rnd;
948 PrecisionStats stats;
949 for (int iter = 0; iter < FLAGS_consistency_iters; ++iter) {
950 rnd.Reset(iter + 1); // Easier to reproduce a specific case.
951 S2Point a0 = ChoosePoint();
952 S1Angle len = S1Angle::Radians(M_PI * pow(1e-20, rnd.RandDouble()));
953 S2Point a1 = S2::InterpolateAtDistance(len, a0, ChoosePoint());
954 if (rnd.OneIn(2)) a1 = -a1;
955 if (a0 == -a1) continue; // Not allowed by API.
956 S2Point n = S2::RobustCrossProd(a0, a1).Normalize();
957 double f = pow(1e-20, rnd.RandDouble());
958 S2Point a = ((1 - f) * a0 + f * a1).Normalize();
959 S1Angle r = S1Angle::Radians(M_PI_2 * pow(1e-20, rnd.RandDouble()));
960 if (rnd.OneIn(2)) r = S1Angle::Radians(M_PI_2) - r;
961 S2Point x = S2::InterpolateAtDistance(r, a, n);
962 if (rnd.OneIn(5)) {
963 // Replace "x" with a random point that is closest to an edge endpoint.
964 do {
965 x = ChoosePoint();
966 } while (CompareEdgeDirections(a0, x, a0, a1) > 0 &&
967 CompareEdgeDirections(x, a1, a0, a1) > 0);
968 r = min(S1Angle(x, a0), S1Angle(x, a1));
969 }
970 Precision prec = TestCompareEdgeDistanceConsistency(x, a0, a1,
971 S1ChordAngle(r));
972 stats.Tally(prec);
973 }
974 S2_LOG(ERROR) << stats.ToString();
975 }
976
977 // Verifies that CompareEdgeDirections(a0, a1, b0, b1) == expected_sign, and
978 // furthermore checks that the minimum required precision is "expected_prec".
TestCompareEdgeDirections(S2Point a0,S2Point a1,S2Point b0,S2Point b1,int expected_sign,Precision expected_prec)979 void TestCompareEdgeDirections(S2Point a0, S2Point a1, S2Point b0, S2Point b1,
980 int expected_sign, Precision expected_prec) {
981 // Don't normalize the arguments unless necessary (to allow testing points
982 // that differ only in magnitude).
983 if (!S2::IsUnitLength(a0)) a0 = a0.Normalize();
984 if (!S2::IsUnitLength(a1)) a1 = a1.Normalize();
985 if (!S2::IsUnitLength(b0)) b0 = b0.Normalize();
986 if (!S2::IsUnitLength(b1)) b1 = b1.Normalize();
987
988 int dbl_sign = TriageCompareEdgeDirections(a0, a1, b0, b1);
989 int ld_sign = TriageCompareEdgeDirections(ToLD(a0), ToLD(a1),
990 ToLD(b0), ToLD(b1));
991 int exact_sign = ExactCompareEdgeDirections(ToExact(a0), ToExact(a1),
992 ToExact(b0), ToExact(b1));
993
994 // Check that the signs are correct (if non-zero), and also that if dbl_sign
995 // is non-zero then so is ld_sign, etc.
996 EXPECT_EQ(expected_sign, exact_sign);
997 if (ld_sign != 0) EXPECT_EQ(exact_sign, ld_sign);
998 if (dbl_sign != 0) EXPECT_EQ(ld_sign, dbl_sign);
999
1000 Precision actual_prec = (dbl_sign ? DOUBLE : ld_sign ? LONG_DOUBLE : EXACT);
1001 EXPECT_EQ(expected_prec, actual_prec);
1002
1003 // Make sure that the top-level function returns the expected result.
1004 EXPECT_EQ(expected_sign, CompareEdgeDirections(a0, a1, b0, b1));
1005
1006 // Check various identities involving swapping or negating arguments.
1007 EXPECT_EQ(expected_sign, CompareEdgeDirections(b0, b1, a0, a1));
1008 EXPECT_EQ(expected_sign, CompareEdgeDirections(-a0, -a1, b0, b1));
1009 EXPECT_EQ(expected_sign, CompareEdgeDirections(a0, a1, -b0, -b1));
1010 EXPECT_EQ(-expected_sign, CompareEdgeDirections(a1, a0, b0, b1));
1011 EXPECT_EQ(-expected_sign, CompareEdgeDirections(a0, a1, b1, b0));
1012 EXPECT_EQ(-expected_sign, CompareEdgeDirections(-a0, a1, b0, b1));
1013 EXPECT_EQ(-expected_sign, CompareEdgeDirections(a0, -a1, b0, b1));
1014 EXPECT_EQ(-expected_sign, CompareEdgeDirections(a0, a1, -b0, b1));
1015 EXPECT_EQ(-expected_sign, CompareEdgeDirections(a0, a1, b0, -b1));
1016 }
1017
TEST(CompareEdgeDirections,Coverage)1018 TEST(CompareEdgeDirections, Coverage) {
1019 TestCompareEdgeDirections(S2Point(1, 0, 0), S2Point(1, 1, 0),
1020 S2Point(1, -1, 0), S2Point(1, 0, 0),
1021 1, DOUBLE);
1022 TestCompareEdgeDirections(S2Point(1, 0, 1.5e-15), S2Point(1, 1, 0),
1023 S2Point(0, -1, 0), S2Point(0, 0, 1),
1024 1, DOUBLE);
1025 TestCompareEdgeDirections(S2Point(1, 0, 1e-18), S2Point(1, 1, 0),
1026 S2Point(0, -1, 0), S2Point(0, 0, 1),
1027 1, LONG_DOUBLE);
1028 TestCompareEdgeDirections(S2Point(1, 0, 1e-50), S2Point(1, 1, 0),
1029 S2Point(0, -1, 0), S2Point(0, 0, 1),
1030 1, EXACT);
1031 TestCompareEdgeDirections(S2Point(1, 0, 0), S2Point(1, 1, 0),
1032 S2Point(0, -1, 0), S2Point(0, 0, 1),
1033 0, EXACT);
1034 }
1035
1036 // Checks that the result at one level of precision is consistent with the
1037 // result at the next higher level of precision. Returns the minimum
1038 // precision that yielded a non-zero result.
TestCompareEdgeDirectionsConsistency(const S2Point & a0,const S2Point & a1,const S2Point & b0,const S2Point & b1)1039 Precision TestCompareEdgeDirectionsConsistency(
1040 const S2Point& a0, const S2Point& a1,
1041 const S2Point& b0, const S2Point& b1) {
1042 int dbl_sign = TriageCompareEdgeDirections(a0, a1, b0, b1);
1043 int ld_sign = TriageCompareEdgeDirections(ToLD(a0), ToLD(a1),
1044 ToLD(b0), ToLD(b1));
1045 int exact_sign = ExactCompareEdgeDirections(ToExact(a0), ToExact(a1),
1046 ToExact(b0), ToExact(b1));
1047 EXPECT_EQ(exact_sign, CompareEdgeDirections(a0, a1, b0, b1));
1048 if (dbl_sign != 0) EXPECT_EQ(ld_sign, dbl_sign);
1049 if (ld_sign != 0) EXPECT_EQ(exact_sign, ld_sign);
1050 return (ld_sign == 0) ? EXACT : (dbl_sign == 0) ? LONG_DOUBLE : DOUBLE;
1051 }
1052
TEST(CompareEdgeDirections,Consistency)1053 TEST(CompareEdgeDirections, Consistency) {
1054 // This test chooses random pairs of edges that are nearly perpendicular,
1055 // then checks that the answer given by a method at one level of precision
1056 // is consistent with the answer given at the next higher level of
1057 // precision. See also the comments in the CompareDistances test.
1058 auto& rnd = S2Testing::rnd;
1059 PrecisionStats stats;
1060 for (int iter = 0; iter < FLAGS_consistency_iters; ++iter) {
1061 rnd.Reset(iter + 1); // Easier to reproduce a specific case.
1062 S2Point a0 = ChoosePoint();
1063 S1Angle a_len = S1Angle::Radians(M_PI * pow(1e-20, rnd.RandDouble()));
1064 S2Point a1 = S2::InterpolateAtDistance(a_len, a0, ChoosePoint());
1065 S2Point a_norm = S2::RobustCrossProd(a0, a1).Normalize();
1066 S2Point b0 = ChoosePoint();
1067 S1Angle b_len = S1Angle::Radians(M_PI * pow(1e-20, rnd.RandDouble()));
1068 S2Point b1 = S2::InterpolateAtDistance(b_len, b0, a_norm);
1069 if (a0 == -a1 || b0 == -b1) continue; // Not allowed by API.
1070 Precision prec = TestCompareEdgeDirectionsConsistency(a0, a1, b0, b1);
1071 // Don't skew the statistics by recording degenerate inputs.
1072 if (a0 == a1 || b0 == b1) {
1073 EXPECT_EQ(EXACT, prec);
1074 } else {
1075 stats.Tally(prec);
1076 }
1077 }
1078 S2_LOG(ERROR) << stats.ToString();
1079 }
1080
1081 // Verifies that EdgeCircumcenterSign(x0, x1, a, b, c) == expected_sign, and
1082 // furthermore checks that the minimum required precision is "expected_prec".
TestEdgeCircumcenterSign(S2Point x0,S2Point x1,S2Point a,S2Point b,S2Point c,int expected_sign,Precision expected_prec)1083 void TestEdgeCircumcenterSign(
1084 S2Point x0, S2Point x1, S2Point a, S2Point b, S2Point c,
1085 int expected_sign, Precision expected_prec) {
1086 // Don't normalize the arguments unless necessary (to allow testing points
1087 // that differ only in magnitude).
1088 if (!S2::IsUnitLength(x0)) x0 = x0.Normalize();
1089 if (!S2::IsUnitLength(x1)) x1 = x1.Normalize();
1090 if (!S2::IsUnitLength(a)) a = a.Normalize();
1091 if (!S2::IsUnitLength(b)) b = b.Normalize();
1092 if (!S2::IsUnitLength(c)) c = c.Normalize();
1093
1094 int abc_sign = Sign(a, b, c);
1095 int dbl_sign = TriageEdgeCircumcenterSign(x0, x1, a, b, c, abc_sign);
1096 int ld_sign = TriageEdgeCircumcenterSign(
1097 ToLD(x0), ToLD(x1), ToLD(a), ToLD(b), ToLD(c), abc_sign);
1098 int exact_sign = ExactEdgeCircumcenterSign(
1099 ToExact(x0), ToExact(x1), ToExact(a), ToExact(b), ToExact(c), abc_sign);
1100 int actual_sign = (exact_sign != 0 ? exact_sign :
1101 SymbolicEdgeCircumcenterSign(x0, x1, a, b, c));
1102
1103 // Check that the signs are correct (if non-zero), and also that if dbl_sign
1104 // is non-zero then so is ld_sign, etc.
1105 EXPECT_EQ(expected_sign, actual_sign);
1106 if (exact_sign != 0) EXPECT_EQ(exact_sign, actual_sign);
1107 if (ld_sign != 0) EXPECT_EQ(exact_sign, ld_sign);
1108 if (dbl_sign != 0) EXPECT_EQ(ld_sign, dbl_sign);
1109
1110 Precision actual_prec = (dbl_sign ? DOUBLE :
1111 ld_sign ? LONG_DOUBLE :
1112 exact_sign ? EXACT : SYMBOLIC);
1113 EXPECT_EQ(expected_prec, actual_prec);
1114
1115 // Make sure that the top-level function returns the expected result.
1116 EXPECT_EQ(expected_sign, EdgeCircumcenterSign(x0, x1, a, b, c));
1117
1118 // Check various identities involving swapping or negating arguments.
1119 EXPECT_EQ(expected_sign, EdgeCircumcenterSign(x0, x1, a, c, b));
1120 EXPECT_EQ(expected_sign, EdgeCircumcenterSign(x0, x1, b, a, c));
1121 EXPECT_EQ(expected_sign, EdgeCircumcenterSign(x0, x1, b, c, a));
1122 EXPECT_EQ(expected_sign, EdgeCircumcenterSign(x0, x1, c, a, b));
1123 EXPECT_EQ(expected_sign, EdgeCircumcenterSign(x0, x1, c, b, a));
1124 EXPECT_EQ(-expected_sign, EdgeCircumcenterSign(x1, x0, a, b, c));
1125 EXPECT_EQ(expected_sign, EdgeCircumcenterSign(-x0, -x1, a, b, c));
1126 if (actual_sign == exact_sign) {
1127 // Negating the input points may not preserve the result when symbolic
1128 // perturbations are used, since -X is not an exact multiple of X.
1129 EXPECT_EQ(-expected_sign, EdgeCircumcenterSign(x0, x1, -a, -b, -c));
1130 }
1131 }
1132
TEST(EdgeCircumcenterSign,Coverage)1133 TEST(EdgeCircumcenterSign, Coverage) {
1134 TestEdgeCircumcenterSign(
1135 S2Point(1, 0, 0), S2Point(1, 1, 0),
1136 S2Point(0, 0, 1), S2Point(1, 0, 1), S2Point(0, 1, 1),
1137 1, DOUBLE);
1138 TestEdgeCircumcenterSign(
1139 S2Point(1, 0, 0), S2Point(1, 1, 0),
1140 S2Point(0, 0, -1), S2Point(1, 0, -1), S2Point(0, 1, -1),
1141 -1, DOUBLE);
1142 TestEdgeCircumcenterSign(
1143 S2Point(1, -1, 0), S2Point(1, 1, 0),
1144 S2Point(1, -1e-5, 1), S2Point(1, 1e-5, -1), S2Point(1, 1 - 1e-5, 1e-5),
1145 -1, DOUBLE);
1146 TestEdgeCircumcenterSign(
1147 S2Point(1, -1, 0), S2Point(1, 1, 0),
1148 S2Point(1, -1e-5, 1), S2Point(1, 1e-5, -1), S2Point(1, 1 - 1e-9, 1e-5),
1149 -1, LONG_DOUBLE);
1150 TestEdgeCircumcenterSign(
1151 S2Point(1, -1, 0), S2Point(1, 1, 0),
1152 S2Point(1, -1e-5, 1), S2Point(1, 1e-5, -1), S2Point(1, 1 - 1e-15, 1e-5),
1153 -1, EXACT);
1154 TestEdgeCircumcenterSign(
1155 S2Point(1, -1, 0), S2Point(1, 1, 0),
1156 S2Point(1, -1e-5, 1), S2Point(1, 1e-5, -1), S2Point(1, 1, 1e-5),
1157 1, SYMBOLIC);
1158
1159 // This test falls back to the second symbolic perturbation:
1160 TestEdgeCircumcenterSign(
1161 S2Point(1, -1, 0), S2Point(1, 1, 0),
1162 S2Point(0, -1, 0), S2Point(0, 0, -1), S2Point(0, 0, 1),
1163 -1, SYMBOLIC);
1164
1165 // This test falls back to the third symbolic perturbation:
1166 TestEdgeCircumcenterSign(
1167 S2Point(0, -1, 1), S2Point(0, 1, 1),
1168 S2Point(0, 1, 0), S2Point(0, -1, 0), S2Point(1, 0, 0),
1169 -1, SYMBOLIC);
1170 }
1171
1172 // Checks that the result at one level of precision is consistent with the
1173 // result at the next higher level of precision. Returns the minimum
1174 // precision that yielded a non-zero result.
TestEdgeCircumcenterSignConsistency(const S2Point & x0,const S2Point & x1,const S2Point & a,const S2Point & b,const S2Point & c)1175 Precision TestEdgeCircumcenterSignConsistency(
1176 const S2Point& x0, const S2Point& x1,
1177 const S2Point& a, const S2Point& b, const S2Point& c) {
1178 int abc_sign = Sign(a, b, c);
1179 int dbl_sign = TriageEdgeCircumcenterSign(x0, x1, a, b, c, abc_sign);
1180 int ld_sign = TriageEdgeCircumcenterSign(
1181 ToLD(x0), ToLD(x1), ToLD(a), ToLD(b), ToLD(c), abc_sign);
1182 int exact_sign = ExactEdgeCircumcenterSign(
1183 ToExact(x0), ToExact(x1), ToExact(a), ToExact(b), ToExact(c), abc_sign);
1184 if (dbl_sign != 0) EXPECT_EQ(ld_sign, dbl_sign);
1185 if (ld_sign != 0) EXPECT_EQ(exact_sign, ld_sign);
1186 if (exact_sign != 0) {
1187 EXPECT_EQ(exact_sign, EdgeCircumcenterSign(x0, x1, a, b, c));
1188 return (ld_sign == 0) ? EXACT : (dbl_sign == 0) ? LONG_DOUBLE : DOUBLE;
1189 } else {
1190 // Unlike the other methods, SymbolicEdgeCircumcenterSign has the
1191 // precondition that the exact sign must be zero.
1192 int symbolic_sign = SymbolicEdgeCircumcenterSign(x0, x1, a, b, c);
1193 EXPECT_EQ(symbolic_sign, EdgeCircumcenterSign(x0, x1, a, b, c));
1194 return SYMBOLIC;
1195 }
1196 }
1197
TEST(EdgeCircumcenterSign,Consistency)1198 TEST(EdgeCircumcenterSign, Consistency) {
1199 // This test chooses random a random edge X, then chooses a random point Z
1200 // on the great circle through X, and finally choose three points A, B, C
1201 // that are nearly equidistant from X. It then checks that the answer given
1202 // by a method at one level of precision is consistent with the answer given
1203 // at the next higher level of precision.
1204 auto& rnd = S2Testing::rnd;
1205 PrecisionStats stats;
1206 for (int iter = 0; iter < FLAGS_consistency_iters; ++iter) {
1207 rnd.Reset(iter + 1); // Easier to reproduce a specific case.
1208 S2Point x0 = ChoosePoint();
1209 S2Point x1 = ChoosePoint();
1210 if (x0 == -x1) continue; // Not allowed by API.
1211 double c0 = (rnd.OneIn(2) ? -1 : 1) * pow(1e-20, rnd.RandDouble());
1212 double c1 = (rnd.OneIn(2) ? -1 : 1) * pow(1e-20, rnd.RandDouble());
1213 S2Point z = (c0 * x0 + c1 * x1).Normalize();
1214 S1Angle r = S1Angle::Radians(M_PI * pow(1e-30, rnd.RandDouble()));
1215 S2Point a = S2::InterpolateAtDistance(r, z, ChoosePoint());
1216 S2Point b = S2::InterpolateAtDistance(r, z, ChoosePoint());
1217 S2Point c = S2::InterpolateAtDistance(r, z, ChoosePoint());
1218 Precision prec = TestEdgeCircumcenterSignConsistency(x0, x1, a, b, c);
1219 // Don't skew the statistics by recording degenerate inputs.
1220 if (x0 == x1) {
1221 // This precision would be SYMBOLIC if we handled this degeneracy.
1222 EXPECT_EQ(EXACT, prec);
1223 } else if (a == b || b == c || c == a) {
1224 EXPECT_EQ(SYMBOLIC, prec);
1225 } else {
1226 stats.Tally(prec);
1227 }
1228 }
1229 S2_LOG(ERROR) << stats.ToString();
1230 }
1231
1232 // Verifies that VoronoiSiteExclusion(a, b, x0, x1, r) == expected_result, and
1233 // furthermore checks that the minimum required precision is "expected_prec".
TestVoronoiSiteExclusion(S2Point a,S2Point b,S2Point x0,S2Point x1,S1ChordAngle r,Excluded expected_result,Precision expected_prec)1234 void TestVoronoiSiteExclusion(
1235 S2Point a, S2Point b, S2Point x0, S2Point x1, S1ChordAngle r,
1236 Excluded expected_result, Precision expected_prec) {
1237 constexpr Excluded UNCERTAIN = Excluded::UNCERTAIN;
1238
1239 // Don't normalize the arguments unless necessary (to allow testing points
1240 // that differ only in magnitude).
1241 if (!S2::IsUnitLength(a)) a = a.Normalize();
1242 if (!S2::IsUnitLength(b)) b = b.Normalize();
1243 if (!S2::IsUnitLength(x0)) x0 = x0.Normalize();
1244 if (!S2::IsUnitLength(x1)) x1 = x1.Normalize();
1245
1246 // The internal methods (Triage, Exact, etc) require that site A is closer
1247 // to X0 and site B is closer to X1. GetVoronoiSiteExclusion has special
1248 // code to handle the case where this is not true. We need to duplicate
1249 // that code here. Essentially, since the API requires site A to be closer
1250 // than site B to X0, then if site A is also closer to X1 then site B must
1251 // be excluded.
1252 if (s2pred::CompareDistances(x1, a, b) < 0) {
1253 EXPECT_EQ(expected_result, Excluded::SECOND);
1254 // We don't know what precision was used by CompareDistances(), but we
1255 // arbitrarily require the test to specify it as DOUBLE.
1256 EXPECT_EQ(expected_prec, DOUBLE);
1257 } else {
1258 Excluded dbl_result = TriageVoronoiSiteExclusion(a, b, x0, x1, r.length2());
1259 Excluded ld_result = TriageVoronoiSiteExclusion(
1260 ToLD(a), ToLD(b), ToLD(x0), ToLD(x1), ToLD(r.length2()));
1261 Excluded exact_result = ExactVoronoiSiteExclusion(
1262 ToExact(a), ToExact(b), ToExact(x0), ToExact(x1), r.length2());
1263
1264 // Check that the results are correct (if not UNCERTAIN), and also that if
1265 // dbl_result is not UNCERTAIN then so is ld_result, etc.
1266 EXPECT_EQ(expected_result, exact_result);
1267 if (ld_result != UNCERTAIN) EXPECT_EQ(exact_result, ld_result);
1268 if (dbl_result != UNCERTAIN) EXPECT_EQ(ld_result, dbl_result);
1269
1270 Precision actual_prec = (dbl_result != UNCERTAIN ? DOUBLE :
1271 ld_result != UNCERTAIN ? LONG_DOUBLE : EXACT);
1272 EXPECT_EQ(expected_prec, actual_prec);
1273 }
1274 // Make sure that the top-level function returns the expected result.
1275 EXPECT_EQ(expected_result, GetVoronoiSiteExclusion(a, b, x0, x1, r));
1276
1277 // If site B is closer to X1, then the same site should be excluded (if any)
1278 // when we swap the sites and the edge direction.
1279 Excluded swapped_result =
1280 expected_result == Excluded::FIRST ? Excluded::SECOND :
1281 expected_result == Excluded::SECOND ? Excluded::FIRST : expected_result;
1282 if (s2pred::CompareDistances(x1, b, a) < 0) {
1283 EXPECT_EQ(swapped_result, GetVoronoiSiteExclusion(b, a, x1, x0, r));
1284 }
1285 }
1286
TEST(VoronoiSiteExclusion,Coverage)1287 TEST(VoronoiSiteExclusion, Coverage) {
1288 // Both sites are closest to edge endpoint X0.
1289 TestVoronoiSiteExclusion(
1290 S2Point(1, -1e-5, 0), S2Point(1, -2e-5, 0),
1291 S2Point(1, 0, 0), S2Point(1, 1, 0), S1ChordAngle::Radians(1e-3),
1292 Excluded::SECOND, DOUBLE);
1293
1294 // Both sites are closest to edge endpoint X1.
1295 TestVoronoiSiteExclusion(
1296 S2Point(1, 1, 1e-30), S2Point(1, 1, -1e-20),
1297 S2Point(1, 0, 0), S2Point(1, 1, 0), S1ChordAngle::Radians(1e-10),
1298 Excluded::SECOND, DOUBLE);
1299
1300 // Test cases where neither site is excluded.
1301 TestVoronoiSiteExclusion(
1302 S2Point(1, -1e-10, 1e-5), S2Point(1, 1e-10, -1e-5),
1303 S2Point(1, -1, 0), S2Point(1, 1, 0), S1ChordAngle::Radians(1e-4),
1304 Excluded::NEITHER, DOUBLE);
1305 TestVoronoiSiteExclusion(
1306 S2Point(1, -1e-10, 1e-5), S2Point(1, 1e-10, -1e-5),
1307 S2Point(1, -1, 0), S2Point(1, 1, 0), S1ChordAngle::Radians(1e-5),
1308 Excluded::NEITHER, LONG_DOUBLE);
1309 TestVoronoiSiteExclusion(
1310 S2Point(1, -1e-17, 1e-5), S2Point(1, 1e-17, -1e-5),
1311 S2Point(1, -1, 0), S2Point(1, 1, 0), S1ChordAngle::Radians(1e-4),
1312 Excluded::NEITHER, LONG_DOUBLE);
1313 TestVoronoiSiteExclusion(
1314 S2Point(1, -1e-20, 1e-5), S2Point(1, 1e-20, -1e-5),
1315 S2Point(1, -1, 0), S2Point(1, 1, 0), S1ChordAngle::Radians(1e-5),
1316 Excluded::NEITHER, EXACT);
1317
1318 // Test cases where the first site is excluded. (Tests where the second
1319 // site is excluded are constructed by TestVoronoiSiteExclusion.)
1320 TestVoronoiSiteExclusion(
1321 S2Point(1, -1e-6, 1.0049999999e-5), S2Point(1, 0, -1e-5),
1322 S2Point(1, -1, 0), S2Point(1, 1, 0), S1ChordAngle::Radians(1.005e-5),
1323 Excluded::FIRST, DOUBLE);
1324 TestVoronoiSiteExclusion(
1325 S2Point(1, -1.00105e-6, 1.0049999999e-5), S2Point(1, 0, -1e-5),
1326 S2Point(1, -1, 0), S2Point(1, 1, 0), S1ChordAngle::Radians(1.005e-5),
1327 Excluded::FIRST, LONG_DOUBLE);
1328 TestVoronoiSiteExclusion(
1329 S2Point(1, -1e-6, 1.005e-5), S2Point(1, 0, -1e-5),
1330 S2Point(1, -1, 0), S2Point(1, 1, 0), S1ChordAngle::Radians(1.005e-5),
1331 Excluded::FIRST, LONG_DOUBLE);
1332 TestVoronoiSiteExclusion(
1333 S2Point(1, -1e-31, 1.005e-30), S2Point(1, 0, -1e-30),
1334 S2Point(1, -1, 0), S2Point(1, 1, 0), S1ChordAngle::Radians(1.005e-30),
1335 Excluded::FIRST, EXACT);
1336 TestVoronoiSiteExclusion(
1337 S2Point(1, -1e-31, 1.005e-30), S2Point(1, 0, -1e-30),
1338 S2Point(1, -1, 0), S2Point(1, 1, 0), S1ChordAngle::Radians(1.005e-30),
1339 Excluded::FIRST, EXACT);
1340
1341 // These two sites are exactly 60 degrees away from the point (1, 1, 0),
1342 // which is the midpoint of edge X. This case requires symbolic
1343 // perturbations to resolve correctly. Site A is closer to every point in
1344 // its coverage interval except for (1, 1, 0), but site B is considered
1345 // closer to that point symbolically.
1346 TestVoronoiSiteExclusion(
1347 S2Point(0, 1, 1), S2Point(1, 0, 1),
1348 S2Point(0, 1, 1), S2Point(1, 0, -1), S1ChordAngle::FromLength2(1),
1349 Excluded::NEITHER, EXACT);
1350
1351 // This test is similar except that site A is considered closer to the
1352 // equidistant point (-1, 1, 0), and therefore site B is excluded.
1353 TestVoronoiSiteExclusion(
1354 S2Point(0, 1, 1), S2Point(-1, 0, 1),
1355 S2Point(0, 1, 1), S2Point(-1, 0, -1), S1ChordAngle::FromLength2(1),
1356 Excluded::SECOND, EXACT);
1357 }
1358
1359 // Checks that the result at one level of precision is consistent with the
1360 // result at the next higher level of precision. Returns the minimum
1361 // precision that yielded a non-zero result.
TestVoronoiSiteExclusionConsistency(const S2Point & a,const S2Point & b,const S2Point & x0,const S2Point & x1,S1ChordAngle r)1362 Precision TestVoronoiSiteExclusionConsistency(
1363 const S2Point& a, const S2Point& b, const S2Point& x0, const S2Point& x1,
1364 S1ChordAngle r) {
1365 constexpr Excluded UNCERTAIN = Excluded::UNCERTAIN;
1366
1367 // The internal methods require this (see TestVoronoiSiteExclusion).
1368 if (s2pred::CompareDistances(x1, a, b) < 0) return DOUBLE;
1369
1370 Excluded dbl_result = TriageVoronoiSiteExclusion(a, b, x0, x1, r.length2());
1371 Excluded ld_result = TriageVoronoiSiteExclusion(
1372 ToLD(a), ToLD(b), ToLD(x0), ToLD(x1), ToLD(r.length2()));
1373 Excluded exact_result = ExactVoronoiSiteExclusion(
1374 ToExact(a), ToExact(b), ToExact(x0), ToExact(x1), r.length2());
1375 EXPECT_EQ(exact_result, GetVoronoiSiteExclusion(a, b, x0, x1, r));
1376
1377 EXPECT_NE(UNCERTAIN, exact_result);
1378 if (ld_result == UNCERTAIN) {
1379 EXPECT_EQ(UNCERTAIN, dbl_result);
1380 return EXACT;
1381 }
1382 EXPECT_EQ(exact_result, ld_result);
1383 if (dbl_result == UNCERTAIN) {
1384 return LONG_DOUBLE;
1385 }
1386 EXPECT_EQ(exact_result, dbl_result);
1387 return DOUBLE;
1388 }
1389
TEST(VoronoiSiteExclusion,Consistency)1390 TEST(VoronoiSiteExclusion, Consistency) {
1391 // This test chooses random a random edge X, a random point P on that edge,
1392 // and a random threshold distance "r". It then choose two sites A and B
1393 // whose distance to P is almost exactly "r". This ensures that the
1394 // coverage intervals for A and B will (almost) share a common endpoint. It
1395 // then checks that the answer given by a method at one level of precision
1396 // is consistent with the answer given at higher levels of precision.
1397 auto& rnd = S2Testing::rnd;
1398 PrecisionStats stats;
1399 for (int iter = 0; iter < FLAGS_consistency_iters; ++iter) {
1400 rnd.Reset(iter + 1); // Easier to reproduce a specific case.
1401 S2Point x0 = ChoosePoint();
1402 S2Point x1 = ChoosePoint();
1403 if (x0 == -x1) continue; // Not allowed by API.
1404 double f = pow(1e-20, rnd.RandDouble());
1405 S2Point p = ((1 - f) * x0 + f * x1).Normalize();
1406 S1Angle r1 = S1Angle::Radians(M_PI_2 * pow(1e-20, rnd.RandDouble()));
1407 S2Point a = S2::InterpolateAtDistance(r1, p, ChoosePoint());
1408 S2Point b = S2::InterpolateAtDistance(r1, p, ChoosePoint());
1409 // Check that the other API requirements are met.
1410 S1ChordAngle r(r1);
1411 if (s2pred::CompareEdgeDistance(a, x0, x1, r) > 0) continue;
1412 if (s2pred::CompareEdgeDistance(b, x0, x1, r) > 0) continue;
1413 if (s2pred::CompareDistances(x0, a, b) > 0) std::swap(a, b);
1414 if (a == b) continue;
1415
1416 Precision prec = TestVoronoiSiteExclusionConsistency(a, b, x0, x1, r);
1417 // Don't skew the statistics by recording degenerate inputs.
1418 if (x0 == x1) {
1419 EXPECT_EQ(DOUBLE, prec);
1420 } else {
1421 stats.Tally(prec);
1422 }
1423 }
1424 S2_LOG(ERROR) << stats.ToString();
1425 }
1426
1427 } // namespace s2pred
1428