4 files changed, 852 insertions, 0 deletions
diff --git a/vendor/github.com/golang/geo/s1/angle.go b/vendor/github.com/golang/geo/s1/angle.go
new file mode 100644
index 000000000..747b23dea
--- /dev/null
+++ b/vendor/github.com/golang/geo/s1/angle.go
@@ -0,0 +1,120 @@
+// Copyright 2014 Google Inc. All rights reserved.
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+package s1
+
+import (
+	"math"
+	"strconv"
+)
+
+// Angle represents a 1D angle. The internal representation is a double precision
+// value in radians, so conversion to and from radians is exact.
+// Conversions between E5, E6, E7, and Degrees are not always
+// exact. For example, Degrees(3.1) is different from E6(3100000) or E7(31000000).
+//
+// The following conversions between degrees and radians are exact:
+//
+//       Degree*180 == Radian*math.Pi
+//   Degree*(180/n) == Radian*(math.Pi/n)     for n == 0..8
+//
+// These identities hold when the arguments are scaled up or down by any power
+// of 2. Some similar identities are also true, for example,
+//
+//   Degree*60 == Radian*(math.Pi/3)
+//
+// But be aware that this type of identity does not hold in general. For example,
+//
+//   Degree*3 != Radian*(math.Pi/60)
+//
+// Similarly, the conversion to radians means that (Angle(x)*Degree).Degrees()
+// does not always equal x. For example,
+//
+//   (Angle(45*n)*Degree).Degrees() == 45*n     for n == 0..8
+//
+// but
+//
+//   (60*Degree).Degrees() != 60
+//
+// When testing for equality, you should allow for numerical errors (ApproxEqual)
+// or convert to discrete E5/E6/E7 values first.
+type Angle float64
+
+// Angle units.
+const (
+	Radian Angle = 1
+	Degree       = (math.Pi / 180) * Radian
+
+	E5 = 1e-5 * Degree
+	E6 = 1e-6 * Degree
+	E7 = 1e-7 * Degree
+)
+
+// Radians returns the angle in radians.
+func (a Angle) Radians() float64 { return float64(a) }
+
+// Degrees returns the angle in degrees.
+func (a Angle) Degrees() float64 { return float64(a / Degree) }
+
+// round returns the value rounded to nearest as an int32.
+// This does not match C++ exactly for the case of x.5.
+func round(val float64) int32 {
+	if val < 0 {
+		return int32(val - 0.5)
+	}
+	return int32(val + 0.5)
+}
+
+// InfAngle returns an angle larger than any finite angle.
+func InfAngle() Angle {
+	return Angle(math.Inf(1))
+}
+
+// isInf reports whether this Angle is infinite.
+func (a Angle) isInf() bool {
+	return math.IsInf(float64(a), 0)
+}
+
+// E5 returns the angle in hundred thousandths of degrees.
+func (a Angle) E5() int32 { return round(a.Degrees() * 1e5) }
+
+// E6 returns the angle in millionths of degrees.
+func (a Angle) E6() int32 { return round(a.Degrees() * 1e6) }
+
+// E7 returns the angle in ten millionths of degrees.
+func (a Angle) E7() int32 { return round(a.Degrees() * 1e7) }
+
+// Abs returns the absolute value of the angle.
+func (a Angle) Abs() Angle { return Angle(math.Abs(float64(a))) }
+
+// Normalized returns an equivalent angle in (-π, π].
+func (a Angle) Normalized() Angle {
+	rad := math.Remainder(float64(a), 2*math.Pi)
+	if rad <= -math.Pi {
+		rad = math.Pi
+	}
+	return Angle(rad)
+}
+
+func (a Angle) String() string {
+	return strconv.FormatFloat(a.Degrees(), 'f', 7, 64) // like "%.7f"
+}
+
+// ApproxEqual reports whether the two angles are the same up to a small tolerance.
+func (a Angle) ApproxEqual(other Angle) bool {
+	return math.Abs(float64(a)-float64(other)) <= epsilon
+}
+
+// BUG(dsymonds): The major differences from the C++ version are:
+//   - no unsigned E5/E6/E7 methods
diff --git a/vendor/github.com/golang/geo/s1/chordangle.go b/vendor/github.com/golang/geo/s1/chordangle.go
new file mode 100644
index 000000000..406c69ef1
--- /dev/null
+++ b/vendor/github.com/golang/geo/s1/chordangle.go
@@ -0,0 +1,250 @@
+// Copyright 2015 Google Inc. All rights reserved.
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+package s1
+
+import (
+	"math"
+)
+
+// ChordAngle represents the angle subtended by a chord (i.e., the straight
+// line segment connecting two points on the sphere). Its representation
+// makes it very efficient for computing and comparing distances, but unlike
+// Angle it is only capable of representing angles between 0 and π radians.
+// Generally, ChordAngle should only be used in loops where many angles need
+// to be calculated and compared. Otherwise it is simpler to use Angle.
+//
+// ChordAngle loses some accuracy as the angle approaches π radians.
+// Specifically, the representation of (π - x) radians has an error of about
+// (1e-15 / x), with a maximum error of about 2e-8 radians (about 13cm on the
+// Earth's surface). For comparison, for angles up to π/2 radians (10000km)
+// the worst-case representation error is about 2e-16 radians (1 nanonmeter),
+// which is about the same as Angle.
+//
+// ChordAngles are represented by the squared chord length, which can
+// range from 0 to 4. Positive infinity represents an infinite squared length.
+type ChordAngle float64
+
+const (
+	// NegativeChordAngle represents a chord angle smaller than the zero angle.
+	// The only valid operations on a NegativeChordAngle are comparisons,
+	// Angle conversions, and Successor/Predecessor.
+	NegativeChordAngle = ChordAngle(-1)
+
+	// RightChordAngle represents a chord angle of 90 degrees (a "right angle").
+	RightChordAngle = ChordAngle(2)
+
+	// StraightChordAngle represents a chord angle of 180 degrees (a "straight angle").
+	// This is the maximum finite chord angle.
+	StraightChordAngle = ChordAngle(4)
+
+	// maxLength2 is the square of the maximum length allowed in a ChordAngle.
+	maxLength2 = 4.0
+)
+
+// ChordAngleFromAngle returns a ChordAngle from the given Angle.
+func ChordAngleFromAngle(a Angle) ChordAngle {
+	if a < 0 {
+		return NegativeChordAngle
+	}
+	if a.isInf() {
+		return InfChordAngle()
+	}
+	l := 2 * math.Sin(0.5*math.Min(math.Pi, a.Radians()))
+	return ChordAngle(l * l)
+}
+
+// ChordAngleFromSquaredLength returns a ChordAngle from the squared chord length.
+// Note that the argument is automatically clamped to a maximum of 4 to
+// handle possible roundoff errors. The argument must be non-negative.
+func ChordAngleFromSquaredLength(length2 float64) ChordAngle {
+	if length2 > maxLength2 {
+		return StraightChordAngle
+	}
+	return ChordAngle(length2)
+}
+
+// Expanded returns a new ChordAngle that has been adjusted by the given error
+// bound (which can be positive or negative). Error should be the value
+// returned by either MaxPointError or MaxAngleError. For example:
+//    a := ChordAngleFromPoints(x, y)
+//    a1 := a.Expanded(a.MaxPointError())
+func (c ChordAngle) Expanded(e float64) ChordAngle {
+	// If the angle is special, don't change it. Otherwise clamp it to the valid range.
+	if c.isSpecial() {
+		return c
+	}
+	return ChordAngle(math.Max(0.0, math.Min(maxLength2, float64(c)+e)))
+}
+
+// Angle converts this ChordAngle to an Angle.
+func (c ChordAngle) Angle() Angle {
+	if c < 0 {
+		return -1 * Radian
+	}
+	if c.isInf() {
+		return InfAngle()
+	}
+	return Angle(2 * math.Asin(0.5*math.Sqrt(float64(c))))
+}
+
+// InfChordAngle returns a chord angle larger than any finite chord angle.
+// The only valid operations on an InfChordAngle are comparisons, Angle
+// conversions, and Successor/Predecessor.
+func InfChordAngle() ChordAngle {
+	return ChordAngle(math.Inf(1))
+}
+
+// isInf reports whether this ChordAngle is infinite.
+func (c ChordAngle) isInf() bool {
+	return math.IsInf(float64(c), 1)
+}
+
+// isSpecial reports whether this ChordAngle is one of the special cases.
+func (c ChordAngle) isSpecial() bool {
+	return c < 0 || c.isInf()
+}
+
+// isValid reports whether this ChordAngle is valid or not.
+func (c ChordAngle) isValid() bool {
+	return (c >= 0 && c <= maxLength2) || c.isSpecial()
+}
+
+// Successor returns the smallest representable ChordAngle larger than this one.
+// This can be used to convert a "<" comparison to a "<=" comparison.
+//
+// Note the following special cases:
+//   NegativeChordAngle.Successor == 0
+//   StraightChordAngle.Successor == InfChordAngle
+//   InfChordAngle.Successor == InfChordAngle
+func (c ChordAngle) Successor() ChordAngle {
+	if c >= maxLength2 {
+		return InfChordAngle()
+	}
+	if c < 0 {
+		return 0
+	}
+	return ChordAngle(math.Nextafter(float64(c), 10.0))
+}
+
+// Predecessor returns the largest representable ChordAngle less than this one.
+//
+// Note the following special cases:
+//   InfChordAngle.Predecessor == StraightChordAngle
+//   ChordAngle(0).Predecessor == NegativeChordAngle
+//   NegativeChordAngle.Predecessor == NegativeChordAngle
+func (c ChordAngle) Predecessor() ChordAngle {
+	if c <= 0 {
+		return NegativeChordAngle
+	}
+	if c > maxLength2 {
+		return StraightChordAngle
+	}
+
+	return ChordAngle(math.Nextafter(float64(c), -10.0))
+}
+
+// MaxPointError returns the maximum error size for a ChordAngle constructed
+// from 2 Points x and y, assuming that x and y are normalized to within the
+// bounds guaranteed by s2.Point.Normalize. The error is defined with respect to
+// the true distance after the points are projected to lie exactly on the sphere.
+func (c ChordAngle) MaxPointError() float64 {
+	// There is a relative error of (2.5*dblEpsilon) when computing the squared
+	// distance, plus a relative error of 2 * dblEpsilon, plus an absolute error
+	// of (16 * dblEpsilon**2) because the lengths of the input points may differ
+	// from 1 by up to (2*dblEpsilon) each. (This is the maximum error in Normalize).
+	return 4.5*dblEpsilon*float64(c) + 16*dblEpsilon*dblEpsilon
+}
+
+// MaxAngleError returns the maximum error for a ChordAngle constructed
+// as an Angle distance.
+func (c ChordAngle) MaxAngleError() float64 {
+	return dblEpsilon * float64(c)
+}
+
+// Add adds the other ChordAngle to this one and returns the resulting value.
+// This method assumes the ChordAngles are not special.
+func (c ChordAngle) Add(other ChordAngle) ChordAngle {
+	// Note that this method (and Sub) is much more efficient than converting
+	// the ChordAngle to an Angle and adding those and converting back. It
+	// requires only one square root plus a few additions and multiplications.
+
+	// Optimization for the common case where b is an error tolerance
+	// parameter that happens to be set to zero.
+	if other == 0 {
+		return c
+	}
+
+	// Clamp the angle sum to at most 180 degrees.
+	if c+other >= maxLength2 {
+		return StraightChordAngle
+	}
+
+	// Let a and b be the (non-squared) chord lengths, and let c = a+b.
+	// Let A, B, and C be the corresponding half-angles (a = 2*sin(A), etc).
+	// Then the formula below can be derived from c = 2 * sin(A+B) and the
+	// relationships   sin(A+B) = sin(A)*cos(B) + sin(B)*cos(A)
+	//                 cos(X) = sqrt(1 - sin^2(X))
+	x := float64(c * (1 - 0.25*other))
+	y := float64(other * (1 - 0.25*c))
+	return ChordAngle(math.Min(maxLength2, x+y+2*math.Sqrt(x*y)))
+}
+
+// Sub subtracts the other ChordAngle from this one and returns the resulting
+// value. This method assumes the ChordAngles are not special.
+func (c ChordAngle) Sub(other ChordAngle) ChordAngle {
+	if other == 0 {
+		return c
+	}
+	if c <= other {
+		return 0
+	}
+	x := float64(c * (1 - 0.25*other))
+	y := float64(other * (1 - 0.25*c))
+	return ChordAngle(math.Max(0.0, x+y-2*math.Sqrt(x*y)))
+}
+
+// Sin returns the sine of this chord angle. This method is more efficient
+// than converting to Angle and performing the computation.
+func (c ChordAngle) Sin() float64 {
+	return math.Sqrt(c.Sin2())
+}
+
+// Sin2 returns the square of the sine of this chord angle.
+// It is more efficient than Sin.
+func (c ChordAngle) Sin2() float64 {
+	// Let a be the (non-squared) chord length, and let A be the corresponding
+	// half-angle (a = 2*sin(A)).  The formula below can be derived from:
+	//   sin(2*A) = 2 * sin(A) * cos(A)
+	//   cos^2(A) = 1 - sin^2(A)
+	// This is much faster than converting to an angle and computing its sine.
+	return float64(c * (1 - 0.25*c))
+}
+
+// Cos returns the cosine of this chord angle. This method is more efficient
+// than converting to Angle and performing the computation.
+func (c ChordAngle) Cos() float64 {
+	// cos(2*A) = cos^2(A) - sin^2(A) = 1 - 2*sin^2(A)
+	return float64(1 - 0.5*c)
+}
+
+// Tan returns the tangent of this chord angle.
+func (c ChordAngle) Tan() float64 {
+	return c.Sin() / c.Cos()
+}
+
+// TODO(roberts): Differences from C++:
+//   Helpers to/from E5/E6/E7
+//   Helpers to/from degrees and radians directly.
+//   FastUpperBoundFrom(angle Angle)
diff --git a/vendor/github.com/golang/geo/s1/doc.go b/vendor/github.com/golang/geo/s1/doc.go
new file mode 100644
index 000000000..52a2c526d
--- /dev/null
+++ b/vendor/github.com/golang/geo/s1/doc.go
@@ -0,0 +1,20 @@
+// Copyright 2014 Google Inc. All rights reserved.
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+/*
+Package s1 implements types and functions for working with geometry in S¹ (circular geometry).
+
+See ../s2 for a more detailed overview.
+*/
+package s1
diff --git a/vendor/github.com/golang/geo/s1/interval.go b/vendor/github.com/golang/geo/s1/interval.go
new file mode 100644
index 000000000..6fea5221f
--- /dev/null
+++ b/vendor/github.com/golang/geo/s1/interval.go
@@ -0,0 +1,462 @@
+// Copyright 2014 Google Inc. All rights reserved.
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+package s1
+
+import (
+	"math"
+	"strconv"
+)
+
+// An Interval represents a closed interval on a unit circle (also known
+// as a 1-dimensional sphere). It is capable of representing the empty
+// interval (containing no points), the full interval (containing all
+// points), and zero-length intervals (containing a single point).
+//
+// Points are represented by the angle they make with the positive x-axis in
+// the range [-π, π]. An interval is represented by its lower and upper
+// bounds (both inclusive, since the interval is closed). The lower bound may
+// be greater than the upper bound, in which case the interval is "inverted"
+// (i.e. it passes through the point (-1, 0)).
+//
+// The point (-1, 0) has two valid representations, π and -π. The
+// normalized representation of this point is π, so that endpoints
+// of normal intervals are in the range (-π, π]. We normalize the latter to
+// the former in IntervalFromEndpoints. However, we take advantage of the point
+// -π to construct two special intervals:
+//   The full interval is [-π, π]
+//   The empty interval is [π, -π].
+//
+// Treat the exported fields as read-only.
+type Interval struct {
+	Lo, Hi float64
+}
+
+// IntervalFromEndpoints constructs a new interval from endpoints.
+// Both arguments must be in the range [-π,π]. This function allows inverted intervals
+// to be created.
+func IntervalFromEndpoints(lo, hi float64) Interval {
+	i := Interval{lo, hi}
+	if lo == -math.Pi && hi != math.Pi {
+		i.Lo = math.Pi
+	}
+	if hi == -math.Pi && lo != math.Pi {
+		i.Hi = math.Pi
+	}
+	return i
+}
+
+// IntervalFromPointPair returns the minimal interval containing the two given points.
+// Both arguments must be in [-π,π].
+func IntervalFromPointPair(a, b float64) Interval {
+	if a == -math.Pi {
+		a = math.Pi
+	}
+	if b == -math.Pi {
+		b = math.Pi
+	}
+	if positiveDistance(a, b) <= math.Pi {
+		return Interval{a, b}
+	}
+	return Interval{b, a}
+}
+
+// EmptyInterval returns an empty interval.
+func EmptyInterval() Interval { return Interval{math.Pi, -math.Pi} }
+
+// FullInterval returns a full interval.
+func FullInterval() Interval { return Interval{-math.Pi, math.Pi} }
+
+// IsValid reports whether the interval is valid.
+func (i Interval) IsValid() bool {
+	return (math.Abs(i.Lo) <= math.Pi && math.Abs(i.Hi) <= math.Pi &&
+		!(i.Lo == -math.Pi && i.Hi != math.Pi) &&
+		!(i.Hi == -math.Pi && i.Lo != math.Pi))
+}
+
+// IsFull reports whether the interval is full.
+func (i Interval) IsFull() bool { return i.Lo == -math.Pi && i.Hi == math.Pi }
+
+// IsEmpty reports whether the interval is empty.
+func (i Interval) IsEmpty() bool { return i.Lo == math.Pi && i.Hi == -math.Pi }
+
+// IsInverted reports whether the interval is inverted; that is, whether Lo > Hi.
+func (i Interval) IsInverted() bool { return i.Lo > i.Hi }
+
+// Invert returns the interval with endpoints swapped.
+func (i Interval) Invert() Interval {
+	return Interval{i.Hi, i.Lo}
+}
+
+// Center returns the midpoint of the interval.
+// It is undefined for full and empty intervals.
+func (i Interval) Center() float64 {
+	c := 0.5 * (i.Lo + i.Hi)
+	if !i.IsInverted() {
+		return c
+	}
+	if c <= 0 {
+		return c + math.Pi
+	}
+	return c - math.Pi
+}
+
+// Length returns the length of the interval.
+// The length of an empty interval is negative.
+func (i Interval) Length() float64 {
+	l := i.Hi - i.Lo
+	if l >= 0 {
+		return l
+	}
+	l += 2 * math.Pi
+	if l > 0 {
+		return l
+	}
+	return -1
+}
+
+// Assumes p ∈ (-π,π].
+func (i Interval) fastContains(p float64) bool {
+	if i.IsInverted() {
+		return (p >= i.Lo || p <= i.Hi) && !i.IsEmpty()
+	}
+	return p >= i.Lo && p <= i.Hi
+}
+
+// Contains returns true iff the interval contains p.
+// Assumes p ∈ [-π,π].
+func (i Interval) Contains(p float64) bool {
+	if p == -math.Pi {
+		p = math.Pi
+	}
+	return i.fastContains(p)
+}
+
+// ContainsInterval returns true iff the interval contains oi.
+func (i Interval) ContainsInterval(oi Interval) bool {
+	if i.IsInverted() {
+		if oi.IsInverted() {
+			return oi.Lo >= i.Lo && oi.Hi <= i.Hi
+		}
+		return (oi.Lo >= i.Lo || oi.Hi <= i.Hi) && !i.IsEmpty()
+	}
+	if oi.IsInverted() {
+		return i.IsFull() || oi.IsEmpty()
+	}
+	return oi.Lo >= i.Lo && oi.Hi <= i.Hi
+}
+
+// InteriorContains returns true iff the interior of the interval contains p.
+// Assumes p ∈ [-π,π].
+func (i Interval) InteriorContains(p float64) bool {
+	if p == -math.Pi {
+		p = math.Pi
+	}
+	if i.IsInverted() {
+		return p > i.Lo || p < i.Hi
+	}
+	return (p > i.Lo && p < i.Hi) || i.IsFull()
+}
+
+// InteriorContainsInterval returns true iff the interior of the interval contains oi.
+func (i Interval) InteriorContainsInterval(oi Interval) bool {
+	if i.IsInverted() {
+		if oi.IsInverted() {
+			return (oi.Lo > i.Lo && oi.Hi < i.Hi) || oi.IsEmpty()
+		}
+		return oi.Lo > i.Lo || oi.Hi < i.Hi
+	}
+	if oi.IsInverted() {
+		return i.IsFull() || oi.IsEmpty()
+	}
+	return (oi.Lo > i.Lo && oi.Hi < i.Hi) || i.IsFull()
+}
+
+// Intersects returns true iff the interval contains any points in common with oi.
+func (i Interval) Intersects(oi Interval) bool {
+	if i.IsEmpty() || oi.IsEmpty() {
+		return false
+	}
+	if i.IsInverted() {
+		return oi.IsInverted() || oi.Lo <= i.Hi || oi.Hi >= i.Lo
+	}
+	if oi.IsInverted() {
+		return oi.Lo <= i.Hi || oi.Hi >= i.Lo
+	}
+	return oi.Lo <= i.Hi && oi.Hi >= i.Lo
+}
+
+// InteriorIntersects returns true iff the interior of the interval contains any points in common with oi, including the latter's boundary.
+func (i Interval) InteriorIntersects(oi Interval) bool {
+	if i.IsEmpty() || oi.IsEmpty() || i.Lo == i.Hi {
+		return false
+	}
+	if i.IsInverted() {
+		return oi.IsInverted() || oi.Lo < i.Hi || oi.Hi > i.Lo
+	}
+	if oi.IsInverted() {
+		return oi.Lo < i.Hi || oi.Hi > i.Lo
+	}
+	return (oi.Lo < i.Hi && oi.Hi > i.Lo) || i.IsFull()
+}
+
+// Compute distance from a to b in [0,2π], in a numerically stable way.
+func positiveDistance(a, b float64) float64 {
+	d := b - a
+	if d >= 0 {
+		return d
+	}
+	return (b + math.Pi) - (a - math.Pi)
+}
+
+// Union returns the smallest interval that contains both the interval and oi.
+func (i Interval) Union(oi Interval) Interval {
+	if oi.IsEmpty() {
+		return i
+	}
+	if i.fastContains(oi.Lo) {
+		if i.fastContains(oi.Hi) {
+			// Either oi ⊂ i, or i ∪ oi is the full interval.
+			if i.ContainsInterval(oi) {
+				return i
+			}
+			return FullInterval()
+		}
+		return Interval{i.Lo, oi.Hi}
+	}
+	if i.fastContains(oi.Hi) {
+		return Interval{oi.Lo, i.Hi}
+	}
+
+	// Neither endpoint of oi is in i. Either i ⊂ oi, or i and oi are disjoint.
+	if i.IsEmpty() || oi.fastContains(i.Lo) {
+		return oi
+	}
+
+	// This is the only hard case where we need to find the closest pair of endpoints.
+	if positiveDistance(oi.Hi, i.Lo) < positiveDistance(i.Hi, oi.Lo) {
+		return Interval{oi.Lo, i.Hi}
+	}
+	return Interval{i.Lo, oi.Hi}
+}
+
+// Intersection returns the smallest interval that contains the intersection of the interval and oi.
+func (i Interval) Intersection(oi Interval) Interval {
+	if oi.IsEmpty() {
+		return EmptyInterval()
+	}
+	if i.fastContains(oi.Lo) {
+		if i.fastContains(oi.Hi) {
+			// Either oi ⊂ i, or i and oi intersect twice. Neither are empty.
+			// In the first case we want to return i (which is shorter than oi).
+			// In the second case one of them is inverted, and the smallest interval
+			// that covers the two disjoint pieces is the shorter of i and oi.
+			// We thus want to pick the shorter of i and oi in both cases.
+			if oi.Length() < i.Length() {
+				return oi
+			}
+			return i
+		}
+		return Interval{oi.Lo, i.Hi}
+	}
+	if i.fastContains(oi.Hi) {
+		return Interval{i.Lo, oi.Hi}
+	}
+
+	// Neither endpoint of oi is in i. Either i ⊂ oi, or i and oi are disjoint.
+	if oi.fastContains(i.Lo) {
+		return i
+	}
+	return EmptyInterval()
+}
+
+// AddPoint returns the interval expanded by the minimum amount necessary such
+// that it contains the given point "p" (an angle in the range [-π, π]).
+func (i Interval) AddPoint(p float64) Interval {
+	if math.Abs(p) > math.Pi {
+		return i
+	}
+	if p == -math.Pi {
+		p = math.Pi
+	}
+	if i.fastContains(p) {
+		return i
+	}
+	if i.IsEmpty() {
+		return Interval{p, p}
+	}
+	if positiveDistance(p, i.Lo) < positiveDistance(i.Hi, p) {
+		return Interval{p, i.Hi}
+	}
+	return Interval{i.Lo, p}
+}
+
+// Define the maximum rounding error for arithmetic operations. Depending on the
+// platform the mantissa precision may be different than others, so we choose to
+// use specific values to be consistent across all.
+// The values come from the C++ implementation.
+var (
+	// epsilon is a small number that represents a reasonable level of noise between two
+	// values that can be considered to be equal.
+	epsilon = 1e-15
+	// dblEpsilon is a smaller number for values that require more precision.
+	dblEpsilon = 2.220446049e-16
+)
+
+// Expanded returns an interval that has been expanded on each side by margin.
+// If margin is negative, then the function shrinks the interval on
+// each side by margin instead. The resulting interval may be empty or
+// full. Any expansion (positive or negative) of a full interval remains
+// full, and any expansion of an empty interval remains empty.
+func (i Interval) Expanded(margin float64) Interval {
+	if margin >= 0 {
+		if i.IsEmpty() {
+			return i
+		}
+		// Check whether this interval will be full after expansion, allowing
+		// for a rounding error when computing each endpoint.
+		if i.Length()+2*margin+2*dblEpsilon >= 2*math.Pi {
+			return FullInterval()
+		}
+	} else {
+		if i.IsFull() {
+			return i
+		}
+		// Check whether this interval will be empty after expansion, allowing
+		// for a rounding error when computing each endpoint.
+		if i.Length()+2*margin-2*dblEpsilon <= 0 {
+			return EmptyInterval()
+		}
+	}
+	result := IntervalFromEndpoints(
+		math.Remainder(i.Lo-margin, 2*math.Pi),
+		math.Remainder(i.Hi+margin, 2*math.Pi),
+	)
+	if result.Lo <= -math.Pi {
+		result.Lo = math.Pi
+	}
+	return result
+}
+
+// ApproxEqual reports whether this interval can be transformed into the given
+// interval by moving each endpoint by at most ε, without the
+// endpoints crossing (which would invert the interval). Empty and full
+// intervals are considered to start at an arbitrary point on the unit circle,
+// so any interval with (length <= 2*ε) matches the empty interval, and
+// any interval with (length >= 2*π - 2*ε) matches the full interval.
+func (i Interval) ApproxEqual(other Interval) bool {
+	// Full and empty intervals require special cases because the endpoints
+	// are considered to be positioned arbitrarily.
+	if i.IsEmpty() {
+		return other.Length() <= 2*epsilon
+	}
+	if other.IsEmpty() {
+		return i.Length() <= 2*epsilon
+	}
+	if i.IsFull() {
+		return other.Length() >= 2*(math.Pi-epsilon)
+	}
+	if other.IsFull() {
+		return i.Length() >= 2*(math.Pi-epsilon)
+	}
+
+	// The purpose of the last test below is to verify that moving the endpoints
+	// does not invert the interval, e.g. [-1e20, 1e20] vs. [1e20, -1e20].
+	return (math.Abs(math.Remainder(other.Lo-i.Lo, 2*math.Pi)) <= epsilon &&
+		math.Abs(math.Remainder(other.Hi-i.Hi, 2*math.Pi)) <= epsilon &&
+		math.Abs(i.Length()-other.Length()) <= 2*epsilon)
+
+}
+
+func (i Interval) String() string {
+	// like "[%.7f, %.7f]"
+	return "[" + strconv.FormatFloat(i.Lo, 'f', 7, 64) + ", " + strconv.FormatFloat(i.Hi, 'f', 7, 64) + "]"
+}
+
+// Complement returns the complement of the interior of the interval. An interval and
+// its complement have the same boundary but do not share any interior
+// values. The complement operator is not a bijection, since the complement
+// of a singleton interval (containing a single value) is the same as the
+// complement of an empty interval.
+func (i Interval) Complement() Interval {
+	if i.Lo == i.Hi {
+		// Singleton. The interval just contains a single point.
+		return FullInterval()
+	}
+	// Handles empty and full.
+	return Interval{i.Hi, i.Lo}
+}
+
+// ComplementCenter returns the midpoint of the complement of the interval. For full and empty
+// intervals, the result is arbitrary. For a singleton interval (containing a
+// single point), the result is its antipodal point on S1.
+func (i Interval) ComplementCenter() float64 {
+	if i.Lo != i.Hi {
+		return i.Complement().Center()
+	}
+	// Singleton. The interval just contains a single point.
+	if i.Hi <= 0 {
+		return i.Hi + math.Pi
+	}
+	return i.Hi - math.Pi
+}
+
+// DirectedHausdorffDistance returns the Hausdorff distance to the given interval.
+// For two intervals i and y, this distance is defined by
+//     h(i, y) = max_{p in i} min_{q in y} d(p, q),
+// where d(.,.) is measured along S1.
+func (i Interval) DirectedHausdorffDistance(y Interval) Angle {
+	if y.ContainsInterval(i) {
+		return 0 // This includes the case i is empty.
+	}
+	if y.IsEmpty() {
+		return Angle(math.Pi) // maximum possible distance on s1.
+	}
+	yComplementCenter := y.ComplementCenter()
+	if i.Contains(yComplementCenter) {
+		return Angle(positiveDistance(y.Hi, yComplementCenter))
+	}
+
+	// The Hausdorff distance is realized by either two i.Hi endpoints or two
+	// i.Lo endpoints, whichever is farther apart.
+	hiHi := 0.0
+	if IntervalFromEndpoints(y.Hi, yComplementCenter).Contains(i.Hi) {
+		hiHi = positiveDistance(y.Hi, i.Hi)
+	}
+
+	loLo := 0.0
+	if IntervalFromEndpoints(yComplementCenter, y.Lo).Contains(i.Lo) {
+		loLo = positiveDistance(i.Lo, y.Lo)
+	}
+
+	return Angle(math.Max(hiHi, loLo))
+}
+
+// Project returns the closest point in the interval to the given point p.
+// The interval must be non-empty.
+func (i Interval) Project(p float64) float64 {
+	if p == -math.Pi {
+		p = math.Pi
+	}
+	if i.fastContains(p) {
+		return p
+	}
+	// Compute distance from p to each endpoint.
+	dlo := positiveDistance(p, i.Lo)
+	dhi := positiveDistance(i.Hi, p)
+	if dlo < dhi {
+		return i.Lo
+	}
+	return i.Hi
+}