<|instruction|>
Solve the below problem. The solution should start with an abridged problem statement. Then key observations. Then full solution based on the observations. Then C++ and Python implementations with comments.

351. A Mission for a Scout
Time limit per test: 0.25 second(s)
Memory limit: 65536 kilobytes
input: standard
output: standard



The scout Sidorov is facing a delicate situation: he has to get from A to B. The embarrassment of the situation apparently lies in the news that the scout received lately. It appears that a bomb has been detonated at the neighborhood recently, and the explosion caused a circle crater. Fortunately, Sidorov has an antigravitator to fly over the craters, but the fuel in the tank is hardly enough for the distance D. Sidorov wrote out all the necessary data he had and sent it to you by secret post. He is looking forward to hearing from you, for nobody else can calculate the length of the shortest route from A to B, taking into consideration the craters and the fuel left. He also would like to get the route itself. Clearly, the route will consist only of segments and arcs of the circle formed by the crater.

Input
The first line describes four numbers XA, YA, XB, YB — coordinates of A and B accordingly. The second line describes three numbers XC, YC, RC — coordinates of the center and the radius of the crater. The third line describes D — the distance the antigravitator can fly over the crater. All coordinates in the input file are integers from 0 to 10\,000. RC and D are integer numbers from 1 to 10\,000. A doesn't coincide with B and both points are outside the crater.

Output
Write the number of the parts T of the desired route to the first line of the output file. After that write the description of route parts one per line. For each part of the route print character "S" (for segment) or "A" (for arc) followed by the coordinates of their endpoints. Coordinates should be written with 8 digits after decimal point. T should not exceed 100. The first point in the route should be point A, the last point — B, and the second endpoint of each part should coincide with the first endpoint of the next part in the route.
Example(s)
sample input
sample output
1 1 101 1
50 1 40
80
1
S 1.000000000 1.000000000 101.000000000 1.000000000

sample input
sample output
1 1 101 1
50 1 40
1
4
S 1.00000000 1.00000000 17.34693878 24.10362726
S 17.34693878 24.10362726 17.93468841 24.91267013
A 17.93468841 24.91267013 81.37254902 25.81457572
S 81.37254902 25.81457572 101.00000000 1.00000000

<|response|>
## 1. Abridged problem statement

Given two points `A` and `B` outside a circular crater with center `C` and radius `R`, find and output a shortest route from `A` to `B`.

A route may consist of:

- straight segments;
- arcs of the crater boundary circle.

The route is allowed to pass through the crater, but the total length of all parts lying inside the circle must not exceed `D`.

Output any shortest valid route using at most `100` parts.

---

## 2. Key observations needed to solve the problem

### Observation 1 — Direct segment may already be optimal

The shortest possible Euclidean route is the straight segment `AB`.

If `AB` does not enter the crater, it is valid and optimal.

If `AB` enters the crater, then the part inside the circle is a chord. If that chord length is at most `D`, the straight segment is still valid and optimal.

---

### Observation 2 — Otherwise, optimal routes touch the circle boundary

If direct travel is impossible, the path must avoid part of the crater. Around a circle, shortest outside paths are made of:

```text
segment from A to tangent point
arc along circle
segment from tangent point to B
```

There are two such routes: clockwise and counterclockwise.

---

### Observation 3 — Fuel can replace boundary arc by one chord

Passing through the crater between two boundary points means using a chord.

A chord of length `x` in a circle of radius `R` subtends central angle:

```text
theta = 2 * asin(x / (2R))
```

Replacing an arc by a chord saves length:

```text
R * theta - x
```

It is never worse to spend all usable fuel on one chord rather than splitting it into several chords.

So for a fixed side of the circle, we can try replacing part of the circular arc by one chord of length `D`.

---

### Observation 4 — Two possible chord situations

For one side of the circle, suppose the tangent-to-tangent arc angle is `Phi`.

Let:

```text
thetaD = 2 * asin(D / (2R))
```

There are two cases.

#### Case A: `thetaD <= Phi`

The chord fits inside the tangent arc.

Then the optimal path for this side may be:

```text
A -> tangent from A -> chord of length D -> remaining arc -> tangent to B -> B
```

The chord can be placed anywhere on the arc, because all chords of length `D` save the same amount. We place it immediately after the tangent point from `A`.

#### Case B: `thetaD > Phi`

The chord is too large to fit between the two tangent points.

Then the optimal path may have no boundary arc:

```text
A -> P -> Q -> B
```

where:

- `P` and `Q` lie on the circle;
- `PQ` is a chord of length `D`;
- the angle between `P` and `Q` is `thetaD`;
- segments `AP` and `QB` must not enter the crater.

Now only one variable remains: the angle of `P`.

We minimize:

```text
|A - P| + D + |Q - B|
```

over all feasible positions of `P`.

This is a one-dimensional unimodal function, so ternary search is enough.

---

## 3. Full solution approach

### Step 1 — Basic geometry utilities

We need:

- vector addition/subtraction;
- dot product;
- cross product;
- distance;
- angle of a vector;
- point on the circle by angle.

For angle `alpha`, the point on the crater boundary is:

```text
C + R * (cos(alpha), sin(alpha))
```

---

### Step 2 — Check whether a segment enters the crater

For segment `P0P1`, find the point on the segment closest to the center `C`.

Let:

```text
v = P1 - P0
t = dot(C - P0, v) / |v|²
```

Clamp `t` to `[0, 1]`.

The closest point is:

```text
P0 + t * v
```

If its distance from `C` is at least `R`, the segment does not enter the crater.

---

### Step 3 — Try direct route

If `AB` does not enter the circle, output:

```text
S A B
```

Otherwise, compute the distance `h` from `C` to the line `AB`.

The chord cut by the circle has length:

```text
2 * sqrt(R² - h²)
```

If this is at most `D`, output the direct segment.

---

### Step 4 — Compute tangent information

For an outside point `P`:

```text
d = |P - C|
psi = angle(C -> P)
gamma = acos(R / d)
```

The two tangent points have angular coordinates:

```text
psi + gamma
psi - gamma
```

For each side `s = +1` or `s = -1`, choose compatible tangent angles:

```text
from A: a1 = psiA + s * gammaA
to B:   a2 = psiB - s * gammaB
```

The arc angle from `a1` to `a2` in direction `s` is:

```text
Phi = s * (a2 - a1) mod 2π
```

---

### Step 5 — Consider candidates for both sides

For each side:

#### Candidate 1 — Pure outside detour

```text
A -> tangent A -> arc -> tangent B -> B
```

Length:

```text
tangentA + R * Phi + tangentB
```

---

#### Candidate 2 — Tangent path with one chord

If `thetaD <= Phi`, use:

```text
A -> tangent A -> chord -> remaining arc -> tangent B -> B
```

Length:

```text
tangentA + D + R * (Phi - thetaD) + tangentB
```

---

#### Candidate 3 — Free chord path

Try:

```text
A -> P -> Q -> B
```

where:

```text
angle(Q) = angle(P) + s * thetaD
```

The point `P` must be visible from `A`, and `Q` must be visible from `B`.

Visibility from an outside point is exactly the interval between its two tangent points.

So we intersect:

```text
P angle in A-visible interval
Q angle in B-visible interval
```

Then minimize:

```text
|A - P| + D + |Q - B|
```

with ternary search.

---

### Step 6 — Output best candidate

Keep the shortest valid route among all candidates and print its parts.

The number of parts is at most `4`.

Complexity:

```text
O(1)
```

---
## 4. C++ implementation

```cpp
#include <bits/stdc++.h>
// #include <coding_library/geometry/point.hpp>

using namespace std;

template<typename T1, typename T2>
ostream& operator<<(ostream& out, const pair<T1, T2>& x) {
    return out << x.first << ' ' << x.second;
}

template<typename T1, typename T2>
istream& operator>>(istream& in, pair<T1, T2>& x) {
    return in >> x.first >> x.second;
}

template<typename T>
istream& operator>>(istream& in, vector<T>& a) {
    for(auto& x: a) {
        in >> x;
    }
    return in;
};

template<typename T>
ostream& operator<<(ostream& out, const vector<T>& a) {
    for(auto x: a) {
        out << x << ' ';
    }
    return out;
};

using coord_t = double;

struct Point {
    static constexpr coord_t eps = 1e-9;
    static inline const coord_t PI = acos((coord_t)-1.0);

    coord_t x, y;
    Point(coord_t x = 0, coord_t y = 0) : x(x), y(y) {}

    Point operator+(const Point& p) const { return Point(x + p.x, y + p.y); }
    Point operator-(const Point& p) const { return Point(x - p.x, y - p.y); }
    Point operator*(coord_t c) const { return Point(x * c, y * c); }
    Point operator/(coord_t c) const { return Point(x / c, y / c); }

    coord_t operator*(const Point& p) const { return x * p.x + y * p.y; }
    coord_t operator^(const Point& p) const { return x * p.y - y * p.x; }

    bool operator==(const Point& p) const { return x == p.x && y == p.y; }
    bool operator!=(const Point& p) const { return x != p.x || y != p.y; }
    bool operator<(const Point& p) const {
        return x != p.x ? x < p.x : y < p.y;
    }
    bool operator>(const Point& p) const {
        return x != p.x ? x > p.x : y > p.y;
    }
    bool operator<=(const Point& p) const {
        return x != p.x ? x < p.x : y <= p.y;
    }
    bool operator>=(const Point& p) const {
        return x != p.x ? x > p.x : y >= p.y;
    }

    coord_t norm2() const { return x * x + y * y; }
    coord_t norm() const { return sqrt(norm2()); }
    coord_t angle() const { return atan2(y, x); }

    Point rotate(coord_t a) const {
        return Point(x * cos(a) - y * sin(a), x * sin(a) + y * cos(a));
    }

    Point perp() const { return Point(-y, x); }
    Point unit() const { return *this / norm(); }
    Point normal() const { return perp().unit(); }
    Point project(const Point& p) const {
        return *this * (*this * p) / norm2();
    }
    Point reflect(const Point& p) const {
        return *this * 2 * (*this * p) / norm2() - p;
    }

    friend ostream& operator<<(ostream& os, const Point& p) {
        return os << p.x << ' ' << p.y;
    }
    friend istream& operator>>(istream& is, Point& p) {
        return is >> p.x >> p.y;
    }

    friend int ccw(const Point& a, const Point& b, const Point& c) {
        coord_t v = (b - a) ^ (c - a);
        if(-eps <= v && v <= eps) {
            return 0;
        } else if(v > 0) {
            return 1;
        } else {
            return -1;
        }
    }

    friend bool point_on_segment(
        const Point& a, const Point& b, const Point& p
    ) {
        return ccw(a, b, p) == 0 && p.x >= min(a.x, b.x) - eps &&
               p.x <= max(a.x, b.x) + eps && p.y >= min(a.y, b.y) - eps &&
               p.y <= max(a.y, b.y) + eps;
    }

    friend bool point_in_triangle(
        const Point& a, const Point& b, const Point& c, const Point& p
    ) {
        int d1 = ccw(a, b, p);
        int d2 = ccw(b, c, p);
        int d3 = ccw(c, a, p);
        return (d1 >= 0 && d2 >= 0 && d3 >= 0) ||
               (d1 <= 0 && d2 <= 0 && d3 <= 0);
    }

    friend Point line_line_intersection(
        const Point& a1, const Point& b1, const Point& a2, const Point& b2
    ) {
        return a1 +
               (b1 - a1) * ((a2 - a1) ^ (b2 - a2)) / ((b1 - a1) ^ (b2 - a2));
    }

    friend bool collinear(const Point& a, const Point& b) {
        return abs(a ^ b) < eps;
    }

    friend Point circumcenter(const Point& a, const Point& b, const Point& c) {
        Point mid_ab = (a + b) / 2.0;
        Point mid_ac = (a + c) / 2.0;
        Point perp_ab = (b - a).perp();
        Point perp_ac = (c - a).perp();
        return line_line_intersection(
            mid_ab, mid_ab + perp_ab, mid_ac, mid_ac + perp_ac
        );
    }

    friend coord_t arc_area(
        const Point& center, coord_t r, const Point& p1, const Point& p2
    ) {
        coord_t theta1 = (p1 - center).angle();
        coord_t theta2 = (p2 - center).angle();
        if(theta2 < theta1 - eps) {
            theta2 += 2 * PI;
        }

        coord_t d_theta = theta2 - theta1;
        coord_t cx = center.x, cy = center.y;
        coord_t area = r * cx * (sin(theta2) - sin(theta1)) -
                       r * cy * (cos(theta2) - cos(theta1)) + r * r * d_theta;
        return area / 2.0;
    }

    friend vector<Point> intersect_circles(
        const Point& c1, coord_t r1, const Point& c2, coord_t r2
    ) {
        Point d = c2 - c1;
        coord_t dist = d.norm();

        if(dist > r1 + r2 + eps || dist < abs(r1 - r2) - eps || dist < eps) {
            return {};
        }

        coord_t a = (r1 * r1 - r2 * r2 + dist * dist) / (2 * dist);
        coord_t h_sq = r1 * r1 - a * a;
        if(h_sq < -eps) {
            return {};
        }
        if(h_sq < 0) {
            h_sq = 0;
        }
        coord_t h = sqrt(h_sq);

        Point mid = c1 + d.unit() * a;
        Point perp_dir = d.perp().unit();

        if(h < eps) {
            return {mid};
        }
        return {mid + perp_dir * h, mid - perp_dir * h};
    }

    friend optional<Point> intersect_ray_segment(
        const Point& ray_start, const Point& ray_through, const Point& seg_a,
        const Point& seg_b
    ) {
        Point ray_dir = ray_through - ray_start;
        if(ray_dir.norm2() < Point::eps) {
            return {};
        }
        Point seg_dir = seg_b - seg_a;
        coord_t denom = ray_dir ^ seg_dir;
        if(fabs(denom) < eps) {
            return {};
        }
        coord_t t = ((seg_a - ray_start) ^ seg_dir) / denom;
        if(t < eps) {
            return {};
        }
        coord_t s = ((seg_a - ray_start) ^ ray_dir) / denom;
        if(s < eps || s > 1 - eps) {
            return {};
        }
        return ray_start + ray_dir * t;
    }
};

Point A, B, center;
coord_t R, D;

void read() { cin >> A >> B >> center >> R >> D; }

struct Part {
    char type;
    Point a, b;
};

void solve() {
    // We travel from A to B; the disc may be entered, but the total length of
    // route lying inside it must not exceed the fuel budget D. If the straight
    // segment A-B stays outside the disc, or the chord it cuts through the disc
    // has length at most D, then that single segment is already the answer.
    //
    // Otherwise the optimum is a taut path that detours around the crater on
    // one of the two sides. Around a convex obstacle a taut path leaves A along
    // a tangent, hugs the boundary, and reaches B along a tangent. The fuel
    // lets us replace one sub-arc of the boundary by a straight chord. A chord
    // of length D always subtends the same central angle 2*asin(D/2R), so the
    // saving (arc length minus the chord) is identical wherever the chord sits,
    // and one chord is never worse than splitting the budget into several. Thus
    // while the cut-out arc still fits inside the tangent-to-tangent arc Phi,
    // the route is tangent, chord of length D, remaining arc, tangent, and we
    // are free to place the chord right after A's tangent point.
    //
    // When D is large enough that the chord would subtend more than Phi, the
    // boundary no longer pins the touch points: the path degenerates to
    // A -> P -> Q -> B where P and Q sit on the circle exactly the chord angle
    // apart and the two outer segments must not re-enter the disc. There we
    // minimise |A-P| + D + |Q-B| over the single free angle by a dense scan
    // refined with ternary search, treating angles whose outer segments clip
    // the disc as infeasible.
    //
    // We build all three shapes (pure detour, tangent+chord+arc, and the free
    // two-tangent chord) for both sides and keep the shortest. Each emitted arc
    // spans less than 180 degrees, so naming only its endpoints is unambiguous.

    const coord_t INF = 1e18;

    auto on_circle = [&](coord_t ang) {
        return center + Point(cos(ang), sin(ang)) * R;
    };

    auto seg_ok = [&](const Point& p0, const Point& p1) {
        Point d = p1 - p0;
        coord_t dn = d.norm2();
        coord_t t = 0;
        if(dn > Point::eps) {
            t = ((center - p0) * d) / dn;
            t = max((coord_t)0, min((coord_t)1, t));
        }
        Point cl = p0 + d * t;
        return (cl - center).norm() >= R - 1e-7;
    };

    auto print = [&](const vector<Part>& parts) {
        vector<Part> p;
        for(const auto& q: parts) {
            if((q.a - q.b).norm() > 1e-9) {
                p.push_back(q);
            }
        }

        cout << fixed << setprecision(8);
        cout << p.size() << '\n';
        for(const auto& q: p) {
            cout << q.type << ' ' << q.a.x << ' ' << q.a.y << ' ' << q.b.x
                 << ' ' << q.b.y << '\n';
        }
    };

    if(seg_ok(A, B)) {
        print({{'S', A, B}});
        return;
    }

    coord_t h = fabs((B - A) ^ (center - A)) / (B - A).norm();
    coord_t chord_ab = 2 * sqrt(max((coord_t)0, R * R - h * h));
    if(chord_ab <= D + 1e-9) {
        print({{'S', A, B}});
        return;
    }

    coord_t dA = (A - center).norm(), dB = (B - center).norm();
    coord_t tA = sqrt(dA * dA - R * R), tB = sqrt(dB * dB - R * R);
    coord_t psiA = (A - center).angle(), psiB = (B - center).angle();
    coord_t gA = acos(R / dA), gB = acos(R / dB);
    coord_t theta_d = 2 * asin(min((coord_t)1, D / (2 * R)));

    coord_t best_len = INF;
    vector<Part> best;
    auto consider = [&](coord_t len, const vector<Part>& parts) {
        if(len < best_len) {
            best_len = len;
            best = parts;
        }
    };

    for(int s: {1, -1}) {
        coord_t a1 = psiA + s * gA;
        coord_t a2 = psiB - s * gB;
        coord_t Phi = fmod(s * (a2 - a1), 2 * Point::PI);
        if(Phi < 0) {
            Phi += 2 * Point::PI;
        }

        Point Ta = on_circle(a1), Tb = on_circle(a2);

        consider(
            tA + tB + R * Phi, {{'S', A, Ta}, {'A', Ta, Tb}, {'S', Tb, B}}
        );

        if(theta_d <= Phi + 1e-12) {
            Point P2 = on_circle(a1 + s * theta_d);
            consider(
                tA + tB + R * (Phi - theta_d) + D,
                {{'S', A, Ta}, {'S', Ta, P2}, {'A', P2, Tb}, {'S', Tb, B}}
            );
        }

        auto g = [&](coord_t al) -> coord_t {
            Point P = on_circle(al), Q = on_circle(al + s * theta_d);
            if(!seg_ok(A, P) || !seg_ok(Q, B)) {
                return INF;
            }
            return (A - P).norm() + D + (Q - B).norm();
        };

        coord_t cB = psiB - s * theta_d;
        coord_t d = atan2(sin(cB - psiA), cos(cB - psiA));
        coord_t lo = max(-gA, d - gB), hi = min(gA, d + gB);
        if(lo <= hi) {
            lo += psiA;
            hi += psiA;
            for(int it = 0; it < 200; it++) {
                coord_t m1 = lo + (hi - lo) / 3, m2 = hi - (hi - lo) / 3;
                if(g(m1) < g(m2)) {
                    hi = m2;
                } else {
                    lo = m1;
                }
            }

            coord_t al = (lo + hi) / 2, c = g(al);
            if(c < INF / 2) {
                Point P = on_circle(al), Q = on_circle(al + s * theta_d);
                consider(c, {{'S', A, P}, {'S', P, Q}, {'S', Q, B}});
            }
        }
    }

    print(best);
}

int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(nullptr);

    int T = 1;
    // cin >> T;
    for(int test = 1; test <= T; test++) {
        read();
        solve();
    }

    return 0;
}
```

---

## 5. Python implementation with detailed comments

```python
import sys
import math


EPS = 1e-9
PI = math.pi
INF = 10 ** 100


class Point:
    """Simple 2D point/vector class."""

    def __init__(self, x=0.0, y=0.0):
        self.x = float(x)
        self.y = float(y)

    def __add__(self, other):
        return Point(self.x + other.x, self.y + other.y)

    def __sub__(self, other):
        return Point(self.x - other.x, self.y - other.y)

    def __mul__(self, value):
        """
        If value is another Point, return dot product.
        Otherwise, multiply by scalar.
        """
        if isinstance(value, Point):
            return self.x * value.x + self.y * value.y
        return Point(self.x * value, self.y * value)

    def __rmul__(self, value):
        return self.__mul__(value)

    def __truediv__(self, value):
        return Point(self.x / value, self.y / value)

    def cross(self, other):
        return self.x * other.y - self.y * other.x

    def norm2(self):
        return self.x * self.x + self.y * self.y

    def norm(self):
        return math.sqrt(self.norm2())

    def angle(self):
        return math.atan2(self.y, self.x)


def dist(a, b):
    return (a - b).norm()


def main():
    data = list(map(float, sys.stdin.read().split()))

    xa, ya, xb, yb = data[0:4]
    xc, yc, r = data[4:7]
    D = data[7]

    A = Point(xa, ya)
    B = Point(xb, yb)
    C = Point(xc, yc)
    R = r

    def on_circle(angle):
        """
        Point on the crater boundary with polar angle `angle`
        around center C.
        """
        return C + Point(math.cos(angle), math.sin(angle)) * R

    def segment_outside_circle(p0, p1):
        """
        Return True if segment p0-p1 does not enter the open disk.
        Touching the boundary is allowed.
        """
        v = p1 - p0
        length2 = v.norm2()

        t = 0.0

        if length2 > EPS:
            t = ((C - p0) * v) / length2
            t = max(0.0, min(1.0, t))

        closest = p0 + v * t

        return dist(closest, C) >= R - 1e-7

    def print_route(parts):
        """
        Print route parts, removing zero-length parts.

        Each part is:
            (type, start_point, end_point)
        """
        filtered = []

        for typ, p, q in parts:
            if dist(p, q) > 1e-9:
                filtered.append((typ, p, q))

        print(len(filtered))

        for typ, p, q in filtered:
            print(
                f"{typ} "
                f"{p.x:.8f} {p.y:.8f} "
                f"{q.x:.8f} {q.y:.8f}"
            )

    # ------------------------------------------------------------
    # Step 1: Try direct route A -> B.
    # ------------------------------------------------------------

    if segment_outside_circle(A, B):
        print_route([("S", A, B)])
        return

    # Distance from C to line AB.
    line_distance = abs((B - A).cross(C - A)) / dist(A, B)

    # Chord length inside the crater.
    chord_ab = 2.0 * math.sqrt(
        max(0.0, R * R - line_distance * line_distance)
    )

    if chord_ab <= D + 1e-9:
        print_route([("S", A, B)])
        return

    # ------------------------------------------------------------
    # Step 2: Compute tangent-related values.
    # ------------------------------------------------------------

    dA = dist(A, C)
    dB = dist(B, C)

    tangentA = math.sqrt(dA * dA - R * R)
    tangentB = math.sqrt(dB * dB - R * R)

    psiA = (A - C).angle()
    psiB = (B - C).angle()

    gammaA = math.acos(R / dA)
    gammaB = math.acos(R / dB)

    # Chord of length D subtends this central angle.
    thetaD = 2.0 * math.asin(min(1.0, D / (2.0 * R)))

    best_len = INF
    best_route = []

    def consider(length, route):
        nonlocal best_len, best_route

        if length < best_len:
            best_len = length
            best_route = route

    # ------------------------------------------------------------
    # Step 3: Try both directions around the circle.
    # ------------------------------------------------------------

    for s in (1, -1):
        # Tangent point from A on this side.
        angleA = psiA + s * gammaA

        # Tangent point to B on this side.
        angleB = psiB - s * gammaB

        # Angular length of arc from angleA to angleB in direction s.
        Phi = math.fmod(s * (angleB - angleA), 2.0 * PI)

        if Phi < 0.0:
            Phi += 2.0 * PI

        TA = on_circle(angleA)
        TB = on_circle(angleB)

        # --------------------------------------------------------
        # Candidate 1: pure outside detour.
        # A -> tangent -> arc -> tangent -> B
        # --------------------------------------------------------

        consider(
            tangentA + R * Phi + tangentB,
            [
                ("S", A, TA),
                ("A", TA, TB),
                ("S", TB, B),
            ],
        )

        # --------------------------------------------------------
        # Candidate 2: replace part of the arc by one chord.
        # --------------------------------------------------------

        if thetaD <= Phi + 1e-12:
            after_chord = on_circle(angleA + s * thetaD)

            consider(
                tangentA + D + R * (Phi - thetaD) + tangentB,
                [
                    ("S", A, TA),
                    ("S", TA, after_chord),
                    ("A", after_chord, TB),
                    ("S", TB, B),
                ],
            )

        # --------------------------------------------------------
        # Candidate 3:
        # A -> P -> Q -> B
        #
        # P and Q are on the boundary.
        # PQ is a chord of length D.
        # angle(Q) = angle(P) + s * thetaD.
        # --------------------------------------------------------

        def free_chord_cost(alpha):
            P = on_circle(alpha)
            Q = on_circle(alpha + s * thetaD)

            if not segment_outside_circle(A, P):
                return INF

            if not segment_outside_circle(Q, B):
                return INF

            return dist(A, P) + D + dist(Q, B)

        """
        Feasible angle interval for P.

        P must be visible from A:
            alpha in [psiA - gammaA, psiA + gammaA]

        Q must be visible from B:
            alpha + s * thetaD in [psiB - gammaB, psiB + gammaB]

        Therefore alpha must also lie in:
            [psiB - s * thetaD - gammaB,
             psiB - s * thetaD + gammaB]
        """

        shifted_center_b = psiB - s * thetaD

        # Bring shifted_center_b close to psiA to avoid 2*pi wrap issues.
        delta = math.atan2(
            math.sin(shifted_center_b - psiA),
            math.cos(shifted_center_b - psiA),
        )

        lo = max(-gammaA, delta - gammaB)
        hi = min(gammaA, delta + gammaB)

        if lo <= hi:
            lo += psiA
            hi += psiA

            # Ternary search over the feasible interval.
            for _ in range(200):
                m1 = lo + (hi - lo) / 3.0
                m2 = hi - (hi - lo) / 3.0

                if free_chord_cost(m1) < free_chord_cost(m2):
                    hi = m2
                else:
                    lo = m1

            alpha = (lo + hi) / 2.0
            cost = free_chord_cost(alpha)

            if cost < INF / 2:
                P = on_circle(alpha)
                Q = on_circle(alpha + s * thetaD)

                consider(
                    cost,
                    [
                        ("S", A, P),
                        ("S", P, Q),
                        ("S", Q, B),
                    ],
                )

    print_route(best_route)


if __name__ == "__main__":
    main()
```