Test/UC/3-37-M.test.cpp

• View this file on GitHub
• Last update: 2025-05-21 20:07:18+09:00
• Problem: https://onlinejudge.u-aizu.ac.jp/courses/lesson/2/ITP1/1/ITP1_1_A
• Link: https://contest.ucup.ac/contest/2025/problem/10748
• Link: https://contest.ucup.ac/submission/1067790

Depends on

• Src/GeometryR3/Plane.hpp
• Src/GeometryR3/Point.hpp
• 標準データ型のエイリアス (Src/Template/TypeAlias.hpp)

Code

// #define PROBLEM "https://contest.ucup.ac/contest/2025/problem/10748"
#define PROBLEM "https://onlinejudge.u-aizu.ac.jp/courses/lesson/2/ITP1/1/ITP1_1_A"

/*
 * The 3rd Universal Cup. Stage 37: Wuhan M. Flight Tracker
 * https://contest.ucup.ac/submission/1067790
 */

#include "../../Src/GeometryR3/Point.hpp"
#include "../../Src/GeometryR3/Plane.hpp"

#include <algorithm>
#include <cmath>
#include <iostream>
#include <iomanip>

using namespace zawa::geometryR3;
using point = Point<long double>;

long double R;
point P, S, T;

long double CenterAngle(const point& s, const point& t) {
    const long double cosT = 1.0l - DistanceSquare(s, t) / (2 * R * R);
    return acosl(std::clamp(cosT, -1.0l, 1.0l));
}

long double ArcLength(const point& p, const point& q) {
    return R * CenterAngle(p, q);
}

long double solve() {
    P = R * P.normalized(); 
    S = R * S.normalized(); 
    T = R * T.normalized(); 
    const long double theta = CenterAngle(S, T);
    auto ST = [&](long double r) -> point {
        const point n = Plane{S, T}.normalVector();
        return S.rotated(n, theta * r);
    };
    long double l = 0, r = 1;
    for (int t = 0 ; t < 50 ; t++) {
        const long double m1 = (l + l + r) / 3, m2 = (l + r + r) / 3;
        const long double d1 = ArcLength(P, ST(m1)), d2 = ArcLength(P, ST(m2));
        if (d1 <= d2) r = m2;
        else l = m1;
    }
    return ArcLength(P, ST(l));
}

int main() {
#ifdef ONLINE_JUDGE
    std::cin.tie(nullptr);
    std::cout.tie(nullptr);
    std::ios::sync_with_stdio(false);
    int t;
    std::cin >> t;
    std::cout << std::fixed << std::setprecision(7);
    while (t--) {
        std::cin >> R >> P >> S >> T;
        std::cout << solve() << '\n';
    }
#else
    std::cout << "Hello World\n";
#endif
}

#line 1 "Test/UC/3-37-M.test.cpp"
// #define PROBLEM "https://contest.ucup.ac/contest/2025/problem/10748"
#define PROBLEM "https://onlinejudge.u-aizu.ac.jp/courses/lesson/2/ITP1/1/ITP1_1_A"

/*
 * The 3rd Universal Cup. Stage 37: Wuhan M. Flight Tracker
 * https://contest.ucup.ac/submission/1067790
 */

#line 2 "Src/GeometryR3/Point.hpp"

#include <array>
#include <cassert>
#include <cmath>
#include <concepts>
#include <iostream>
#include <utility>

#line 2 "Src/Template/TypeAlias.hpp"

#include <cstdint>
#include <cstddef>

namespace zawa {

using i16 = std::int16_t;
using i32 = std::int32_t;
using i64 = std::int64_t;
using i128 = __int128_t;

using u8 = std::uint8_t;
using u16 = std::uint16_t;
using u32 = std::uint32_t;
using u64 = std::uint64_t;

using usize = std::size_t;

} // namespace zawa
#line 11 "Src/GeometryR3/Point.hpp"

namespace zawa {

namespace geometryR3 {

template <std::floating_point T>
class Point {
public:

    Point() = default;

    Point(T x, T y, T z) : m_p{x, y, z} {}

    explicit Point(const std::array<T, 3>& p) : m_p{p} {}

    explicit Point(std::array<T, 3>&& p) : m_p{std::move(p)} {}

    const T& x() const {
        return m_p[0];
    }

    T& x() {
        return m_p[0];
    }

    const T& y() const {
        return m_p[1];
    }

    T& y() {
        return m_p[1];
    }

    const T& z() const {
        return m_p[2];
    }

    T& z() {
        return m_p[2];
    }

    const T& operator[](usize i) const {
        assert(i < 3);
        return m_p[i];
    }

    T& operator[](usize i) {
        assert(i < 3);
        return m_p[i];
    }

    Point<T>& operator+=(const Point<T>& p) {
        m_p[0] += p[0];
        m_p[1] += p[1];
        m_p[2] += p[2];
        return *this;
    }

    Point<T>& operator-=(const Point<T>& p) {
        m_p[0] -= p[0];
        m_p[1] -= p[1];
        m_p[2] -= p[2];
        return *this;
    }

    Point<T>& operator*=(T k) {
        m_p[0] *= k;
        m_p[1] *= k;
        m_p[2] *= k;
        return *this;
    }

    Point<T>& operator/=(T k) {
        m_p[0] /= k;
        m_p[1] /= k;
        m_p[2] /= k;
        return *this;
    }

    T normSquare() const {
        return m_p[0]*m_p[0] + m_p[1]*m_p[1] + m_p[2]*m_p[2];
    }

    T norm() const {
        return sqrtl(normSquare());
    }

    Point<T> normalized() const {
        return Point<T>{*this} /= norm();
    }

    // ロドリゲスの回転公式 https://manabitimes.jp/math/2649
    Point<T> rotated(const Point<T>& axis, T theta) const {
        const T cosT = cosl(theta), sinT = sinl(theta);
        const Point<T> cp{*this};
        return cosT * cp + (1 - cosT) * DotProduct(axis, cp) * axis + sinT * CrossProduct(axis, cp);
    }

private:

    std::array<T, 3> m_p{T{0}, T{0}, T{0}};

};

template <std::floating_point T>
Point<T> operator+(const Point<T>& lhs, const Point<T>& rhs) {
    return Point<T>{lhs} += rhs;
}

template <std::floating_point T>
Point<T> operator-(const Point<T>& lhs, const Point<T>& rhs) {
    return Point<T>{lhs} -= rhs;
}

template <std::floating_point T>
Point<T> operator*(const T k, const Point<T>& p) {
    return Point<T>{p} *= k;
}

template <std::floating_point T>
std::ostream& operator<<(std::ostream& os, const Point<T>& p) {
    os << '(' << p[0] << ',' << p[1] << ',' << p[2] << ')';
    return os;
}

template <std::floating_point T>
std::istream& operator>>(std::istream& is, Point<T>& p) {
    is >> p[0] >> p[1] >> p[2];
    return is;
}

template <std::floating_point T>
T DistanceSquare(const Point<T>& lhs, const Point<T>& rhs) {
    return Point{lhs - rhs}.normSquare();
}

template <std::floating_point T>
T Distance(const Point<T>& lhs, const Point<T>& rhs) {
    return sqrtl(DistanceSquare(lhs, rhs));
}

template <std::floating_point T>
T DotProduct(const Point<T>& lhs, const Point<T>& rhs) {
    return lhs[0]*rhs[0] + lhs[1]*rhs[1] + lhs[2]*rhs[2];
}

template <std::floating_point T>
Point<T> CrossProduct(const Point<T>& lhs, const Point<T>& rhs) {
    return {
        lhs[1]*rhs[2]-lhs[2]*rhs[1],
        lhs[2]*rhs[0]-lhs[0]*rhs[2],
        lhs[0]*rhs[1]-lhs[1]*rhs[0],
    };
}

template <std::floating_point T>
using Vector = Point<T>;

} // namespace geometryR3

} // namespace zawa
#line 2 "Src/GeometryR3/Plane.hpp"

#line 4 "Src/GeometryR3/Plane.hpp"

#line 6 "Src/GeometryR3/Plane.hpp"

namespace zawa {

namespace geometryR3 {

template <std::floating_point T>
class Plane {
public: 

    Plane() = default;

    Plane(const Vector<T>& p, const Vector<T>& q) : m_p{p}, m_q{q} {}

    Plane(const Point<T>& p, const Point<T>& q, const Point<T>& r) : m_p{q - p}, m_q{r - p} {}

    Vector<T> normalVector() const {
        return CrossProduct(m_p, m_q).normalized();
    }

private:

    Vector<T> m_p, m_q;
};

} // namespace geometryR3

} // namespace zawa
#line 11 "Test/UC/3-37-M.test.cpp"

#include <algorithm>
#line 15 "Test/UC/3-37-M.test.cpp"
#include <iomanip>

using namespace zawa::geometryR3;
using point = Point<long double>;

long double R;
point P, S, T;

long double CenterAngle(const point& s, const point& t) {
    const long double cosT = 1.0l - DistanceSquare(s, t) / (2 * R * R);
    return acosl(std::clamp(cosT, -1.0l, 1.0l));
}

long double ArcLength(const point& p, const point& q) {
    return R * CenterAngle(p, q);
}

long double solve() {
    P = R * P.normalized(); 
    S = R * S.normalized(); 
    T = R * T.normalized(); 
    const long double theta = CenterAngle(S, T);
    auto ST = [&](long double r) -> point {
        const point n = Plane{S, T}.normalVector();
        return S.rotated(n, theta * r);
    };
    long double l = 0, r = 1;
    for (int t = 0 ; t < 50 ; t++) {
        const long double m1 = (l + l + r) / 3, m2 = (l + r + r) / 3;
        const long double d1 = ArcLength(P, ST(m1)), d2 = ArcLength(P, ST(m2));
        if (d1 <= d2) r = m2;
        else l = m1;
    }
    return ArcLength(P, ST(l));
}

int main() {
#ifdef ONLINE_JUDGE
    std::cin.tie(nullptr);
    std::cout.tie(nullptr);
    std::ios::sync_with_stdio(false);
    int t;
    std::cin >> t;
    std::cout << std::fixed << std::setprecision(7);
    while (t--) {
        std::cin >> R >> P >> S >> T;
        std::cout << solve() << '\n';
    }
#else
    std::cout << "Hello World\n";
#endif
}

Back to top page