POJ 3845 Fractal (向量旋转，好题)

#include <iostream>
#include <stdio.h>
#include <string.h>
#include <algorithm>
#include <queue>
#include <map>
#include <vector>
#include <set>
#include <string>
#include <math.h>

#define LL long long
#define mem(i, j) memset(i, j, sizeof(i))
#define rep(i, j, k) for(int i = j; i <= k; i++)
#define dep(i, j, k) for(int i = k; i >= j; i--)
#define pb push_back
#define make make_pair
#define INF INT_MAX
#define inf LLONG_MAX
#define PI acos(-1)
#define fir first
#define sec second
using namespace std;

const int N = 110;
const double eps = 1e-8;
const double maxL = 10.0;

struct Point {
    double x, y;
    Point(double x = 0, double y = 0) : x(x), y(y) { }
};

int dcmp(double x) {
    if(fabs(x) < eps) return 0; else return x < 0 ? -1 : 1;
}

Point operator + (Point A, Point B) { return Point(A.x + B.x, A.y + B.y); }
Point operator - (Point A, Point B) { return Point(A.x - B.x, A.y - B.y); }
Point operator * (Point A, double p) { return Point(A.x * p, A.y * p); }
Point operator / (Point A, double p) { return Point(A.x / p, A.y / p); }

double Cross(Point A, Point B) { return A.x * B.y - A.y * B.x; }
double Dot(Point A, Point B) { return A.x * B.x + A.y * B.y; }
double Length(Point A) { return sqrt(Dot(A, A)); }

Point Rotate(Point A, double rad) { /// 向量逆时针旋转 rad (弧度)
    return Point(A.x*cos(rad) - A.y*sin(rad), A.x*sin(rad) + A.y*cos(rad));
}


/* 有向直线，它的左边就是对应的半平面 */
struct Line {
    Point p; /// 直线任意一点
    Point v; /// 方向向量
    double ang; /// 极角，即从x正半轴旋转到向量v所需要的角（弧度）
    Line() { }
    Line(Point p, Point v) : p(p), v(v) { ang = atan2(v.y, v.x); }
    bool operator < (const Line& L) const {
        return ang < L.ang;
    }
};

Point P[N], ans;

Line L[N];

double sum[N], tmp, Mul[N];

int n, d;
double f;

void dfs(Point A, Point B, double len, int depth) {
    if(depth == 0) { 
        /// (B - A) 是向量长度，乘以比例再加A，就是点A沿着向量方向前进len之后的那个点
        ans = (B - A) * len / Length(B - A) + A;
        return ;
    }

    rep(i, 1, n - 1) {
        /// 找到第一个大于 len 的点
        if(dcmp(sum[i] * Length(B - A) / Length(P[n - 1] - P[0]) * Mul[depth] - len) < 0) continue;
        ///最初的极角
        double Angle1 = atan2(P[n - 1].y - P[0].y, P[n - 1].x - P[0].x);
        ///当前边的极角
        double Angle2 = atan2(B.y - A.y, B.x - A.x);
        ///旋转
        Point st = Rotate(P[i - 1] - P[0], Angle2 - Angle1);
        ///缩减
        st = st * Length(B - A) / Length(P[n - 1] - P[0]) + A;
        Point ed = Rotate(P[i] - P[0], Angle2 - Angle1);
        ed = ed * Length(B - A) / Length(P[n - 1] - P[0]) + A;
        dfs(st, ed, len - sum[i - 1] * Mul[depth] * Length(B - A) / Length(P[n - 1] - P[0]), depth - 1);
        return ;
    }

}

void solve() {


    scanf("%d", &n);
    rep(i, 0, n - 1) scanf("%lf %lf", &P[i].x, &P[i].y);
    scanf("%d %lf", &d, &f);

    rep(i, 1, n - 1) sum[i] = sum[i - 1] + Length(P[i - 1] - P[i]);

    tmp = sum[n - 1] / Length(P[n - 1] - P[0]);
    Mul[1] = 1;
    rep(i, 2, d) Mul[i] = Mul[i - 1] * tmp;
    
    /// sum[n-1] * Mul[d] 是折线变了 d 次后的总长
    dfs(P[0], P[n - 1], sum[n - 1] * Mul[d] * f, d);

    printf("(%.10f,%.10f)\n", ans.x, ans.y);

}

int main() {

    int _; scanf("%d", &_);
    while(_--) solve();

    return 0;
}