Project Euler 370

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

Project Euler 370

题目

Geometric triangles

Let us define a geometric triangle as an integer sided triangle with sides \(a \le b \le c\) so that its sides form a geometric progression, i.e. \(b^2=a\cdot c\).

An example of such a geometric triangle is the triangle with sides \(a = 144, b = 156\) and \(c = 169\).

There are \(861805\) geometric triangles with perimeter \(\le 10^6\) .

How many geometric triangles exist with perimeter \(\le 2.5\cdot10^{13}\) ?

解决方案

Project Euler 370 题解：Geometric triangles

题目要求整数边三角形 \(a\le b\le c\)，且三边成等比数列，即\(b^2=a\cdot c.\)令 \(g=\gcd(a,c)\)，写成 \(a=gu,c=gv,\gcd(u,v)=1.\)代回 \(b^2=ac\) 得\(b^2=g^2uv.\)

因此 \(u,v\) 是完全平方数。由于 \(\gcd(u,v)=1\)，只能是\(u=n_0^2, v=x_0^2, \gcd(n_0,x_0)=1,\)从而\(a=g,n_0^2, b=g,n_0x_0, c=g,x_0^2.\)

接下来把 \(g\) 拆成无平方因子部分 \(\times\) 平方数”：\(g=st^2,\)其中 \(s\) 为无平方因子（squarefree），\(t\ge 1\) 为整数。于是\(a=s(tn_0)^2, b=s(tn_0)(tx_0), c=s(tx_0)^2.\)令\(n=tn_0, x=tx_0,\)得到统一形式

\[ (a,b,c)=(s n^2,s n x,s x^2), \]

其中 \(s\) 是无平方因子数，且 \(n,x\) 为正整数。

唯一性说明：\(g=\gcd(a,c)\) 唯一确定；\(g=s t^2\) 的分解（\(s\) 取无平方因子部分）唯一；\(u=a/g\) 与 \(v=c/g\) 唯一，进而 \(n_0=\sqrt u,x_0=\sqrt v\) 唯一；于是 \(n=tn_0,x=tx_0\) 唯一。因此每个几何三角形对应唯一三元组 \((n,x,s)\)。

对边\(a=s n^2, b=s n x, c=s x^2\)只需检查 \(a+b>c\)（其余两条显然成立），化简得\(n^2+n x>x^2\Longleftrightarrow\left(\dfrac{x}{n}\right)^2-\left(\dfrac{x}{n}\right)-1<0.\)

设黄金比例\(\varphi=\dfrac{1+\sqrt5}{2},\)则上述二次不等式等价于\(1\le \dfrac{x}{n}<\varphi,\)也就是\(n\le x<\varphi n.\)

那么周长为\(a+b+c=s(n^2+n x+x^2).\)记\(D(n,x)=n^2+n x+x^2,\)题目要求\(s\cdot D(n,x)\le P.\) 对固定的 \((n,x)\)，允许的 \(s\)（无平方因子）满足\(1\le s\le \left\lfloor \dfrac{P}{D(n,x)}\right\rfloor.\)因此总答案是

\[ \sum_{\substack{n\ge 1\\ n\le x<\varphi n\\D(n,x)\le P}} Q\left(\left\lfloor\frac{P}{D(n,x)}\right\rfloor\right), \]

其中 \(Q(m)\) 表示 \(\le m\) 的无平方因子数个数（包含 \(1\)）。

基于容斥原理，我们可以通过莫比乌斯函数\(\mu\)得到\(\displaystyle{Q(m)=\sum_{d=1}^{\lfloor\sqrt m\rfloor}\mu(d)\left\lfloor\frac{m}{d^2}\right\rfloor.}\)

当进行计算时，直接遍历所有 \((n,x)\) 不可行，因为 \(n\) 最大可达约 \(\sqrt{P/3}\)。

选一个分界 \(N\)，分两部分计数：

对于对每个 \(n\le N\)，枚举所有\(x\in [n,\lfloor\varphi n\rfloor]\) （严格不等式 \(x<\varphi n\) 因为 \(\varphi n\) 不会是整数，可直接用 \(\lfloor\varphi n\rfloor\)）。统计每个值\(D=n^2+nx+x^2\) 出现次数 \(\text{cnt}[D]\)，则这一部分贡献为 \[ \sum_D \text{cnt}[D]\cdot Q\left(\left\lfloor\frac{P}{D}\right\rfloor\right). \]

5.2 大 \(n>N\)：改为枚举无平方因子 \(s\)

当 \(n>N\) 时，因 \(x\ge n\)，\(\displaystyle{D(n,x)\ge 3n^2\Rightarrow s\le \left\lfloor\frac{P}{3n^2}\right\rfloor.}\)令\(\displaystyle{S_{\max}=\left\lfloor\frac{P}{3(N+1)^2}\right\rfloor},\) 所以对所有 \(n>N\) 的三角形，其 \(s\) 一定不超过 \(S_{\max}\)。

于是我们改为遍历每个无平方因子 \(s\le S_{\max}\)，令\(B=\left\lfloor\dfrac{P}{s}\right\rfloor,\)需要统计满足\(n>N, n\le x<\varphi n, D(n,x)\le B\)的整数对 \((n,x)\) 个数，记为 \(F(B)\)，并累加

\[ \sum_{\substack{1\le s\le S_{\max} \\ s\text{ squarefree}}} F\left(\left\lfloor\frac{P}{s}\right\rfloor\right). \]

对固定 \(B\)，对每个 \(n\)，不等式\(x^2+n x+(n^2-B)\le 0\)给出上界\(x\le \left\lfloor\dfrac{-n+\sqrt{4B-3n^2}}{2}\right\rfloor.\)再与 \(x<\varphi n\) 的上界取最小即可。

实现时，为了避免在 \(n\) 循环里频繁开方，用判别式单调递减的性质维护\(\Delta(n)=4B-3n^2\)及其整数平方根 \(\lfloor\sqrt{\Delta(n)}\rfloor\)，每次 \(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
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
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
#include <bits/stdc++.h>
#include "tools.h"
using namespace std;
typedef long long ll;
const ll P = 25000000000000ll;

inline ll icbrt_ll(ll x)
{
    long double d = pow((long double)x, 1.0L / 3.0L);
    ll r = (ll)d;
    auto cube = [&](ll t) -> ll
    { return (ll)t * t * t; };
    while (cube(r + 1) <= x)
        ++r;
    while (cube(r) > x)
        --r;
    return r;
}

// 线性筛 Möbius：mu[1]=1，mu[n]=0/±1
void mobius_sieve(int LIM, vector<int8_t> &mu, vector<int> &mertens, vector<int> &sqfree_prefix)
{
    mu.assign(LIM + 1, 0);
    mertens.assign(LIM + 1, 0);
    sqfree_prefix.assign(LIM + 1, 0);

    vector<int> primes;
    vector<int> lp(LIM + 1, 0);
    mu[1] = 1;

    for (int i = 2; i <= LIM; ++i)
    {
        if (lp[i] == 0)
        {
            lp[i] = i;
            primes.push_back(i);
            mu[i] = -1;
        }
        for (int p : primes)
        {
            ll v = 1LL * i * p;
            if (v > LIM)
                break;
            lp[(int)v] = p;
            if (i % p == 0)
            {
                mu[(int)v] = 0;
                break;
            }
            else
            {
                mu[(int)v] = (int8_t)(-mu[i]);
            }
        }
    }

    // Mertens 与 squarefree 前缀
    for (int i = 1; i <= LIM; ++i)
    {
        mertens[i] = mertens[i - 1] + (int)mu[i];
        sqfree_prefix[i] = sqfree_prefix[i - 1] + (mu[i] != 0);
    }
}

// ------- 去掉结构体：用全局数组 + 普通函数实现 Q(x) -------
int PREF_LIMIT = 0;
vector<int8_t> MU;
vector<int> MERTENS;
vector<int> SQFREE_PREFIX;

ll count_sqfree(ll x)
{
    if (x <= (ll)PREF_LIMIT)
        return (ll)SQFREE_PREFIX[(int)x];

    ll r = int_square(x);
    ll B = icbrt_ll(x);
    if (B > r)
        B = r;

    ll res = 0;
    for (ll i = 1; i <= B; ++i)
    {
        res += (ll)MU[(int)i] * (ll)(x / (i * i));
    }

    ll maxv = x / ((B + 1) * (B + 1));
    for (ll v = 1; v <= maxv; ++v)
    {
        ll R = int_square(x / v);
        ll L = int_square(x / (v + 1)) + 1;
        if (R <= B)
            continue;
        if (L <= B)
            L = B + 1;
        if (L <= R)
        {
            res += (ll)v * (ll)(MERTENS[(int)R] - MERTENS[(int)L - 1]);
        }
    }
    return (ll)res;
}

int main()
{
    // 分界 N：可调整
    const int N = 3000;

    // 预处理上限
    const ll nMaxAll = int_square(P / 3);                  // n 的理论最大值（x>=n 时）
    const ll Smax = P / (3ll * (ll)(N + 1) * (ll)(N + 1)); // n>N 时 s 的最大可能值

    // Squarefree 计数的查表上限
    PREF_LIMIT = 10000000;

    // Möbius / 前缀需要到 max(PREF_LIMIT, sqrt(P/3))
    const int MU_LIMIT = max<int>(PREF_LIMIT, (int)(nMaxAll + 5));

    mobius_sieve(MU_LIMIT, MU, MERTENS, SQFREE_PREFIX);

    // 预计算 x_phi[n] 与 prefix_phi[n]
    vector<int> xphi((size_t)nMaxAll + 1, 0);
    vector<ll> prefix_phi((size_t)nMaxAll + 1, 0);

    const long double phi = (1.0L + sqrtl(5.0L)) / 2.0L;

    for (ll n = 1; n <= nMaxAll; ++n)
    {
        ll x = (ll)(phi * (long double)n);
        auto ok = [&](ll xx) -> bool
        {
            // 需要 n^2 + n*xx > xx^2
            ll nn = (ll)n;
            ll xh = (ll)xx;
            return nn * nn + nn * xh > xh * xh;
        };
        while (ok(x + 1))
            ++x;
        while (!ok(x))
            --x;

        xphi[n] = (int)x;

        ll cnt = (x >= n) ? (x - n + 1) : 0;
        prefix_phi[n] = prefix_phi[n - 1] + cnt;
    }

    // -------- Part 1: n <= N，按 D 聚合 ----------
    const int D_LIMIT = 6 * N * N + 5;
    vector<int> cntD(D_LIMIT + 1, 0);

    for (int n = 1; n <= N; ++n)
    {
        int xhi = xphi[n];
        for (int x = n; x <= xhi; ++x)
        {
            int D = n * n + n * x + x * x;
            ++cntD[D];
        }
    }

    ll ans = 0;
    for (int D = 1; D <= D_LIMIT; ++D)
    {
        int c = cntD[D];
        if (!c)
            continue;
        ll q = P / (ll)D;
        ans += (ll)c * count_sqfree(q);
    }

    // -------- Part 2: n > N，枚举 squarefree s ----------
    const long double C = 3.0L + sqrtl(5.0L); // 1 + phi + phi^2 = 3 + sqrt(5)

    auto D_at_phi = [&](ll n) -> ll
    {
        ll x = (ll)xphi[n];
        return (ll)n * n + (ll)n * x + (ll)x * x;
    };

    for (ll s = 1; s <= Smax; ++s)
    {
        if (MU[(int)s] == 0)
            continue; // 非无平方因子，跳过

        ll B = P / s;
        ll nmax = int_square(B / 3ll);
        if (nmax <= (ll)N)
            continue;

        // fll_end：在 n<=fll_end 时，x 仅受 x<phi n 约束就一定满足 D<=B
        ll fll_end = (ll)sqrtl((long double)B / C);
        if (fll_end > nmax)
            fll_end = nmax;

        while (fll_end + 1 <= nmax && D_at_phi(fll_end + 1) <= (ll)B)
            ++fll_end;
        while (fll_end >= 1 && D_at_phi(fll_end) > (ll)B)
            --fll_end;

        ll add = 0;
        if (fll_end > (ll)N)
        {
            add += prefix_phi[fll_end] - prefix_phi[N];
        }

        ll start = max<ll>((ll)N + 1, fll_end + 1);
        if (start <= nmax)
        {
            ll n = start;
            ll disc = 4ll * B - 3ll * n * n; // Δ(n)=4B-3n^2
            ll y = int_square(disc);         // y=floor(sqrt(Δ))

            for (; n <= nmax; ++n)
            {
                if (n != start)
                {
                    // Δ(n) = Δ(n-1) - 3(2n-1)
                    disc -= 3ll * (2ll * n - 1ll);
                    while ((ll)y * y > (ll)disc)
                        --y;
                }

                ll xq = (ll)((y - n) / 2ll);
                ll xup = xq;
                int xph = xphi[n];
                if (xup > xph)
                    xup = xph;
                if (xup >= (ll)n)
                    add += (ll)(xup - (ll)n + 1);
            }
        }

        ans += add;
    }

    cout << ans << "\n";
    return 0;
}