Project Euler 395

题目

Pythagorean tree

The Pythagorean tree is a fractal generated by the following procedure: Start with a unit square. Then, calling one of the sides its base (in the animation, the bottom side is the base):

1. Attach a right triangle to the side opposite the base, with the hypotenuse coinciding with that side and with the sides in a 3-4-5 ratio. Note that the smaller side of the triangle must be on the ‘right’ side with respect to the base (see animation).
2. Attach a square to each leg of the right triangle, with one of its sides coinciding with that leg.
3. Repeat this procedure for both squares, considering as their bases the sides touching the triangle.

The resulting figure, after an infinite number of iterations, is the Pythagorean tree.

It can be shown that there exists at least one rectangle, whose sides are parallel to the largest square of the Pythagorean tree, which encloses the Pythagorean tree completely.

Find the smallest area possible for such a bounding rectangle, and give your answer rounded to \(10\) decimal places.

解决方案

将图形进行多次迭代后，得到了以下大概的一个轮廓：

考虑使用广度优先搜索的过程进行迭代。

我设计了一个阈值：在迭代过程中，如果正方形的长度小于阈值，那么停止迭代。需要注意的是，阈值的选择需要得当，否则勾股树的左下方那一部分长度可能会被忽略。

在计算点的过程中还需要比较小心。以第一轮为例，正方形上方这将会产生一个直角三角形。直角顶点的坐标的计算方式是将直角三角形的斜边绕左边的点逆时针旋转约\(37°\)，然后再将边长缩减到原来的\(0.8\)倍得到。

代码

# include <bits/stdc++.h>
using namespace std;
const double eps = 1e-15;
double minx = 0, miny = 0, maxx = 1, maxy = 1;
double leftRotate = acos(0.8);
//直接将正方形的下边旋转约180-53度。
double rightRotate = -acos(0.6);
double tag = 5;
struct Square {
    double len, x, y, rad;
};
int main() {

    queue<Square> q;
    q.push({1, 0, 0, 0});
    while (!q.empty()) {
        Square prm = q.front();q.pop();
        //正方形左下角的点。
        double len = prm.len, ldx = prm.x, ldy = prm.y, rad = prm.rad;

        minx = min(minx, ldx);
        maxx = max(maxx, ldx);
        miny = min(miny, ldy);
        maxy = max(maxy, ldy);

        double cosrad = cos(rad);
        double sinrad = sin(rad);
        double dx = len * cosrad;
        double dy = len * sinrad;
        //正方形左上角的点。
        double lux = ldx - dy, luy = ldy + dx;
        //正方形位于上边上面的新点。
        double tx = lux + ((dx * 0.8) - (dy * 0.6)) * 0.8,ty = luy + ((dx * 0.6) + (dy * 0.8)) * 0.8;

        if (len * 0.8 > eps) {
            double pl = len * tag * 0.8;
            if (lux > maxx - pl || lux < minx + pl || luy > maxy - pl || luy < miny + pl) {
                q.push({len * 0.8, lux, luy, rad + leftRotate});
            }

            if (len * 0.6 > eps) {
                pl = len * tag * 0.6;
                if (tx > maxx - pl || tx < minx + pl || ty > maxy - pl || ty < miny + pl) {
                    q.push({len * 0.6, tx, ty, rad + rightRotate});
                }
            }
        }
    }
    double ans = (maxx - minx) * (maxy - miny);
    printf("%.10f\n", ans);
}