Project Euler 572

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

Project Euler 572

题目

Idempotent matrices

A matrix \(M\) is called idempotent if \(M^2 = M\).

Let \(M\) be a three by three matrix :

\(M=\begin{pmatrix} a & b & c\ d & e & f\\ g & h &i\\ \end{pmatrix}\).

Let \(C(n)\) be the number of idempotent three by three matrices \(M\) with integer elements such that

\(-n \le a,b,c,d,e,f,g,h,i \le n\).

\(C(1)=164\) and \(C(2)=848\).

Find \(C(200)\).

解决方案

直接计算 \(M^2\) 并与 \(M\) 对应元素相等，得到 9 条二次方程：

对角线 3 条：

\[ \begin{aligned} a^2+bd+cg &= a,\\ e^2+bd+fh &= e,\\ i^2+cg+fh &= i. \end{aligned} \]

非对角线 6 条：

\[ \begin{aligned} ab+be+ch = b,ac+bf+ci = c,da+ed+fg = d,\\ dc+ef+fi = f,ga+hd+ig = g,gb+he+ih = h. \end{aligned} \]

令\(A=a-a^2, E=e-e^2, I=i-i^2.\)

则对角线三式等价于线性系统：\(bd+cg = A,bd+fh = E,cg+fh= I.\)

解得：\(bd = \dfrac{A+E-I}{2},cg = \dfrac{A+I-E}{2},fh = \dfrac{E+I-A}{2}.\)

因此若 \(A+E+I\) 为奇数，则 \(bd,cg,fh\) 不是整数，直接无解。

定义三个乘子：\(x=1-a-e, y=1-a-i, z=1-e-i.\)把 6 条非对角线方程移项后，可以整理成非常对称的形式：

\[ \begin{aligned} ch = bx,fg = dx,bf = cy,\\ hd = gy,cd = fz,bg = hz. \end{aligned} \]

接下来整个计数就围绕 \(x,y,z\) 是否为零分情况讨论。

\(xyz\ne 0\)（通用情形）

假设\(xyz\ne 0\)，此时从上面的 6 个乘积式可以推出一个“核心三元乘积”恒等式。

由\(bf=cy, cd=fz\)消去 \(c\) 得\(c=\dfrac{bf}{y}, d=\dfrac{fz}{c}=\dfrac{fz}{bf/y}=\dfrac{yz}{b}.\)

再由\(fg=dx\)代入 \(d=\dfrac{yz}{b}\) 得\(fg=\dfrac{yz}{b}x\iff bfg=xyz.\)

因此在通用情形必有：\(bfg=xyz.\)同理还可推出另外三个乘积恒等式：\(bd=yz,cg=xz,fh=xy.\)

于是得到三条必要一致性条件，即把 \(xy, xz, yz\) 消掉后的形式：

\[fh\cdot z^2 = bd\cdot cg,bd\cdot x^2 = fh\cdot cg,cg\cdot y^2 = bd\cdot fh.\]

这三式不成立时，通用情形无解，可以直接剪枝。

在通用情形中，对角线已经唯一确定了三个乘积：\(bd, cg, fh.\)

我们做如下构造：

  1. 枚举所有整数因子对 \((b,d)\) 满足\(b\cdot d = bd, |b|,|d|\le n.\)
  2. 枚举所有整数因子对 \((f,h)\) 满足\(f\cdot h = fh, |f|,|h|\le n.\)
  3. 由核心恒等式 \(bfg=xyz\) 确定\(g=\dfrac{xyz}{bf}.\)这一步要求 \(bf\mid xyz\)，并且 \(g\ne 0\) 且 \(|g|\le n\)。
  4. 由 \(cg=xz\) 以及 \(cg=cg\)（已知）确定\(c=\dfrac{cg}{g},\)需要 \(g\mid cg\) 且 \(|c|\le n\) 且 \(c\ne 0\)。

在这些条件都满足时，矩阵就被完全确定，并且必然满足 \(M^2=M\)。

所以通用情形的计数本质上就是：枚举 \((b,d)\) 与 \((f,h)\) 的受限因子对，再检查能否得到合法的 \(g,c\)。

为了加速，必须缓存所有 \(t\) 的因子对集合\(D_n(t)\)（特别是 \(t=0\) 要特殊处理），也就是\(D_n(t)=\{(u,v)\in \mathbb{Z}^2:uv=t,-n\le u,v\le n\}\)。

\(z=0\)（恰有一个为零）

不妨先固定在标准形：\(z=1-e-i=0\iff e+i=1, x\ne 0,y\ne 0.\)

此时两条式子变成：\(cd=fz=0, bg=hz=0.\)

而对角线解出来的乘积在这个情形会强迫：\(bd=0, cg=0.\)

因此 \(b,d,c,g\) 中会出现大量自由度，不能再用通用情形的 “\(bfg=xyz\)” 直接锁死全部变量。

这一类情况可以拆成两种情况：

\(b=c=0\)

此时 \(cg=0\) 已由对角线条件自动推出，因此不再需要额外约束。剩下只需令 \((f,h)\) 满足\(fh=\dfrac{E+I-A}{2},\)即可保证与对角线方程相容。

另一方面，由非对角项方程 \(fg=dx\) 可得\(g=\dfrac{dx}{f}.\)

因此问题转化为计数：对每一组合法的 \((f,h)\)（满足 \(|f|,|h|\le n\) 且 \(fh\) 固定），统计有多少个\(d\in[-n,n]\)使得 \(g\) 为整数，并且同时满足 \(|g|\le n\)。

这一计数可以通过一次\(\gcd(f,x)\)化简，把整除条件 + 取值范围压缩成对某个整数参数的区间计数，从而避免逐一枚举所有 \(d\)。

\(d=g=0\) 且 \(b\ne 0\)

此时 \(bd=0\) 自动满足，且 \(fh\) 仍由 \((f,h)\) 决定。由 \(bf=cy\) 得\(c=\dfrac{bf}{y},\)

所以只要统计对每个 \((f,h)\)，有多少个 \(b\ne 0\) 使得 \(y\mid bf\) 且 \(|c|\le n\)。

同样可以用一次 \(\gcd(f,y)\) 化简来快速得到可行个数。

最终，如果 \(x=0\) 或 \(y=0\)，通过对角线置换（对称性）可以转化为 \(z=0\) 的同型问题，因此只需实现 \(z=0\) 的计数函数即可。

\(x=0,y=0,z\ne 0\)（恰有两个为零）

当恰有两个乘子为零时（例如 \(x=0,y=0,z\ne 0\)），它只会出现在非常特殊的对角线配置上。

以迹为 \(1\) 的情况为例：\(x=0\iff a+e=1, y=0\iff a+i=1\Rightarrow e=i,a+2e=1\Rightarrow(a,e,i)=(1,0,0)\)，其余是该对角线的排列。

在这种对角线之下，对角线三式给出：\(bd=cg=fh=0.\)

并且由 \(M^2=M\) 的非对角线方程可以化简到：\(f=dc,h=gb\)并且约束只剩下\(bd=0, cg=0, |dc|\le n, |gb|\le n.\)

也就是说，只要选定 \(b,c,d,g\)，就自动确定 \(f,h\)，其余方程全部成立。

设\(\displaystyle{H=\sum_{k=1}^n\left\lfloor \frac{n}{k}\right\rfloor.}\)这正是正整数对 \((u,v)\) 满足 \(uv\le n\) 的个数。

在在此情形中，解集可以分成三类互不相交的情况：

  1. \(dc\ne 0\)（即 \(d\ne 0,c\ne 0\)），此时必须有 \(b=0,g=0\)（否则 \(bd=0\) 与 \(cg=0\) 会矛盾），并且需要\(|dc|\le n.\)正数对数量为 \(H\)，带符号共有 \(4H\) 种。
  2. \(gb\ne 0\)（即 \(g\ne 0,b\ne 0\)），同理得到 \(d=0,c=0\)，并且需要 \(|gb|\le n\)，带符号也是 \(4H\)。
  3. \(dc=0\) 且 \(gb=0\)，这等价于\((d=0\lor c=0)\land (g=0\lor b=0),\)在同时满足 \(|b|,|c|,|d|,|g|\le n\) 下可数出总数为\(8n(n+1)+1.\)

因此这种情况下的总数为：

\[T(n)=8H+8n(n+1)+1.\]

可见，幂等矩阵的特征多项式满足\(M(M-I)=0,\)在特征为 \(0\) 的域上其最小多项式无重根，因此所有特征值只可能为 \(0\) 或 \(1\)。于是\(\operatorname{tr}(M)\in\{0,1,2,3\}.\)

其中：

  • \(\operatorname{tr}(M)=0\) 只对应零矩阵 \(M=0\)；
  • \(\operatorname{tr}(M)=3\) 只对应单位矩阵 \(M=I\)；

所以最终只需要枚举：\(\operatorname{tr}(M)=1,2\) 并在最后补上这两个特殊解，总和为：

\[ C(n)=2+\sum_{\operatorname{tr}=1,2}\sum_{a,e}\#(a,e,i=\operatorname{tr}-a-e). \]

实现中外层就是遍历 \(a,e\)，确定\(i=\operatorname{tr}-a-e\)再根据 \(x,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
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
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
#include <bits/stdc++.h>
using namespace std;

// ------------------------------------------------------------
// Project Euler 572 - Idempotent 3x3 integer matrices
// Count matrices M with entries in [-LIM, LIM] such that M^2 = M.
// ------------------------------------------------------------

// --------- Global parameters ----------
typedef long long ll;
int LIM = 200;
ll LIM2; // LIM^2

// Cache for factor pairs of t (u,v) with u*v=t, u,v in [-LIM,LIM]
unordered_map<ll, vector<pair<int, int>>> div_cache;

// ---------- helpers ----------
bool ok(ll x) { return -LIM <= x && x <= LIM; }

// All integer pairs (u,v) with u*v=t, u,v in [-LIM, LIM]
const vector<pair<int, int>> &factorPairs(ll t)
{
    auto it = div_cache.find(t);
    if (it != div_cache.end())
        return it->second;

    vector<pair<int, int>> res;

    if (t == 0)
    {
        // (u,0) for all u, and (0,v) for all v != 0, include (0,0) once
        res.reserve(2 * (2 * LIM + 1) - 1);
        for (int u = -LIM; u <= LIM; ++u)
            res.push_back({u, 0});
        for (int v = -LIM; v <= LIM; ++v)
            if (v != 0)
                res.push_back({0, v});
        div_cache.emplace(t, std::move(res));
        return div_cache.find(t)->second;
    }

    ll at = abs(t);
    // If |t| > LIM^2, no (u,v) within bounds can satisfy u*v=t
    if (at > 1LL * LIM * LIM)
    {
        div_cache.emplace(t, std::move(res));
        return div_cache.find(t)->second;
    }

    ll rt = (ll)sqrt((long double)at);
    for (ll d = 1; d <= rt; ++d)
    {
        if (at % d)
            continue;
        ll q = at / d;

        // generate signed pairs (u,v) such that u*v=t
        // and also include swapped (v,u)
        auto addPair = [&](ll u, ll v)
        {
            if (ok(u) && ok(v))
            {
                res.push_back({(int)u, (int)v});
            }
        };

        if (t > 0)
        {
            addPair(d, q);
            addPair(-d, -q);
            if (d != q)
            {
                addPair(q, d);
                addPair(-q, -d);
            }
        }
        else
        {
            addPair(d, -q);
            addPair(-d, q);
            if (d != q)
            {
                addPair(q, -d);
                addPair(-q, d);
            }
        }
    }

    div_cache.emplace(t, std::move(res));
    return div_cache.find(t)->second;
}

// Precompute TwoZeroCount for given LIM:
// counts solutions for the (canonical) diagonal pattern with exactly two zero multipliers,
// which always ends up with bd=cg=fh=0 and two of x,y,z are 0.
// The final count only depends on LIM, not on the specific diagonal permutation.

ll TWO_ZERO_COUNT = 0;

// -------- Case 0: x,y,z all nonzero --------
ll countCommon(ll bd, ll cg, ll fh,
               ll x, ll y, ll z)
{
    // Consistency conditions (necessary and sufficient)
    // If any fails -> no solutions.
    // Use 128-bit for safety (though values are small for LIM=200).
    ll BD = bd, CG = cg, FH = fh, X = x, Y = y, Z = z;

    if (FH * Z * Z != BD * CG || BD * X * X != FH * CG || CG * Y * Y != BD * FH)
        return 0;

    const auto &BDPairs = factorPairs(bd);
    const auto &FHPairs = factorPairs(fh);

    ll mul = x * y * z;
    ll ans = 0;

    for (const auto &p1 : BDPairs)
    {
        ll b = p1.first;
        if (b == 0)
            continue; // bd != 0 in this case
        for (const auto &p2 : FHPairs)
        {
            ll f = p2.first;
            if (f == 0)
                continue; // fh != 0 in this case
            ll denom = b * f;
            if (mul % denom != 0)
                continue;
            ll g = mul / denom;
            if (!ok(g) || g == 0 || cg % g != 0)
                continue;
            ll c = cg / g;
            if (!ok(c) || c == 0)
                continue;

            // Unique matrix determined by (a,e,i,b,c,d,f,g,h)
            ans++;
        }
    }
    return ans;
}

// -------- Case 1: exactly one of x,y,z is zero --------
// We normalize by permuting diagonal so that z=0 (i.e. e+i=1).
ll countOneZero_Z0(int a, int e, int i)
{
    // Here z=0, and x,y !=0
    ll x = 1LL - a - e, y = 1LL - a - i, z = 1LL - e - i;
    if (z != 0 || x == 0 || y == 0)
        return 0;

    ll A = 1LL * a - 1LL * a * a, E = 1LL * e - 1LL * e * e, I = 1LL * i - 1LL * i * i;

    if (((A + E + I) & 1) != 0)
        return 0;
    ll bd = (A + E - I) / 2, cg = (A + I - E) / 2, fh = (E + I - A) / 2;

    // In this situation, we must have bd=cg=0, otherwise contradictions arise.
    if (bd != 0 || cg != 0)
        return 0;

    // fh must be representable as product of two bounded ints
    const auto &FHPairs = factorPairs(fh);
    if (FHPairs.empty())
        return 0;

    ll ans = 0;

    // Family A: b=0, c=0. Variables (f,h) factor of fh, then choose d, g is determined.
    // We can count d without looping by reducing divisibility.
    for (const auto &pr : FHPairs)
    {
        ll f = pr.first;
        ll h = pr.second;
        if (f == 0 || h == 0)
            continue; // fh != 0 in this case

        // Need d in [-LIM,LIM] such that g = d*x/f is integer and within [-LIM,LIM].
        // Let g0 = gcd(|f|,|x|), f = g0*f1, x = g0*x1 (gcd(f1,x1)=1).
        // Then d must be multiple of |f1|, write d = f1*t, and g = t*x1.
        ll af = abs(f);
        ll ax = abs(x);
        ll g0 = __gcd(af, ax);
        ll f1 = af / g0;
        ll x1 = ax / g0;

        ll T = min((ll)(LIM / f1), (ll)(LIM / x1));
        // t in [-T..T]  => count = 2T+1
        ans += (2 * T + 1);

        // Family B: d=0, g=0, but b != 0.
        // Need c = b*f / y integer and |c|<=LIM.
        // Reduce: let g1=gcd(|f|,|y|), f=g1*f2, y=g1*y2 (gcd(f2,y2)=1)
        // then b must be multiple of |y2|, b=y2*t (t!=0),
        // and c=t*f2, so |t| <= min(LIM/|y2|, LIM/|f2|).
        ll ay = abs(y);
        ll g1 = __gcd(af, ay);
        ll f2 = af / g1;
        ll y2 = ay / g1;

        ll T2 = min((ll)(LIM / y2), (ll)(LIM / f2));
        // t != 0 => count = 2*T2
        ans += (2 * T2);
    }

    return ans;
}

// Main per-diagonal counter
ll countDiagonal(int a, int e, int i)
{
    ll x = 1LL - a - e, y = 1LL - a - i, z = 1LL - e - i;

    int zero = int(x == 0) + int(y == 0) + int(z == 0);

    ll A = 1LL * a - 1LL * a * a, E = 1LL * e - 1LL * e * e, I = 1LL * i - 1LL * i * i;

    if (((A + E + I) & 1) != 0)
        return 0;
    ll bd = (A + E - I) / 2, cg = (A + I - E) / 2, fh = (E + I - A) / 2;

    if (zero == 0)
    {
        // bd,cg,fh must be representable within bounds
        if (abs(bd) > LIM2 || abs(cg) > LIM2 || abs(fh) > LIM2)
            return 0;
        return countCommon(bd, cg, fh, x, y, z);
    }
    else if (zero == 1)
    {
        // normalize to z=0 by diagonal permutation (counts are invariant under permutation conjugation)
        if (z == 0)
        {
            return countOneZero_Z0(a, e, i);
        }
        else if (x == 0)
        {
            // (a,e,i) -> (i,a,e) makes z' = a+e? actually e'+i' = a+e = 1
            return countOneZero_Z0(i, a, e);
        }
        else
        {
            // y==0: (a,e,i) -> (e,a,i) makes e'+i' = a+i = 1
            return countOneZero_Z0(e, a, i);
        }
    }
    else if (zero == 2)
    {
        // Only possible (for trace 1 or 2) on the special 0/1 diagonals, and then bd=cg=fh=0.
        if (bd != 0 || cg != 0 || fh != 0)
            return 0;
        return TWO_ZERO_COUNT;
    }
    else
    {
        // x=y=z=0 impossible for integer diagonal with trace 1 or 2
        return 0;
    }
}

int main()
{
    LIM2 = 1LL * LIM * LIM;
    // Let H = sum_{k=1..LIM} floor(LIM/k)
    ll H = 0;
    for (int k = 1; k <= LIM; k++)
        H += LIM / k;
    // Derived closed-form:
    // TwoZeroCount = 8*H + 8*LIM*(LIM+1) + 1
    // (works for x=±1 canonical, but independent of sign)
    TWO_ZERO_COUNT = 8LL * H + 8LL * LIM * (LIM + 1) + 1;

    ll total = 1;
    if (LIM >= 1)
        ++total;
    // Only trace 1 or 2 are possible for non-zero idempotent 3x3 matrices
    for (int tr = 1; tr <= 2; ++tr)
    {
        for (int a = -LIM; a <= LIM; ++a)
        {
            for (int e = -LIM; e <= LIM; ++e)
            {
                int i = tr - a - e;
                if (i < -LIM || i > LIM)
                    continue;
                total += countDiagonal(a, e, i);
            }
        }
    }
    cout << total << "\n";
}