TopCoder SRM 500 Div1 Medium FractalPicture

問題

次のようにしてフラクタルな図形を描く。

  • 第一世代: (0, 0)から(0, 81)の世代
  • 第二世代: 前の線分の上1/3を三つに枝分かれさせる((0, 54)から、(0, 81), (-27, 54), (27, 54)に)
  • 第三世代: 枝分かれした三つの枝を、更に同様に三つずつに枝分かれさせる
  • 以下500世代まで繰り返し

この図形のうち、(x1, y1), (x2, y2)の長方形の領域に含まれる部分の長さを求めよ。

制約条件

入力は整数。

方針

第n世代の枝の部分の長さを前計算で計算しておく。

再帰により図形の長さを求める。
単純に再帰すると全然終わらないが、
図形は、辺の比が1:3の長方形の内部にすっぽり納まることを使って、
図形が求める領域の中にすっぽり入るときはその長さを返してしまうようにする、
また、図形を覆う長方形と領域が共有点を持たないときは即座に0を返すようにすると、
繰り返し回数が少なくなってすぐに終わる。


長方形の座標が整数なので実はEPSとかは要らなかったりする。

ソースコード

typedef long double ld;
typedef complex<ld> P;
int x, X, y, Y;
ld len[510];
inline ld calc(const P &s, const P &t){
	G g;
	g.pb(s); g.pb(t);
	const int xs[] = {x, X}, ys[] = {y, Y};
	P p[4];
	rep(i, 2) rep(j, 2) p[i * 2 + j] = P(xs[i], ys[i ? j : 1 - j]);
	rep(i, 4) g = convex_cut(g, L(p[i], p[(i + 1) % 4]));
	
	if(g.empty()) return 0;
	return abs(g[0] - g[1]);
}
ld rec(int n, const P &s, const P &t){
	
	P m = (s + t * 2.0l) / 3.0l;
	P d = (t - s) / 3.0l;
	
	P l = m + d * P(0, 1), r = m - d * P(0, 1);
	ld x1 = min(min(l.real(), r.real()), min(s.real(), t.real()));
	ld x2 = max(max(l.real(), r.real()), max(s.real(), t.real()));
	ld y1 = min(min(l.imag(), r.imag()), min(s.imag(), t.imag()));
	ld y2 = max(max(l.imag(), r.imag()), max(s.imag(), t.imag()));
	
	if(x2 <= x || x1 >= X ||
		y2 <= y || y1 >= Y ) return 0;
	
	if(x < x1 + EPS && x2 < X + EPS &&
		y < y1 + EPS && y2 < Y + EPS) return len[n];
	
	if(n == 499){
		return calc(s, t);
	}
	
	ld res = calc(s, m);
	
	res += rec(n + 1, m, l);
	res += rec(n + 1, m, t);
	res += rec(n + 1, m, r);
	return res;
}

class FractalPicture {
	public:
	double getLength(int x1, int y1, int x2, int y2) 
	{
		rep(i, 500){
			len[i] = 81 + 54 * (499 - i);
			rep(j, i) len[i] /= 3;
		}
		x = x1; X = x2;
		y = y1; Y = y2;
		
		return (double)rec(0, P(0, 0), P(0, 81));
	}
};