Project Euler 583

发表于 分类于 阅读次数： 本文字数： 3.3k 阅读时长 ≈ 3 分钟

Project Euler 583

题目

Heron Envelopes

A standard envelope shape is a convex figure consisting of an isosceles triangle (the flap) placed on top of a rectangle. An example of an envelope with integral sides is shown below. Note that to form a sensible envelope, the perpendicular height of the flap (\(BCD\)) must be smaller than the height of the rectangle (\(ABDE\)).

In the envelope illustrated, not only are all the sides integral, but also all the diagonals (\(AC, AD, BD, BE\) and \(CE\)) are integral too. Let us call an envelope with these properties a Heron envelope.

Let \(S(p)\) be the sum of the perimeters of all the Heron envelopes with a perimeter less than or equal to \(p\).

You are given that \(S(10^4) = 884680\). Find \(S(10^7)\).

解决方案

将信封放在平面直角坐标系中，令矩形 \(ABDE\) 的底边为 \(AE\)，上边为 \(BD\)，翻盖为等腰三角形 \(BCD\)。

设\(a=AE=BD,b=AB=D,c=BC=CD,h\) 为翻盖高（点 \(C\) 到直线 \(BD\) 的垂直距离）。取坐标\(A(0,0),E(a,0),B(0,b),D(a,b),C\left(\dfrac a2,b+h\right).\)

题目要求：五条对角线 \(AC,AD,BD,BE,CE\) 全为整数；外边 \(AB,BC,CD,DE,EA\) 全为整数且翻盖高必须小于矩形高，即\(h<b.\)外周长为\(P=a+2b+2c.\)由于对称性：\(BD=a, AD=BE, AC=CE,\)所以只需确保下列三类长度为整数：\(AD=\sqrt{a^2+b^2},AC=\sqrt{\left(\dfrac a2\right)^2+(b+h)^2}\)，以及翻盖腰 \(c\)（等价于 \(\left(\dfrac a2\right)^2+h^2=c^2\)）

由翻盖等腰三角形 \(BCD\)（高落在底边中点）得：\(c^2=\left(\dfrac a2\right)^2+h^2\Longrightarrow4c^2=a^2+(2h)^2.\)

对角线 \(AC\)满足：\(AC^2=\left(\dfrac a2\right)^2+(b+h)^2\Longrightarrow(b+h)^2=AC^2-\left(\dfrac a2\right)^2.\)右侧是整数减去有理数，因此 \((b+h)^2\in\mathbb Q\)，从而 \(b+h\in\mathbb Q\)（取非负根），故\(h\in\mathbb Q.\)

再看\((2h)^2=4c^2-a^2\in\mathbb Z.\)由 \(2h\in\mathbb Q\) 且 \((2h)^2\in\mathbb Z\) 可知 \(2h\in\mathbb Z\)，记\(m=2h\in\mathbb Z.\)将 \(4c^2=a^2+(2h)^2\) 写成\(a^2+m^2=(2c)^2.\)右侧必被 \(4\) 整除，故左侧 \(a^2+m^2\equiv 0\pmod 4\)。平方数模 \(4\) 只能是 \(0\) 或 \(1\)，要使和为 \(0\)，只能二者都为 \(0\)：\(a\equiv 0\pmod 2, m\equiv 0\pmod 2.\) 于是 \(a\) 为偶数，且 \(m\) 为偶数。由 \(m=2h\) 且 \(m\) 偶数，可得\(h=\dfrac m2\in\mathbb Z.\)

因此任何解必满足 \(a\) 为偶数、\(h\) 为整数。

引入 \(m=2h\)，则：矩形对角线整数\(a^2+b^2=AD^2.\)翻盖腰为整数：\(a^2+m^2=(2c)^2, m=2h.\)对角线 \(AC\) 为整数：\(4AC^2=a^2+(2b+m)^2\Longrightarrow a^2+(2b+m)^2=(2AC)^2.\)再加上信封可用条件：\(h<b\iff \dfrac m2<b\iff 2b>m.\)周长为\(P=a+2b+2c=a+2b+\sqrt{a^2+m^2}.\)

因此问题等价为：找所有整数 \(a,b,m\)（其中 \(a\) 偶数、\(m\) 偶数）使得\((a,b),(a,m),(a,2b+m)\)分别都能作为某个勾股三元组的两条直角边，并满足 \(2b>m\) 与 \(P\le N\)。

对固定的 \(a\)，考虑集合 \(L(a)\) 为所有满足 \(a^2+s^2\) 是完全平方数的正整数 \(s\) 的集合。那么三条条件等价于\(b\in L(a), m\in L(a), 2b+m\in L(a).\)于是对每个 \(a\)，只要拿到 \(L(a)\)，就可以在集合内部构造：选 \(m\in L(a)\)，且 \(m\) 必为偶数；选 \(t\in L(a)\)，且 \(t=2b+m\) 必为偶数，那么得\(b=\dfrac{t-m}{2}.\)检查 \(b\in L(a)\) 且 \(b>\dfrac m2\)，并且周长用\(P=a+2b+\sqrt{a^2+m^2}\)判断是否 \(\le N\)。

所有本原勾股三元组由如下公式给出：\(x=u^2-v^2, y=2uv, z=u^2+v^2,\)其中 \(u>v\)、\(\gcd(u,v)=1\)、且 \(u,v\) 一奇一偶。任意非本原三元组是其倍数 \(k(x,y,z)\)。

因为已证明 \(a\) 必为偶数，所以只需要把偶数直角边作为首边记录。对每个三元组 \((x,y,z)\) 及其倍数：若 \(x\) 为偶数，则记录一对 \((a=x,s=y)\)；若 \(y\) 为偶数，则记录一对 \((a=y,s=x)\)。把所有 \((a,s)\) 记录排序并按 \(a\) 分组，就得到每个 \(a\) 的集合 \(L(a)\)。

根据前面的整性结论，设 \(m=2h\)，则 \(m\) 必为偶数，并且 \(2b+m\) 也必为偶数。因此在 \(L(a)\) 中只需要关注偶数元素来承担 \(m\) 与 \(2b+m\) 的角色。

首先选取 \(m\in L(a)\)（\(m\) 为偶数）满足\(a^2+m^2=(2c)^2,\)从而\(2c=\sqrt{a^2+m^2}.\)当 \(a\) 与 \(m\) 固定后，周长可以写成\(P=a+2b+2c=a+2b+\sqrt{a^2+m^2}.\)因此对于该 \(m\)，周长中除了 \(2b\) 以外的部分已被完全确定。

接着在 \(L(a)\) 中选取另一个偶数 \(t\) 来对应\(t=2b+m,\)于是\(b=\frac{t-m}{2}.\)可折叠条件 \(h<b\) 等价于\(b>\frac m2,\)再代入 \(t=2b+m\) 得到\(t>2m.\)因此只需要枚举所有满足 \(t>2m\) 的候选即可。

得到 \(b\) 之后，还必须满足矩形对角线为整数，即\(a^2+b^2\)是完全平方数。这等价于 \(b\in L(a)\)。当该条件成立时，三条关键长度 \(AD,,2c,,2AC\) 同时为整数，从而构成一个合法的解。

最后利用单调性剪枝：当 \(m\) 固定时，随着 \(t\) 增大，\(b=(t-m)/2\) 单调增大，从而周长\(P=a+2b+\sqrt{a^2+m^2}\)也单调增大。一旦某个 \(t\) 使 \(P>N\)，后续更大的 \(t\) 不可能再满足周长限制，因此可立即停止该 \(m\) 下的枚举。这样对每个 \(a\) 的所有组合进行遍历，便能无重无漏地枚举所有解，并将其周长累加得到 \(S(N)\)。

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
#include <bits/stdc++.h>
#include "tools.h"
using namespace std;
const int N = 10000000;
typedef long long ll;

struct Key
{
    int a;
    int s;
    bool operator<(const Key &o) const
    {
        if (a != o.a)
            return a < o.a;
        return s < o.s;
    }
    bool operator==(const Key &o) const
    {
        return a == o.a && s == o.s;
    }
};

ll solve(int N)
{
    vector<Key> keys;
    int lim = int_square(N) + 2;
    for (int m = 2; m <= lim; ++m)
    {
        ll mm = 1LL * m * m;
        for (int n = 1; n < m; ++n)
        {
            ll nn = 1LL * n * n;
            ll c0 = mm + nn;
            if (c0 > N)
                break;
            if (((m + n) & 1) == 0)
                continue;
            if (__gcd(m, n) != 1)
                continue;

            ll a0 = mm - nn;
            ll b0 = 2LL * m * n;

            for (ll k = 1; k * c0 <= N; ++k)
            {
                ll x = k * a0;
                ll y = k * b0;
                if ((x & 1) == 0)
                    keys.push_back(Key{(int)x, (int)y});
                if ((y & 1) == 0)
                    keys.push_back(Key{(int)y, (int)x});
            }
        }
    }

    sort(keys.begin(), keys.end());
    keys.erase(unique(keys.begin(), keys.end()), keys.end());

    ll ans = 0;
    size_t i = 0;

    vector<int> all;
    vector<int> ev;

    while (i < keys.size())
    {
        int a = keys[i].a;
        size_t j = i;
        while (j < keys.size() && keys[j].a == a)
            ++j;

        all.clear();
        all.reserve(j - i);
        for (size_t t = i; t < j; ++t)
            all.push_back(keys[t].s);

        ev.clear();
        ev.reserve(all.size());
        for (auto s : all)
            if ((s & 1) == 0)
                ev.push_back(s);

        for (size_t u = 0; u < ev.size(); ++u)
        {
            int m2h = ev[u];
            if (!m2h)
                continue;

            ll ss = 1LL * a * a + 1LL * m2h * m2h;
            ll hyp = int_square(ss);
            if (hyp * hyp != ss)
                continue;

            ll base = (ll)a + hyp;

            int thr = m2h * 2 + 2;
            auto it0 = lower_bound(ev.begin(), ev.end(), thr);

            int bmin = m2h / 2 + 1;
            size_t pb = lower_bound(all.begin(), all.end(), bmin) - all.begin();

            for (auto it = it0; it != ev.end(); ++it)
            {
                int t2 = *it;
                int diff = t2 - m2h;
                if (diff & 1)
                    continue;
                int b = diff >> 1;
                if (b <= m2h / 2)
                    continue;

                ll per = base + 2LL * b;
                if (per > N)
                    break;

                while (pb < all.size() && all[pb] < b)
                    ++pb;
                if (pb < all.size() && all[pb] == b)
                    ans += per;
            }
        }

        i = j;
    }
    return ans;
}
int main()
{
    cout << solve(N) << "\n";
    return 0;
}