p446.ans1
======================
3.14159265359

=================
p446.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;
};

const double EPS = 1e-9;
const double PI = acos(-1.0L);

int n;
vector<pair<double, double>> stars1, stars2;

void read() {
    cin >> n;
    stars1.resize(n);
    stars2.resize(n);
    cin >> stars1;
    cin >> stars2;
}

double normalize_angle(double angle) {
    if(angle < 0) {
        angle += 2 * PI;
    }

    angle = fmod(angle, 2 * PI);
    return angle;
}

pair<double, double> rotate_vector(pair<double, double> v, double angle) {
    double cos_a = cos(angle);
    double sin_a = sin(angle);
    return {v.first * cos_a - v.second * sin_a,
            v.first * sin_a + v.second * cos_a};
}

vector<pair<double, double>> get_all_vectors(
    const vector<pair<double, double>>& stars
) {
    int m = stars.size();
    double cx = 0, cy = 0;
    for(auto& s: stars) {
        cx += s.first;
        cy += s.second;
    }

    cx /= m;
    cy /= m;

    vector<tuple<double, double, pair<double, double>>> star_data;
    for(int i = 0; i < m; i++) {
        double dx = stars[i].first - cx;
        double dy = stars[i].second - cy;

        double angle;
        if(fabs(dx) < EPS && fabs(dy) < EPS) {
            angle = numeric_limits<double>::infinity();
        } else {
            angle = normalize_angle(atan2(dy, dx));
        }

        double distance = sqrt(dx * dx + dy * dy);
        star_data.push_back({angle, distance, {dx, dy}});
    }

    sort(star_data.begin(), star_data.end(), [](const auto& a, const auto& b) {
        if(get<0>(a) != get<0>(b)) {
            return get<0>(a) < get<0>(b);
        }
        return get<1>(a) < get<1>(b);
    });

    vector<pair<double, double>> vectors;
    for(auto& item: star_data) {
        vectors.push_back(get<2>(item));
    }

    return vectors;
}

bool vectors_match(
    const vector<pair<double, double>>& vectors1,
    const vector<pair<double, double>>& vectors2, double rotation_angle,
    int offset
) {
    int m = vectors1.size();
    vector<pair<double, double>> rotated(m);
    for(int i = 0; i < m; i++) {
        rotated[i] = rotate_vector(vectors2[i], rotation_angle);
    }

    for(int i = 0; i < m; i++) {
        auto v1 = vectors1[i];
        auto v2 = rotated[(i + offset) % m];
        if(fabs(v1.first - v2.first) > EPS ||
           fabs(v1.second - v2.second) > EPS) {
            return false;
        }
    }

    return true;
}

void solve() {
    // Each photo is a set of stars whose centroid is translation invariant, so
    // we represent each star by its vector from the centroid. Within a photo we
    // sort these vectors by (polar angle, length), giving a canonical cyclic
    // ordering that is the same up to rotation.
    //
    // A rotation maps vectors1 onto vectors2 with some cyclic offset. We anchor
    // on vectors1[0] and try matching it against each vectors2[i] of equal
    // length, which fixes a candidate rotation angle = angle(v1) - angle(v2).
    // For each candidate we verify the whole set matches under that rotation
    // and cyclic offset (checking both the angle and its 2*pi complement), then
    // keep the smallest valid angle. The answer lies in [0, pi], so for each
    // accepted rotation we also consider 2*pi - angle.

    if(n == 1) {
        cout << fixed << setprecision(11) << 0.0 << '\n';
        return;
    }

    auto vectors1 = get_all_vectors(stars1);
    auto vectors2 = get_all_vectors(stars2);

    double min_angle = numeric_limits<double>::infinity();

    for(int i = 0; i < n; i++) {
        auto v1 = vectors1[0];
        auto v2 = vectors2[i];

        double v1_len = sqrt(v1.first * v1.first + v1.second * v1.second);
        double v2_len = sqrt(v2.first * v2.first + v2.second * v2.second);

        if(v1_len < EPS || v2_len < EPS) {
            continue;
        }

        if(fabs(v1_len - v2_len) > EPS) {
            continue;
        }

        double v1_angle = atan2(v1.second, v1.first);
        double v2_angle = atan2(v2.second, v2.first);
        double rotation_angle = normalize_angle(v1_angle - v2_angle);

        if(vectors_match(vectors1, vectors2, rotation_angle, i) ||
           vectors_match(vectors1, vectors2, 2 * PI - rotation_angle, i)) {
            min_angle = min(min_angle, rotation_angle);
            min_angle = min(min_angle, 2 * PI - rotation_angle);
        }
    }

    cout << fixed << setprecision(11) << min_angle << '\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();
        // cout << "Case #" << test << ": ";
        solve();
    }

    return 0;
}

=================
p446.in1
======================
3
0.0 0.0
1.0 1.0
0.0 1.0
3.0 3.0
2.0 2.0
3.0 2.0

=================
p446.py
======================
import math

EPS = 1e-9


def vector_length(v):
    return math.sqrt(v[0] ** 2 + v[1] ** 2)


def vector_angle(v):
    return math.atan2(v[1], v[0])


def rotate_vector(v, angle):
    cos_a = math.cos(angle)
    sin_a = math.sin(angle)
    return (v[0] * cos_a - v[1] * sin_a, v[0] * sin_a + v[1] * cos_a)


def normalize_angle(angle):
    if angle < 0:
        angle += 2 * math.pi
    angle = angle % (2 * math.pi)
    return angle


def get_all_vectors(stars):
    n = len(stars)
    if n == 0:
        return []

    cx = sum(star[0] for star in stars) / n
    cy = sum(star[1] for star in stars) / n

    star_data = []
    for i in range(n):
        dx = stars[i][0] - cx
        dy = stars[i][1] - cy
        vector = (dx, dy)

        if abs(dx) < EPS and abs(dy) < EPS:
            angle = float("inf")
        else:
            angle = math.atan2(dy, dx)
            angle = normalize_angle(angle)

        distance = math.sqrt(dx * dx + dy * dy)
        star_data.append((angle, distance, vector))

    star_data.sort(key=lambda x: (x[0], x[1]))
    vectors = [item[2] for item in star_data]
    return vectors


def vectors_match(vectors1, vectors2, rotation_angle, offset):
    if len(vectors1) != len(vectors2):
        return False

    n = len(vectors1)
    rotated_vectors2 = [rotate_vector(v, rotation_angle) for v in vectors2]

    for i in range(n):
        v1 = vectors1[i]
        v2 = rotated_vectors2[(i + offset) % n]
        if abs(v1[0] - v2[0]) > EPS or abs(v1[1] - v2[1]) > EPS:
            return False
    return True


def solve_rotation_estimation():
    n = int(input())

    stars1 = []
    for _ in range(n):
        x, y = map(float, input().split())
        stars1.append((x, y))

    stars2 = []
    for _ in range(n):
        x, y = map(float, input().split())
        stars2.append((x, y))

    if n == 1:
        return 0.0

    vectors1 = get_all_vectors(stars1)
    vectors2 = get_all_vectors(stars2)

    min_angle = float("inf")

    for i in range(n):
        v1 = vectors1[0]
        v2 = vectors2[i]

        v1_len = vector_length(v1)
        v2_len = vector_length(v2)

        if v1_len < EPS or v2_len < EPS:
            continue

        if abs(v1_len - v2_len) > EPS:
            continue

        v1_angle = vector_angle(v1)
        v2_angle = vector_angle(v2)
        rotation_angle = v1_angle - v2_angle

        rotation_angle = normalize_angle(rotation_angle)
        if vectors_match(
            vectors1, vectors2, rotation_angle, i
        ) or vectors_match(
            vectors1, vectors2, 2 * math.pi - rotation_angle, i
        ):
            min_angle = min(min_angle, rotation_angle)
            min_angle = min(min_angle, 2 * math.pi - rotation_angle)

    if min_angle == float("inf"):
        raise ValueError("No valid rotation found")

    return min_angle


if __name__ == "__main__":
    try:
        result = solve_rotation_estimation()
        print(f"{result:.11f}")
    except ValueError as e:
        print(f"Error: {e}")
        raise

=================
statement.txt
======================
446. Rotation Estimation
Time limit per test: 0.75 second(s)
Memory limit: 262144 kilobytes
input: standard
output: standard

Mr. Nod is an astrologist and has defined a new constellation. He took two photos of the constellation to foretell a future of his friend. The constellation consists of n stars. The shape of the constellation in these photos are the same, but the angle of them are different because these photos were taken on a different day. He foretells a future by the difference of the angle of them. Your job is to write a program to calculate the difference of the angle of two constellation.
Input
The input is given in the following format:
n
x1,1 y1,1
...
x1,n y1,n
x2,1 y2,1
...
x2,n y2,n
The first line of the input contains a positive integer n (n ≤ 1000). The next n lines contain two real numbers x1,i and y1,i (|x1,i|, |y1,i| ≤ 100), where (x1,i, y1,i) denotes the coordinates of the i-th star of the constellation in the first photo. The next n lines contain two real numbers x2,i and y2,i (|x2,i|, |y2,i| ≤ 100), where (x2,i, y2,i) denotes the coordinates of the i-th star of the constellation in the second photo. Note that the ordering of the stars does not matter for the sameness. It is guaranteed that distance between every pair of stars in each photo is larger than 10-5.
Output
You should print a non-negative real number which is the difference of the angle of the constellation in the first photo and in the second photo. The difference should be in radian, and should not be negative. If there are two or more solutions, you should print the smallest one, i.e. your solution should be an angle between 0 and pi radians, inclusive. The difference may be printed with any number of digits after decimal point, provided the absolute error does not exceed 10-7.
Example(s)
sample input
sample output
3
0.0 0.0
1.0 1.0
0.0 1.0
3.0 3.0
2.0 2.0
3.0 2.0
3.14159265359



=================
