<|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.

266. Berlion
time limit per test: 0.25 sec.
memory limit per test: 65536 KB
input: standard
output: standard



Far far away in the galaxy there is a planet called Berlion. The legend says that the Precursors left rich heritage on the planet. There are a lot of galactic races competing to claim the heritage. By the will of mighty Precursors each civilization has only one chance to land on the planet and explore small area of the surface. Powerful guards left by Precursors keep competitors from violation of the will. No race can compete the guards in a combat.
Recently race of bertonoids has built two super-scanners to scan the surface for artifacts. Bertonoids set up scanners in two locations of the space, where the guards hopefully won't detect them. You can assume that the planet has the form of a sphere and the scanners are points.
You, as intern in Forecast Department, have to find the area of the surface the scanners can explore. The scanners can explore only the directly visible part of the surface (you can assume that there are no any other objects in the sector with Berlion and scanners).
So your task is to find the area of that part of the planet surface, from which at least one of the scanners is visible.

Input
The first line of the input file contains the integer number R -- the radius of the planet (1 <= R <= 10). The second line contains the coordinates of the first scanner relatively to the center of the planet -- 3 numbers -- the distances along the x, y and z axes of an ortogonal normalized coordinate system. The third line contains the coordinates of the second scanner in the same format. All coordinates are integers and do not exceed 100 by the absolute value. Both scanners are located at the distance not less than 10-2 above the planet surface.

Output
Output the visible area of the planet surface. The answer should be rounded to 3 decimal digits.

Sample test(s)

Input
Test #1

1
0 0 100
0 100 0

Test #2

1
0 0 100
0 0 -100

Output
Test #1

9.362

Test #2

12.441
Author:	Ilya V. Elterman
Resource:	ACM ICPC 2004-2005, NEERC, Southern Subregional Contest
Date:	Saratov, October 7, 2004

<|response|>
1. Abridged Problem Statement

We are given a sphere of radius R centered at the origin and two scanner points outside the sphere.

Each scanner can see the part of the sphere surface that is directly visible from its position. We need to compute the total surface area visible from at least one of the two scanners.

Output the answer rounded to 3 decimal digits.

2. Key Observations

Observation 1: One scanner sees a spherical cap
Let a scanner be at position vector p. Let d = |p| and u = p / d, where u is the unit vector from the sphere center toward the scanner. A point x on the sphere is visible from the scanner if the segment from x to the scanner does not go through the sphere. The boundary of visibility is the tangent circle from the scanner to the sphere. On the unit sphere, the visible region is n · u >= R / d where n = x / R. So the scanner sees a spherical cap whose angular radius a satisfies cos(a) = R / d. The area of this cap on a sphere of radius R is 2πR²(1 - cos(a)) = 2πR²(1 - R / d).

Observation 2: We need the union of two spherical caps
Each scanner sees one cap. Let a = angular radius of cap 1, b = angular radius of cap 2, t = angular distance between the cap centers. The answer is area(cap1 ∪ cap2) = area(cap1) + area(cap2) - area(cap1 ∩ cap2). So the main challenge is computing the intersection area of two spherical caps.

Observation 3: There are three possible relationships between the caps
Case 1: Disjoint caps. If t >= a + b, the caps do not overlap, so intersection = 0.
Case 2: One cap fully contains the other. If t + min(a, b) <= max(a, b), then the smaller cap lies completely inside the larger cap, so intersection = area of smaller cap.
Case 3: Partial overlap. Otherwise, the intersection is a spherical lens, computed using spherical geometry.

3. Full Solution Approach

Let the scanner positions be p1 and p2. First compute their distances from the origin: d1 = |p1|, d2 = |p2|. The angular radii of the two visible caps are a = acos(R / d1), b = acos(R / d2). Let ca = cos(a) = R / d1, cb = cos(b) = R / d2; also sa = sin(a), sb = sin(b).

The angle t between the two cap centers is the angle between vectors p1 and p2: cos(t) = (p1 · p2) / (|p1| |p2|), then t = acos(cos(t)).

The individual cap areas are area1 = 2πR²(1 - ca), area2 = 2πR²(1 - cb).

Partial overlap formula
In the partial overlap case, the two caps form a spherical lens. Define ct = cos(t), st = sin(t). Let psi_a be half the angular span of the lens arc along cap 1's boundary (measured around cap 1's center), and similarly psi_b for cap 2. Using spherical trigonometry:
  cos(psi_a) = (cb - ca * ct) / (sa * st)
  cos(psi_b) = (ca - cb * ct) / (sb * st)
Define gamma as the angle between the two boundary tangent directions at an intersection point:
  cos(gamma) = (ct - ca * cb) / (sa * sb)
Then the lens area on the sphere of radius R is:
  intersection = R² * (2(π - gamma) - 2 * psi_a * ca - 2 * psi_b * cb)
Finally, answer = area1 + area2 - intersection.

Because floating-point calculations may produce values like 1.0000000002, every argument passed to acos should be clamped into [-1, 1]. The algorithm is constant time: O(1) time, O(1) memory.

4. C++ Implementation

```cpp
#include <bits/stdc++.h>

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 coord3_t = double;

const coord3_t PI = acos((coord3_t)-1.0);

struct Point3 {
    coord3_t x, y, z;

    coord3_t dot(const Point3& p) const { return x * p.x + y * p.y + z * p.z; }
    coord3_t norm() const { return sqrt(dot(*this)); }

    friend istream& operator>>(istream& in, Point3& p) {
        return in >> p.x >> p.y >> p.z;
    }
};

coord3_t R;
Point3 p1, p2;

void read() { cin >> R >> p1 >> p2; }

void solve() {
    // A scanner at distance d from the center sees a spherical cap. The
    // tangent points form a circle lying in the plane perpendicular to the
    // line center-scanner at distance R^2 / d from the center, so in terms
    // of the unit vector u toward the scanner the cap is
    //
    //     { x on sphere : (x / R) . u >= R / d },
    //
    // i.e. a cap of half-angle a = acos(R / d) and area 2 pi R^2 (1 - R/d).
    //
    // By inclusion-exclusion,
    //
    //     answer = area(cap_1) + area(cap_2) - area(cap_1 cap cap_2),
    //
    // so we only need the intersection. On the sphere, the intersection of
    // two caps is a "lens" bounded by two small-circle arcs meeting at the
    // two tangent-circle crossings. Let a, b be the cap half-angles and t
    // the angle between their axes u_a, u_b.
    //
    //   1) Each arc spans a half-angle psi at its pole. Restricting cap a's
    //      boundary (x . u_a = cos a) to also lie in cap b gives
    //
    //          cos(psi_a) = (cos b - cos a cos t) / (sin a sin t),
    //          cos(psi_b) = (cos a - cos b cos t) / (sin b sin t).
    //
    //   2) Each lens corner P is the apex of the spherical triangle
    //      (P, u_a, u_b), so the spherical law of cosines gives the angle g
    //      between the two boundary tangents at P:
    //
    //          cos(g) = (cos t - cos a cos b) / (sin a sin b).
    //
    //      The interior angle of the lens at P is then pi - g.
    //
    //   3) Apply Gauss-Bonnet on the unit sphere (K = 1). A small circle of
    //      half-angle r has geodesic curvature cot r, and an arc spanning
    //      2 psi at its pole has length 2 psi sin r, so the integral of
    //      k_g ds over each arc is 2 psi cos r. With two corners of
    //      interior angle pi - g (i.e. exterior angle g) the theorem reads
    //
    //          area + 2 psi_a cos a + 2 psi_b cos b + 2 g = 2 pi,
    //
    //      hence on a sphere of radius R
    //
    //          area(lens) = R^2 (2(pi - g) - 2 psi_a cos a - 2 psi_b cos b).
    //
    // The disjoint / nested cases are handled explicitly up front; otherwise
    // we clamp the cosines into [-1, 1] before acos for numerical safety.

    coord3_t d1 = p1.norm(), d2 = p2.norm();
    coord3_t cos_a = R / d1, cos_b = R / d2;
    coord3_t sin_a = sqrt(max(0.0, 1.0 - cos_a * cos_a));
    coord3_t sin_b = sqrt(max(0.0, 1.0 - cos_b * cos_b));
    coord3_t alpha = acos(cos_a), beta = acos(cos_b);

    coord3_t cos_t = p1.dot(p2) / (d1 * d2);
    cos_t = max((coord3_t)-1.0, min((coord3_t)1.0, cos_t));
    coord3_t sin_t = sqrt(max(0.0, 1.0 - cos_t * cos_t));
    coord3_t theta = acos(cos_t);

    coord3_t area1 = 2 * PI * R * R * (1 - cos_a);
    coord3_t area2 = 2 * PI * R * R * (1 - cos_b);

    coord3_t intersection = 0.0;
    if(theta >= alpha + beta) {
        intersection = 0.0;
    } else if(theta + min(alpha, beta) <= max(alpha, beta)) {
        intersection = 2 * PI * R * R * (1 - (alpha < beta ? cos_a : cos_b));
    } else {
        auto clamp_cos = [](coord3_t c) {
            return max((coord3_t)-1.0, min((coord3_t)1.0, c));
        };
        coord3_t psi_a =
            acos(clamp_cos((cos_b - cos_a * cos_t) / (sin_a * sin_t)));
        coord3_t psi_b =
            acos(clamp_cos((cos_a - cos_b * cos_t) / (sin_b * sin_t)));
        coord3_t gamma =
            acos(clamp_cos((cos_t - cos_a * cos_b) / (sin_a * sin_b)));
        intersection =
            R * R * (2 * (PI - gamma) - 2 * psi_a * cos_a - 2 * psi_b * cos_b);
    }

    cout << fixed << setprecision(3) << area1 + area2 - intersection << '\n';
}

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

```python
import math
import sys


EPS = 1e-12


def clamp_cos(x):
    """
    Clamp x into [-1, 1].

    This is necessary before acos because floating-point arithmetic
    may produce values slightly outside the valid range.
    """
    return max(-1.0, min(1.0, x))


def dot(a, b):
    """
    Dot product of two 3D vectors.
    """
    return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]


def norm(a):
    """
    Euclidean length of a 3D vector.
    """
    return math.sqrt(dot(a, a))


def main():
    data = sys.stdin.read().strip().split()
    if not data:
        return

    R = float(data[0])

    p1 = (
        float(data[1]),
        float(data[2]),
        float(data[3]),
    )

    p2 = (
        float(data[4]),
        float(data[5]),
        float(data[6]),
    )

    pi = math.pi

    # Distances from sphere center to scanners.
    d1 = norm(p1)
    d2 = norm(p2)

    # Angular radii of visible caps satisfy:
    # cos(a) = R / d1
    # cos(b) = R / d2
    ca = clamp_cos(R / d1)
    cb = clamp_cos(R / d2)

    a = math.acos(ca)
    b = math.acos(cb)

    sa = math.sqrt(max(0.0, 1.0 - ca * ca))
    sb = math.sqrt(max(0.0, 1.0 - cb * cb))

    # Angle between the two scanner directions.
    ct = dot(p1, p2) / (d1 * d2)
    ct = clamp_cos(ct)

    t = math.acos(ct)
    st = math.sqrt(max(0.0, 1.0 - ct * ct))

    # Area of each spherical cap.
    area1 = 2.0 * pi * R * R * (1.0 - ca)
    area2 = 2.0 * pi * R * R * (1.0 - cb)

    # Compute intersection area of the two caps.
    if t >= a + b - EPS:
        # Disjoint or tangent caps.
        intersection = 0.0

    elif t + min(a, b) <= max(a, b) + EPS:
        # One cap completely contains the other.
        intersection = min(area1, area2)

    else:
        # Partial overlap: spherical lens.

        # Half-span of the arc from cap 1 boundary.
        psi_a = math.acos(
            clamp_cos((cb - ca * ct) / (sa * st))
        )

        # Half-span of the arc from cap 2 boundary.
        psi_b = math.acos(
            clamp_cos((ca - cb * ct) / (sb * st))
        )

        # Angle between boundary tangent directions.
        gamma = math.acos(
            clamp_cos((ct - ca * cb) / (sa * sb))
        )

        # Lens area formula, scaled by R^2.
        intersection = R * R * (
            2.0 * (pi - gamma)
            - 2.0 * psi_a * ca
            - 2.0 * psi_b * cb
        )

        # Numerical safety.
        if -1e-9 < intersection < 0.0:
            intersection = 0.0

    # Union area by inclusion-exclusion.
    answer = area1 + area2 - intersection

    print(f"{answer:.3f}")


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