UTPC 2012 D問題

Tatsuya Yatagawa

2012-12-08

さて次はD問題です。

これは簡単な数学的知識があればなんとなく方針は立つと思います。一般にある図形を拡大・縮小、回転、平行移動だけで写すような変換は、

$$ \begin{bmatrix} a & b \\ c & d \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix} + \begin{bmatrix} tx \\ ty \end{bmatrix} = \begin{bmatrix} x’ \\ y' \end{bmatrix} $$

みたいな形の連立方程式でかけます。拡大・縮小の比率は問題文から0.5であることが分かっているので、適当な方法で回転角度と平行移動距離を求めれば問題は解けそうです。ここで注意したい（僕が間違えた）ところは、回転角度を求める時に変換前のある点と変換後の同じ点の座標の内積をとって、acosで回転角度を求めようとすると失敗するというところです。なぜでしょうか？

上記の方法で回転角を求めると求まる角度は0からπで、ここまでは問題ないのですが、問題は時計回りに回転したのか、反時計回りに回転したのかが分からなくなってしまいます。ですので、正しく0から2πの範囲で回転角を求めるために、atan2を使って求めた角度の差分をとって、それを回転角とする必要があります。

これさえ間違えなければ、あとはC問題よりもはるかに簡単だと思います。

ちなみに、もう少し違った発想で問題を解くこともできて、このような線形変換には不動点定理と呼ばれる定理が成り立つらしく、同じ変換をひたすら繰り返していくと図形の各頂点が同じ点に収束します。下に示したコードはこの方法を使って解いたもので、変換を求めたのち、図形の最初の2点だけを距離が一定以上近くなるまで変換し続けています。こっちの方が書くべきコードの量は少なめですかね。

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <ctime>
#include <cstring>
#include <climits>
#include <algorithm>
#include <map>
#include <set>
#include <stack>
#include <vector>
#include <string>
using namespace std;

#define REP(i,n) for(int i=0; i<n; i++)
#define REPA(i,s,e) for(int i=s; i<=e; i++)
#define REPD(i,s,e) for(int i=s; i>=e; i--)

typedef unsigned long long ULL;
typedef vector<int> VI;
typedef pair<int, int> PII;

template <typename T> T tmax(T t1, T t2) { return t1 > t2 ? t1 : t2; }
template <typename T> T tmin(T t1, T t2) { return t1 < t2 ? t1 : t2; }

class Point {
public:
    double x, y;
    Point() : x(0.0), y(0.0) {}
    Point(double _x, double _y)
        : x(_x), y(_y) {}
    Point(const Point& p)
        : x(p.x), y(p.y) {}

    Point& operator =(const Point& p) {
        this->x = p.x;
        this->y = p.y;
        return (*this);
    }

    Point operator +(const Point& p) {
        Point q;
        q.x = this->x + p.x;
        q.y = this->y + p.y;
        return q;
    }

    Point operator -(const Point& p) {
        Point q;
        q.x = this->x - p.x;
        q.y = this->y - p.y;
        return q;
    }

    double norm() {
        return sqrt(this->x * this->x + this->y * this->y);
    }

    Point& normalize() {
        double l = this->norm();
        this->x /= l;
        this->y /= l;
        return (*this);
    }

    double inner(Point p) {
        return this->x * p.x + this->y * p.y;
    }

    Point rotate(double theta) {
        Point q;
        q.x = this->x * cos(theta) - this->y * sin(theta);
        q.y = this->x * sin(theta) + this->y * cos(theta);
        return q;
    }

    Point scale(double s) {
        Point q;
        q.x = this->x * s;
        q.y = this->y * s;
        return q;
    }
};

const int N_MAX = 20;
Point pt1[N_MAX];
Point pt2[N_MAX];

int main() {
    int N;
    cin >> N;
    REP(i, N) {
        cin >> pt1[i].x;
        cin >> pt1[i].y;
    }

    REP(i, N) {
        cin >> pt2[i].x;
        cin >> pt2[i].y;
    }

    Point e1 = pt1[1] - pt1[0];
    Point e2 = pt2[1] - pt2[0];
    double s = e2.norm() / e1.norm();
    double theta = atan2(e2.y, e2.x) - atan2(e1.y, e1.x);
    Point trans = pt2[0] - pt1[0].rotate(theta).scale(s);

    while((pt1[1] - pt1[0]).norm() > 1.0e-8) {
        pt1[0] = pt1[0].rotate(theta).scale(s) + trans;
        pt1[1] = pt1[1].rotate(theta).scale(s) + trans;
    }
    printf("%f %fn", pt1[0].x, pt1[0].y);
}