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

217. Two Cylinders
time limit per test: 0.5 sec.
memory limit per test: 65536 KB
input: standard
output: standard



In this problem your task is very simple.

Consider two infinite cylinders in three-dimensional space, of radii R1 and R2 respectively, located in such a way that their axes intersect and are perpendicular.

Your task is to find the volume of their intersection.

Input

Input file contains two real numbers R1 and R2 (1 ≤ R1, R2 ≤ 100);

Output

Output the volume of the intersection of the cylinders. Your answer must be accurate up to 10-4.

Sample test(s)

Input
1 1

Output
5.3333

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

Two infinite right circular cylinders in 3D have radii \(R_1\) and \(R_2\). Their axes intersect at a point and are perpendicular. Compute the volume of the region that lies inside **both** cylinders.  
Input: \(R_1, R_2\) (real, 1…100).  
Output: intersection volume with error \(\le 10^{-4}\).

---

## 2) Key observations needed to solve the problem

1. **Convenient coordinates**: Place the origin at the intersection point of the axes, and align axes with coordinate axes.
2. **Cylinder equations**:
   - Cylinder along the \(x\)-axis, radius \(r\): \(y^2 + z^2 \le r^2\)
   - Cylinder along the \(y\)-axis, radius \(R\): \(x^2 + z^2 \le R^2\)
   where \(r=\min(R_1,R_2)\), \(R=\max(R_1,R_2)\).
3. **Slice by planes \(z=\text{const}\)**: For a fixed \(z\), the constraints become independent bounds on \(x\) and \(y\), producing a rectangle in the \(xy\)-plane.
4. **Symmetry**: The solid is symmetric in \(z\) (and also in signs of \(x,y\)), so we can integrate only over \(z\in[0,r]\) and multiply by 8.
5. **Numerical integration is enough**: The integral is smooth on \([0,r]\); Simpson’s rule with a large even number of intervals easily meets \(10^{-4}\) accuracy for a single test case.

---

## 3) Full solution approach

### Step A: Set up equations
Let:
- \(r=\min(R_1,R_2)\)
- \(R=\max(R_1,R_2)\)

Model the cylinders as:
\[
\text{Cyl}_1: y^2 + z^2 \le r^2,\qquad
\text{Cyl}_2: x^2 + z^2 \le R^2
\]

### Step B: Cross-sectional area at fixed \(z\)
Fix \(z\). Then:
- From \(\text{Cyl}_1\): \(y^2 \le r^2 - z^2 \Rightarrow |y|\le \sqrt{r^2-z^2}\)
- From \(\text{Cyl}_2\): \(x^2 \le R^2 - z^2 \Rightarrow |x|\le \sqrt{R^2-z^2}\)

So the intersection at height \(z\) is a rectangle:
- width in \(x\): \(2\sqrt{R^2-z^2}\)
- width in \(y\): \(2\sqrt{r^2-z^2}\)

Area:
\[
A(z)=4\sqrt{(r^2-z^2)(R^2-z^2)}
\]

Valid \(z\) range for being inside the smaller cylinder: \(|z|\le r\).

### Step C: Volume integral
\[
V=\int_{-r}^{r} A(z)\,dz
= \int_{-r}^{r} 4\sqrt{(r^2-z^2)(R^2-z^2)}\,dz
\]
The integrand is even, so:
\[
V = 8\int_{0}^{r}\sqrt{(r^2-z^2)(R^2-z^2)}\,dz
\]

Define:
\[
f(z)=\sqrt{(r^2-z^2)(R^2-z^2)}
\]
Then:
\[
V = 8\int_{0}^{r} f(z)\,dz
\]

### Step D: Simpson’s rule on \([0,r]\)
Choose an even \(n\) (e.g. \(10^6\)), step \(h=r/n\). Simpson’s approximation:
\[
\int_0^r f(z)\,dz \approx \frac{h}{3}\left(
f(0)+f(r)+4\sum_{i\ \text{odd}} f(ih)+2\sum_{i\ \text{even}, i\neq 0,n} f(ih)
\right)
\]
Finally:
\[
V \approx 8 \cdot \text{SimpsonIntegral}
\]

Print with 4 digits after decimal; this matches the required \(10^{-4}\).

---

## 4) C++ implementation (detailed comments)

```cpp
#include <bits/stdc++.h>
using namespace std;

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

    double R1, R2;
    cin >> R1 >> R2;

    // Let r be the smaller radius, R the larger radius.
    double r = min(R1, R2);
    double R = max(R1, R2);

    // Integrand for the reduced integral:
    // f(z) = sqrt((r^2 - z^2) * (R^2 - z^2)), valid for z in [0, r].
    auto f = [&](double z) -> double {
        // Slightly defensive: due to floating errors, clamp to 0 if needed.
        double a = r * r - z * z;
        double b = R * R - z * z;
        if (a < 0) a = 0;
        if (b < 0) b = 0;
        return sqrt(a * b);
    };

    // Simpson's rule parameters.
    // n must be even. 1e6 is fine for one test case and required accuracy.
    const int n = 1'000'000;
    const double h = r / n;

    // Simpson sum: f(0) + f(r) + 4*f(h) + 2*f(2h) + 4*f(3h) + ...
    double sum = f(0.0) + f(r);

    for (int i = 1; i < n; i++) {
        double z = i * h;
        sum += (i % 2 ? 4.0 : 2.0) * f(z);
    }

    // Approximate integral_0^r f(z) dz
    double integral = (h / 3.0) * sum;

    // Multiply by 8 due to symmetry:
    // V = 8 * integral_0^r f(z) dz
    double volume = 8.0 * integral;

    cout << fixed << setprecision(4) << volume << "\n";
    return 0;
}
```

---

## 5) Python implementation (detailed comments)

```python
import sys
import math

def main():
    data = sys.stdin.read().strip().split()
    R1 = float(data[0])
    R2 = float(data[1])

    # r = smaller radius, R = larger radius
    r = min(R1, R2)
    R = max(R1, R2)

    # Integrand f(z) = sqrt((r^2 - z^2)(R^2 - z^2)) for z in [0, r]
    def f(z: float) -> float:
        a = r*r - z*z
        b = R*R - z*z
        # Clamp to 0 to avoid tiny negative due to floating point round-off
        if a < 0.0: a = 0.0
        if b < 0.0: b = 0.0
        return math.sqrt(a * b)

    # Simpson's rule setup (n must be even)
    n = 1_000_000
    h = r / n

    # Simpson weighted sum
    s = f(0.0) + f(r)
    for i in range(1, n):
        z = i * h
        s += (4.0 if i % 2 == 1 else 2.0) * f(z)

    integral = (h / 3.0) * s
    volume = 8.0 * integral

    sys.stdout.write(f"{volume:.4f}\n")

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

Both implementations follow the same derived integral and use Simpson’s rule to meet the \(10^{-4}\) accuracy requirement.