Project Euler 397

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

Project Euler 397

题目

Triangle on parabola

On the parabola \(y = x^2/k\), three points \(A(a, a^2/k), B(b, b^2/k)\) and \(C(c, c^2/k)\) are chosen.

Let \(F(K, X)\) be the number of the integer quadruplets \((k, a, b, c)\) such that at least one angle of the triangle \(ABC\) is \(45\)-degree, with \(1 \le k \le K\) and \(-X \le a < b < c \le X\).

For example, \(F(1, 10) = 41\) and \(F(10, 100) = 12492\).

Find \(F(10^6, 10^9)\).

解决方案

两点\(P_u=\left(u,\dfrac{u^2}{k}\right), P_v=\left(v,\dfrac{v^2}{k}\right)\) 的差向量为\(\overrightarrow{P_uP_v}=\left(v-u,\dfrac{v^2-u^2}{k}\right)=\left(v-u,\dfrac{(v-u)(v+u)}{k}\right).\)

因此距离平方为\(|P_v-P_u|^2=(v-u)^2\left(1+\dfrac{(v+u)^2}{k^2}\right).\)令\(L_{uv}^2=(v-u)^2\left(k^2+(v+u)^2\right),\)则真实距离平方满足 \(|P_v-P_u|^2=L_{uv}^2/k^2\)。由于角度只与比值有关，我们可以用 \(L_{uv}^2\) 替代计算。

若某顶点角为 \(45^\circ\)，则该角的余弦平方满足\(\cos^2\theta=\dfrac12.\)以顶点为例（取某一顶点、对应两条边的夹角），将余弦公式写为\(\cos^2\theta=\dfrac{(u\cdot v)^2}{|u|^2|v|^2}\)，并代入上面的边长平方表达式进行化简，可以得到一个可因式分解的整数多项式条件。

对非退化三角形（\(a<b<c\)）而言，关键约束等价于要求下面两类二次因子之一为 \(0\)：\(ab+ac-ak+bc+bk+c^2+k^2=0\)或\(ab+ac+ak+bc-bk+c^2+k^2=0.\)

这两类方程在整体变换 \((a,b,c)\mapsto(-a,-b,-c)\) 下互相对应，因此只需处理其中一类并在最终计数中自然包含对称贡献即可。下面固定讨论\(ab+ac-ak+bc+bk+c^2+k^2=0.\)

把它视为关于 \(c\) 的二次方程：\(c^2+(a+b)c+(ab-ak+bk+k^2)=0.\)判别式为\(\Delta=(a+b)^2-4(ab-ak+bk+k^2).\)展开整理后可写成\(\Delta=(a-b+2k)^2-8k^2.\)要使 \(c\) 为整数，必须有\((a-b+2k)^2-8k^2=d^2\) 对某个整数 \(d\) 成立，即\(d^2+2(2k)^2=(a-b+2k)^2.\)

设\(u=d, v=2k, w=a-b+2k,\)就得到经典形式\(u^2+2v^2=w^2.\)

该方程的本原解可用如下参数化给出：取互素正整数 \((m,n)=1\) 且 \(m\) 为奇数，则

\[w=m^2+2n^2, u=m^2-2n^2, v=2mn\]

给出一组本原解。

再乘以任意缩放因子 \(l\ge 1\) 可得到全部解：

\[ w=l(m^2+2n^2), u=l(m^2-2n^2), v=2lmn. \]

由 \(v=2k\) 得 \(k=lmn.\)

由定义\(w=a-b+2k,\)因此\(b-a=2k-w.\)在不同符号选择、以及顶点轮换（哪一个角取 \(45^\circ\)）下，会导出若干等价的线性族。将它们整理后，可以归并成两大族（同时，每族再含两个 \(\gamma\) 分支）。

第 1 族

第 1 族取\(w=-l(m^2+2n^2)\)。我们令\(k_0=mn\)，也就是本原解对应的最小 \(k\)。可见，\(b-a=2k-w=l((m+n)^2+n^2)\)，同样取\(b-a\)最小的本原解，那么有\(\beta=(m+n)^2+n^2,\)

另一方面，二次方程求根形式可写作\(2c=-(a+b)\pm u,\)而参数化同时给出\(u=\pm l(m^2-2n^2).\) 将 \(b=a+l\beta\) 代入 \(a+b=2a+l\beta\)，得到\(c=-a-\dfrac{l\beta}{2}\pm\dfrac{l(m^2-2n^2)}{2}.\) 为了避免分数，把 \[ \gamma=\dfrac{-\beta\pm(m^2-2n^2)}{2} \] 记成整数常量（这里 \(m\) 为奇数可保证分子为偶数），两种符号分别给出：

  • 若取 \(+(m^2-2n^2)\)，则\(\gamma=\dfrac{-\beta+(m^2-2n^2)}{2}=\dfrac{-(m^2+2mn+2n^2)+(m^2-2n^2)}{2}=-n(m+2n),\)
  • 若取 \(-(m^2-2n^2)\)，则\(\gamma=\dfrac{-\beta-(m^2-2n^2)}{2}=\dfrac{-(m^2+2mn+2n^2)-m^2+2n^2}{2}=-m(m+n).\)

因此对每个\(\gamma\in\{-m(m+n),-n(m+2n)\}.\)都要考虑。

对每个缩放 \(l=1,2,\dots,\left\lfloor\dfrac{K}{k_0}\right\rfloor\)，令真实 \(k=lk_0\)，并构造 \(b=-a+l\gamma, c=a+l\beta.\)此时只需统计满足\(-X\le a<b<c\le X\)的整数 \(a\) 数量。

由不等式：

  • \(a<b\) 等价于\(2a<l\gamma.\)
  • \(b<c\) 等价于\(-a+l\gamma<a+l\beta\iff 2a>l(\gamma-\beta).\)
  • \(c\le X\) 等价于\(a\le X-l\beta.\)
  • \(b\ge -X\) 等价于\(-a+l\gamma\ge -X\iff a\le X+l\gamma.\)

最终得到 \(a\) 的整数区间：\(\left\lfloor\dfrac{l(\gamma-\beta)}{2}\right\rfloor+1\le a\le\left\lfloor\dfrac{l\gamma-1}{2}\right\rfloor,\)并与盒约束合并：\(a\ge -X,a\le X-l\beta,a\le X+l\gamma.\)最终得到

\[\max\left\{-X,\left\lfloor\dfrac{l(\gamma-\beta)}{2}\right\rfloor+1\right\}\le a\le\min\left\{\left\lfloor\dfrac{l\gamma-1}{2}\right\rfloor,X-l\beta,X+l\gamma\right\}\]

第 2 族

第 2 族取\(w=+l(m^2+2n^2).\)则\(b-a=2k-w=2lmn-l(m^2+2n^2)=l(2mn-m^2-2n^2)=-l((m-n)^2+n^2),\) 于是自然令\(\beta= -((m-n)^2+n^2),\)从而\(b-a=l\beta.\)

和第 1族 一样，我们得到\(c=-a-\dfrac{l\beta}{2}\pm\dfrac{l(m^2-2n^2)}{2}\)，并记\(\gamma=\dfrac{-\beta\pm(m^2-2n^2)}{2},\)

这里 \(m\) 为奇数可保证分子为偶数，从而 \(\gamma\) 为整数。两种符号分别给出：

  • 若取 \(+(m^2-2n^2)\)，则\(\gamma=\dfrac{-\beta+(m^2-2n^2)}{2}=\dfrac{(m^2-2mn+2n^2)+(m^2-2n^2)}{2}=m(m-n),\)
  • 若取 \(-(m^2-2n^2)\)，则\(\gamma=\dfrac{-\beta-(m^2-2n^2)}{2}=\dfrac{(m^2-2mn+2n^2)-(m^2-2n^2)}{2}=n(2n-m).\)

因此对每个\(\gamma\in\{m(m-n),n(2n-m)\}\)都要考虑。

与第 1 族相同，为了让计数关于 \(a\) 更对称（并方便一次性乘 \(2\) 处理对称解），使用同样的线性构造形式：\(b=-a+l\gamma, c=a+l\beta.\)

对每个缩放\(l=1,2,\dots,\left\lfloor\dfrac{K}{k_0}\right\rfloor\)我们只需统计满足\(-X\le a<b<c\le X\)的整数 \(a\) 数量。

先把三点都落在 \([-X,X]\) 的限制写出来：

  • \(a\in[-X,X]\) 给出\(-X\le a\le X.\)
  • \(b\in[-X,X]\) 且 \(b=-a+l\gamma\) 给出\(-X\le -a+l\gamma\le X\iff -X+l\gamma\le a\le X+l\gamma.\)
  • \(c\in[-X,X]\) 且 \(c=a+l\beta\) 给出\(-X\le a+l\beta\le X\iff -X-l\beta\le a\le X-l\beta.\)

因此把三组限制取交集，得到第 2 族的 基本可行区间：

\[ \max\{-X-l\beta,-X+l\gamma\} \le a\le \min\{X,X+l\gamma,X-l\beta\}. \]

由于本族中 \(\beta<0\)，于是 \(X-l\beta\ge X\)，通常右端的 \(X-l\beta\) 不会成为真正限制，因此常写成更紧凑的形式：

\[ \max\{-X-l\beta,-X+l\gamma\} \le a\le \min\{X,X+l\gamma\}. \]

与第 1 族不同，第 2 族由 \(\cos^2\theta=\dfrac12\) 推出的解 可能对应 \(\theta=45^\circ\) 也可能对应 \(\theta=135^\circ\)，两者的 \(\cos^2\) 相同，但 余弦符号相反。因此必须额外要求该角满足\(\cos\theta\ge 0.\)

对这一族的线性构造 \(b=-a+l\gamma,c=a+l\beta\)，可以把余弦符号化简为一个只与 \(a\) 相关的不等式（常数因子在本族固定为正），等价于：\((2a-l\gamma)(2a+l(\beta-\gamma))\ge 0.\)

注意\(2a+l(\beta-\gamma)=2a-l(\gamma-\beta),\)于是该条件等价于\((2a-l\gamma)(2a-l(\gamma-\beta))\ge 0,\)也就是说 \(2a\) 必须落在两个根的 同侧，因此不允许处在它们之间。换句话说，必须剔除中间带：

\[ \frac{l\gamma}{2}\le a\le \frac{l(\gamma-\beta)}{2}. \]

把它写成整数端点，就是剔除区间\(\left\lceil\dfrac{l\gamma}{2}\right\rceil\le a\le\left\lfloor\dfrac{l(\gamma-\beta)}{2}\right\rfloor.\)

再与上一步的基本可行区间取交，得到最终需要 从总计数中减去 的那一段： \[ \max\left\{ -X-l\beta,-X+l\gamma, \left\lceil\frac{l\gamma}{2}\right\rceil \right\} \le a\le \min\left\{ X,X+l\gamma, \left\lfloor\frac{l(\gamma-\beta)}{2}\right\rfloor \right\}. \]

综上，第 2 族对固定 \((m,n)\)、固定 \(\gamma\)、固定 \(l\) 时：

  1. 先把基本区间内所有 \(a\) 计入：\(\max\{-X-l\beta,-X+l\gamma\}\le a\le\min\{X,X+l\gamma\};\)
  2. 再剔除中间不合法区间：\(\left\lceil\dfrac{l\gamma}{2}\right\rceil\le a\le\left\lfloor\dfrac{l(\gamma-\beta)}{2}\right\rfloor.\)

因此第 2 族对 \(a\) 的有效取值集合实际上是一个“两段式”的并集：

\[ [\max\{-X-l\beta,-X+l\gamma\},\min\{X,X+l\gamma\}]\setminus\left[\left\lceil\frac{l\gamma}{2}\right\rceil,\left\lfloor\frac{l(\gamma-\beta)}{2}\right\rfloor\right], \]

并在实现中用总数减去剔除段的方式完成计数。

若同一个三角形出现两个 \(45^\circ\) 角，则它是一个 \(45^\circ-45^\circ-!90^\circ\) 的等腰直角三角形。 这类三角形会被不同的生成方式重复计入：例如我们用“第 1 族/第 2 族”的构造保证某一个顶点角为 \(45^\circ\)，但如果另一个顶点角恰好也为 \(45^\circ\)，那么同一个三角形也可以被“把另一个顶点当作 \(45^\circ\) 顶点”的方式再次生成，因此需要对这种重复进行扣除。

回想一下，假设点 \(P_u(u, u^2/k)\)、\(P_v(v, v^2/k), P_w(w,w^2/k)\) 在抛物线 \(y=x^2/k\) 上，向量\(\overrightarrow{P_uP_v}=\left(v-u,\dfrac{v^2-u^2}{k}\right)=(v-u)\left(1,\dfrac{v+u}{k}\right).\)

因此两条边 \(\overrightarrow{P_uP_v}\cdot \overrightarrow{P_uP_w}\) 的点积为\(\overrightarrow{PQ}\cdot\overrightarrow{PR}=(v-u)(w-u)\left(1+\dfrac{(v+u)(w+u)}{k^2}\right).\)在非退化情形下 \(u\neq v, u\neq w\)，于是\(\angle QPR=90^\circ\iff1+\dfrac{(v+u)(w+u)}{k^2}=0\iff(v+u)(w+u)=-k^2.\)

这就是我们用来“锁死”两角 \(45^\circ\) 的核心： 若出现两个 \(45^\circ\)，则第三角为 \(90^\circ\)，从而必须满足上式。

对固定的 \((m,n,l,\beta,\gamma)\)，我们在计数时使用的线性构造为\(k=lk_0, b=-a+l\gamma, c=a+l\beta.\) 于是有三个“和”的简化表达：\(a+b=l\gamma, a+c=2a+l\beta, b+c=l(\beta+\gamma).\)

如果同一个三角形出现两个 \(45^\circ\)，则第三角为 \(90^\circ\)。由于三角形只有三个顶点，所以“直角顶点”只可能落在两个位置（对应两种重复来源），分别给出两个候选的 \(a\)。

直角在顶点 \(A\)

即 \(\angle BAC=90^\circ\)。由上一小节的直角条件\((a+b)(a+c)=-k^2.\)代入\(a+b=l\gamma, a+c=2a+l\beta, k=lk_0\)，得到\(l\gamma(2a+l\beta)=-(lk_0)^2.\)

两边除以 \(l\) 并整理：\(2a\gamma=-l(k_0^2+\beta\gamma),\)，从而有\(a=-\dfrac{l(k_0^2+\beta\gamma)}{2\gamma}.\)

直角在顶点 \(C\)

即 \(\angle BCA=90^\circ\)。直角条件为\((c+a)(c+b)=-k^2.\)注意这里\(c+a=a+c=2a+l\beta, c+b=b+c=l(\beta+\gamma),\)代入同样得到\((2a+l\beta)\cdot l(\beta+\gamma)=-(lk_0)^2.\)

两边除以 \(l\) 并整理：\(2a(\beta+\gamma)=-l\bigl(k_0^2+\beta(\beta+\gamma)\bigr)=-l(k_0^2+\beta^2+\beta\gamma),\)从而\(a=-\dfrac{l(k_0^2+\beta^2+\beta\gamma)}{2(\beta+\gamma)}.\)

综上，在固定 \((m,n,l,\beta,\gamma)\) 的前提下，上面两种情况分别给出一个一次方程，因此每种情况至多产生一个候选 \(a\)。也就是说：

  • 要么没有整数解；
  • 要么恰好出现一个整数 \(a\) 落在我们已计数区间内；
  • 若出现，则这一个点对应的三角形会被重复计数一次，需要减去 \(1\)。

因此最终只需检查以下两种候选：

\[ a=-\frac{l(k_0^2+\beta\gamma)}{2\gamma},a=-\frac{l(k_0^2+\beta^2+\beta\gamma)}{2(\beta+\gamma)}. \]

最终我们需要遍历所有互素正整数对 \((m,n)\)。此时我们使用 Stern–Brocot 树 可无须显式 \(\gcd\)进行判定：

从节点 \(\dfrac{n}{m}\) 生成两个子节点：\(\dfrac{n}{n+m},\dfrac{m}{n+m}\).只要满足\((m+n)m\le K, (m+n)n\le K\)就继续扩展。

同时对每个节点同时计算 \((m,n)\) 与 \((n,m)\) 两个方向的贡献。最后由于全局对称会产生双计数，将总和除以 \(2\) 即得最终 \(F(K,X)\)。

代码

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
#include <bits/stdc++.h>
using namespace std;

typedef long long ll;
typedef __int128_t i128;

ll K = 1000000LL, X = 1000000000LL;

ll floordiv(ll a, ll b)
{
    ll q = a / b, r = a % b;
    if (r && ((r > 0) != (b > 0)))
        --q;
    return q;
}

ll count_mn(ll m, ll n, ll K, ll X)
{
    if ((m & 1) == 0)
        return 0;

    ll k0 = m * n;
    ll maxl = K / k0;
    ll kk = k0 * k0;
    ll res = 0;

    auto sub_iso = [&](ll ba, ll ca, ll l, ll min_a, ll max_a)
    {
        if (ca)
        {
            i128 num = (i128)l * ((i128)kk + (i128)ba * (i128)ca);
            i128 den = (i128)2 * ca;
            if (num % den == 0)
            {
                ll a0 = (ll)(-(num / den));
                if (min_a <= a0 && a0 <= max_a)
                    --res;
            }
        }
        ll bc = ba + ca;
        if (bc)
        {
            i128 num = (i128)l * ((i128)kk + (i128)ba * ba + (i128)ba * ca);
            i128 den = (i128)2 * bc;
            if (num % den == 0)
            {
                ll a0 = (ll)(-(num / den));
                if (min_a <= a0 && a0 <= max_a)
                    --res;
            }
        }
    };

    auto add_seg = [&](ll min_a, ll max_a, ll ba, ll ca, ll l)
    {
        if (max_a < min_a)
            return;
        res += 2 * (max_a - min_a + 1);
        sub_iso(ba, ca, l, min_a, max_a);
    };

    {
        ll s = m + n;
        ll ba = s * s + n * n;
        array<ll, 2> cas = {-m * (m + n), -n * (m + 2 * n)};

        for (ll ca : cas)
        {
            for (ll l = 1; l <= maxl; ++l)
            {
                ll lca = l * ca;
                ll lba = l * ba;

                ll min_a = floordiv(l * (ca - ba), 2) + 1;
                min_a = max(min_a, -X);

                ll max_a = floordiv(lca, 2) - ((lca & 1) == 0);
                max_a = min(max_a, X - lba);
                max_a = min(max_a, X + lca);

                add_seg(min_a, max_a, ba, ca, l);
            }
        }
    }

    {
        ll d = m - n;
        ll ba = -(d * d + n * n);
        array<ll, 2> cas = {m * (m - n), n * (2 * n - m)};

        for (ll ca : cas)
        {
            for (ll l = 1; l <= maxl; ++l)
            {
                ll lca = l * ca;
                ll lba = l * ba;

                ll min_a = max(-X - lba, -X + lca);
                ll max_a = min(X, X + lca);

                if (max_a < min_a)
                    continue;

                add_seg(min_a, max_a, ba, ca, l);

                ll max_a2 = min(max_a, floordiv(l * (ca - ba), 2));
                ll min_a2 = max(min_a, floordiv(lca, 2) + (lca & 1));
                if (max_a2 >= min_a2)
                    res -= 2 * (max_a2 - min_a2 + 1);
            }
        }
    }

    return res;
}

ll solve(ll K, ll X)
{
    ll res = count_mn(1, 1, K, X);

    vector<pair<ll, ll>> st;
    st.reserve(1 << 20);
    st.push_back({2, 1});

    while (!st.empty())
    {
        auto [m, n] = st.back();
        st.pop_back();

        ll s = m + n;
        if (s * m <= K)
            st.push_back({s, m});
        if (s * n <= K)
            st.push_back({s, n});

        res += count_mn(m, n, K, X);
        res += count_mn(n, m, K, X);
    }

    return res / 2;
}

int main()
{
    cout << solve(K, X) << "\n";
    return 0;
}