Project Euler 283

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

Project Euler 283

题目

Integer sided triangles for which the area/perimeter ratio is integral

Consider the triangle with sides \(6, 8\) and \(10\). It can be seen that the perimeter and the area are both equal to \(24\).

So the area/perimeter ratio is equal to \(1\).

Consider also the triangle with sides \(13, 14\) and \(15\). The perimeter equals \(42\) while the area is equal to \(84\).

So for this triangle the area/perimeter ratio is equal to \(2\).

Find the sum of the perimeters of all integer sided triangles for which the area/perimeter ratios are equal to positive integers not exceeding \(1000\).

解决方案

设\(M=1000\)，三边为整数 \(a,b,c\)，周长 \(P=a+b+c\)，半周长\(s=\dfrac{P}{2}\)，面积 \(A\)。

海伦公式的等价形式为

\[16A^2=(a+b+c)(-a+b+c)(a-b+c)(a+b-c).\]

题目要求面积周长比为正整数且不超过 \(1000\)。令\(\dfrac{A}{P}=m, 1\le m\le 1000\)

则 \(A=mP=m(a+b+c)\)。代入上式并约去一项 \((a+b+c)\) 得

\[ (-a+b+c)(a-b+c)(a+b-c)=16m^2(a+b+c). \tag{1} \]

观察三个因子\(X_1=-a+b+c, X_2=a-b+c, X_3=a+b-c\)，它们两两之差为\(X_2-X_1=2a, X_3-X_2=2b, X_3-X_1=2b\)，因此 \(X_1,X_2,X_3\) 同奇偶。

而式 \((1)\) 右侧含 \(16m^2\)，必为偶数，所以左侧乘积为偶数，从而每个因子都必须为偶数。于是可令

\[2x=-a+b+c,2y=a-b+c,2z=a+b-c\tag{2}\]

把 \((2)\) 相加两两消元可解出三边：

\[a=y+z, b=x+z, c=x+y. \tag{3}\]

并且周长为

\[P=a+b+c=2(x+y+z), s=x+y+z. \tag{4}\]

注意：只要 \(x,y,z>0\)，由 \((3)\) 得到的 \(a,b,c\) 自动满足三角形不等式（例如 \(a=y+z<(x+z)+(x+y)=2x+y+z\) 恒成立），因此不会产生退化或非法三角形。

由海伦公式 \(A^2=s(s-a)(s-b)(s-c)\)，结合 \((3),(4)\)，可得\(s-a=x, s-b=y, s-c=z\)。

因此

\[ A^2=s\cdot x\cdot y\cdot z=(x+y+z)xyz. \tag{5} \]

另一方面 \(A=mP=2m(x+y+z)\)，所以

\[ A^2=4m^2(x+y+z)^2. \tag{6} \]

联立 \((5),(6)\)，并约去 \(x+y+z>0\)，得到核心丢番图方程：

\[ xyz=4m^2(x+y+z). \tag{7} \]

令\(n=4m^2\)，则 \((7)\) 写成

\[ xyz=n(x+y+z). \tag{8} \]

可以得到：

\[ xyz-ny-nz=nx. \tag{9} \]

注意左边形如 “关于 \(n\) 的二次式差一点就能配方”。对 \((9)\) 先乘 \(x\) 再加 \(n^2\)：

左右分别得到：

\[ \begin{aligned} x(xyz-ny-nz)+n^2&=x^2yz-nx(y+z)+n^2\\ &=n^2-n(xy+xz)+x^2yz\\ &=(n-xy)(n-xz)\\ &=(xy-n)(xz-n)\\ x\cdot nx+n^2&=nx^2+n^2\\ &=n(n+x^2) \end{aligned} \]

于是得到最终方程为

\[ (xy-n)(xz-n)=n(n+x^2). \tag{10} \]

由于方程 \((7)\) 对 \(x,y,z\) 对称。为避免同一三角形被不同排列重复计数，可强制排序\(x\le y\le z\)。

由 \((7)\) 与 \(x\le y\le z\) 得 \(x+y+z\le 3z\)，所以

\[ xyz=n(x+y+z)\le 3nz \Longrightarrow xy\le 3n. \tag{11} \]

再由 \(x\le y\) 得 \(x^2\le xy\le 3n\)，于是\(x\le \sqrt{3n}=\sqrt{12}m\)

这给出枚举 \(x\) 的有限范围。

固定 \(m\)（因此固定 \(n=4m^2\)）与某个 \(x\)，令\(R=n(n+x^2)\)。由 \((10)\) 可设\(d_1=xy-n, d_2=xz-n\)，则\(d_1d_2=R\)。

因此只要枚举 \(R\) 的正约数对 \((d_1,d_2)\)（满足 \(d_1d_2=R\)），即可反解\(y=\dfrac{d_1+n}{x}, z=\dfrac{d_2+n}{x}\)。

只保留满足以下条件的解：

  • \(y,z\) 为正整数（即 \(x\mid(d_1+n)\) 且 \(x\mid(d_2+n)\)）；
  • 排序约束 \(x\le y\le z\)（用于去重）。

一旦得到 \((x,y,z)\)，由 \((3)\) 可得三边\((a,b,c)=(y+z,x+z,x+y)\)

周长由 \((4)\) 直接为\(P=2(x+y+z)\)。

对每个 \(m\le M\)，有

而由于 \(x^2\le 3n\)，因此\(n+x^2\le 4n\le 16m^2\le 16\cdot 10^6.\)

也就是说，\((10)\) 里的两个因子 \(n\) 与 \(n+x^2\) 都不超过 \(16m^2\) 的量级。于是可以预处理一个“最小质因子”数组到上界，快速分解每个 \(n\) 与 \(n+x^2\)，再合并得到 \(R=n(n+x^2)\) 的质因数分解，从而生成全部约数用于产生 \(y,z\)。

（这正对应你给的实现备注：用最小质因子表来提取所需分解与约数。）

代码

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
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
#include <bits/stdc++.h>
#include "tools.h"
using namespace std;

// ------------------------------------------------------------
// 题意背景（从你的代码推断）：
// 你在枚举 m 与 x，并利用数论变形：
//   (x y - n)(x z - n) = n (n + x^2)
// 其中 n = 4 m^2。
// 通过令 g = gcd(x, n)，把公共因子约掉，得到更“原始”的等价形式：
//   (x1*y - n1)(x1*z - n1) = n1 (n1 + x1 * x)
// 其中 x = g*x1, n = g*n1。
// 然后枚举右侧 R = n1 (n1 + x1*x) 的所有因子 d1，
// 令 d2 = R / d1，则可恢复：
//   y = (d1 + n1) / x1
//   z = (d2 + n1) / x1
// 还需要整除性检查： (d1 + n1) % x1 == 0
// 最后加上约束 x <= y <= z 来避免重复计数（以及满足题目中隐含顺序条件）。
//
// 你使用 Factorizer 做快速质因数分解，并通过 merge_add_exp 合并两边因子指数，
// 再用 divisors(fc, false) 生成所有约数（false 表示不排序/或不做某种额外处理）。
// ------------------------------------------------------------

const int M = 1000;  // m 的上界（外层枚举到 1000）

// ------------------------------------------------------------
// merge_add_exp
// 输入：两个“质因数分解结果” fa, fb
//   fa: vector<pair<int,int>> 形式，每项是 (prime, exponent)
//   fb: 同上
// 且假设 fa 与 fb 都按 prime 升序排列（factor_small 通常会保证）。
//
// 输出：fc，表示两者指数相加后的结果：
//   对同一个质数 p，指数 e = ea + eb
//   若 e > 0 则写入 fc
//
// 用途：你要分解的目标数是：
//   mul = n1 * (n1 + x1*x)
// 先分别分解 n1 与 (n1 + x1*x)，得到 fa 与 fb，
// 合并指数即可得到 mul 的整体分解 fc。
// ------------------------------------------------------------
vector<pl> merge_add_exp(const vector<pi> &fa, const vector<pi> &fb)
{
    vector<pl> fc;
    // 预留空间：最坏情况下两个分解没有公共质因子，长度相加
    fc.reserve(fa.size() + fb.size());

    size_t i = 0, j = 0;
    // 双指针归并（类似 merge in merge-sort）
    while (i < fa.size() && j < fb.size())
    {
        if (fa[i].first == fb[j].first)
        {
            // 同一个质数：指数相加
            int e = fa[i].second + fb[j].second;
            if (e)
                fc.push_back({fa[i].first, e});
            ++i;
            ++j;
        }
        else if (fa[i].first < fb[j].first)
        {
            // fa 的当前质数更小：直接加入 fa[i]
            if (fa[i].second)
                fc.emplace_back(fa[i]); // 等价于 push_back，但可原地构造
            ++i;
        }
        else
        {
            // fb 的当前质数更小：直接加入 fb[j]
            if (fb[j].second)
                fc.emplace_back(fb[j]);
            ++j;
        }
    }

    // 收尾：fa 剩余部分
    while (i < fa.size())
    {
        if (fa[i].second)
            fc.emplace_back(fa[i]);
        ++i;
    }

    // 收尾：fb 剩余部分
    while (j < fb.size())
    {
        if (fb[j].second)
            fc.emplace_back(fb[j]);
        ++j;
    }

    return fc;
}

int main()
{
    // 计时起点（你目前没打印耗时，但 start 已经取了）
    auto start = std::chrono::high_resolution_clock::now();

    // Factorizer：用于对 <= 16*M*M 的整数做快速质因数分解
    // 因为 n = 4m^2 <= 4*M^2，且还会分解 n1 + x1*x，大致也在同数量级。
    Factorizer fact = Factorizer(M * M * 16);

    long long ans = 0; // 最终累加答案（2*(x+y+z) 的总和）

    // ------------------------------------------------------------
    // 枚举 m：决定 n = 4m^2
    // ------------------------------------------------------------
    for (int m = 1; m <= M; m++)
    {
        // 输出进度（注意：这本身会极大拖慢运行速度）
        printf("%d\n", m);

        int n = m * m * 4; // n = 4m^2

        // ------------------------------------------------------------
        // 枚举 x：约束是 x^2 <= 3n
        // 这等价于 x <= floor(sqrt(3n))
        // ------------------------------------------------------------
        for (int x = 1; x * x <= 3 * n; x++)
        {
            // --------------------------------------------------------
            // 数论化简：
            // 原式：(x y - n)(x z - n) = n (n + x^2)
            //
            // 令 g = gcd(x, n)，写成 x = g*x1, n = g*n1
            // 则
            //   (g*x1*y - g*n1)(g*x1*z - g*n1) = g*n1 (g*n1 + g^2*x1^2)
            // 左边有 g^2，右边也有 g^2，可以整体约去 g^2，得到等价式：
            //   (x1*y - n1)(x1*z - n1) = n1 (n1 + x1*x)
            // （因为 x = g*x1，所以 x^2 = g^2*x1^2，右边可整理成 n1(n1 + x1*x)）
            //
            // 这个化简的好处：
            // - 把公共因子 g 去掉，减少无意义的重复
            // - 右侧变成一个固定的乘积 R = n1 (n1 + x1*x)
            //   只需枚举其因子对 (d1, d2) 就能恢复 (y, z)
            // --------------------------------------------------------
            int g = __gcd(x, n);
            int n1 = n / g, x1 = x / g;

            // 分解 n1
            vector<pi> fa = fact.factor_small(n1);

            // 分解 (n1 + x1*x)
            // 注意：这里用的是 x1*x 而不是 x1^2
            // 因为化简后右侧是 n1 (n1 + x1*x)
            vector<pi> fb = fact.factor_small(n1 + x1 * x);

            // 合并分解：得到 mul = n1 * (n1 + x1*x) 的整体质因数分解
            vector<pl> fc = merge_add_exp(fa, fb);

            // 根据质因数分解生成所有约数 d1
            // divisors(fc, false)：
            // - fc 给出 prime^exp
            // - 生成所有 ∏ p_i^{k_i} 形式的数
            // - false 很可能表示“不排序”或“不做额外处理”
            vector<ll> divisors_c = divisors(fc, false);

            // mul = n1 * (n1 + x1*x)，用于算 d2 = mul/d1
            ll mul = n1 * (n1 + 1ll * x1 * x);

            // --------------------------------------------------------
            // 枚举 d1 | mul：
            // 设 d2 = mul/d1
            // 我们希望把它们对应到：
            //   d1 = (x1*y - n1)
            //   d2 = (x1*z - n1)
            // 则：
            //   x1*y = d1 + n1   => y = (d1 + n1)/x1
            //   x1*z = d2 + n1   => z = (d2 + n1)/x1
            // 需要保证整除：
            //   (d1 + n1) % x1 == 0
            // 同理 z 其实也会自动为整数（因为 d2 = mul/d1，结构对称）
            // 但为了保险也可以检查 (d2 + n1) % x1 == 0（你这里没写）
            // --------------------------------------------------------
            for (ll d1 : divisors_c)
            {
                // y 必须是整数
                if ((d1 + n1) % x1 != 0)
                    continue;

                ll d2 = mul / d1;
                ll y = (d1 + n1) / x1;
                ll z = (d2 + n1) / x1;

                // 只统计满足顺序的解，避免 (y,z) 与 (z,y) 重复计数
                // 同时也符合题目里隐含的 x <= y <= z 约束（常见于无序三元组）
                if (x <= y && y <= z)
                {
                    // 题目最终要加 2*(x+y+z)
                    // “2*” 通常来自对称计数（例如 ±、镜像、或另一侧的对应解）
                    // 具体原因取决于原题，但你这里与 Java 版本保持一致。
                    ans += 2 * (x + y + z);
                }
            }
        }
    }

    // 输出最终答案
    printf("%lld\n", ans);

    // 你如果想输出耗时，可以加：
    // auto end = std::chrono::high_resolution_clock::now();
    // auto ms = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
    // cerr << "time(ms)=" << ms << "\n";
}