Project Euler 324

Project Euler 324

题目

Building a tower

Let \(f(n)\) represent the number of ways one can fill a \(3\times3\times n\) tower with blocks of \(2\times1\times1\).

You’re allowed to rotate the blocks in any way you like; however, rotations, reflections etc of the tower itself are counted as distinct.

For example (with \(q = 100000007\)) :

\(\begin{aligned} &f(2) = 229,\\ &f(4) = 117805,\\ &f(10) \bmod q = 96149360,\\ &f(10^3) \bmod q = 24806056,\\ &f(10^6) \bmod q = 30808124.\\ \end{aligned}\)

Find \(f(10^{10000}) \bmod 100000007\).

解决方案

\(N=10^{10000}\)

\(f(2)\)\(f(4)\)放入OEIS中查找,发现结果为A028452

找到FORMULA一栏,发现如下信息:

1
a(n) = 679a(n-1) -76177a(n-2) +3519127a(n-3) -85911555a(n-4) +1235863045a(n-5) -11123194131a(n-6) +65256474997a(n-7) -257866595482a(n-8) +705239311926a(n-9) -1363115167354a(n-10) +1888426032982a(n-11) -1888426032982a(n-12) +1363115167354a(n-13) -705239311926a(n-14) +257866595482a(n-15) -65256474997a(n-16) +11123194131a(n-17) -1235863045a(n-18) +85911555a(n-19) -3519127a(n-20) +76177a(n-21) -679a(n-22) +a(n-23).

其中,这里的\(a(n)=f(2n)\)。因为很明显知道,当\(m\)为奇数时,\(f(m)=0\)

因此此时就变成了求\(a\left(\dfrac{N}{2}\right)\)的值。

注意到\(a\)的递推式有\(M=23\)个项,因此使用矩阵快速幂加速时,所使用的矩阵大小\(M\)。因此,考虑使用第258题中的哈密顿-凯莱定理,将计算递推式的时间复杂度进一步降至\(O(M^2 \log N)\)

代码

此处使用GMP库对\(N\)进行处理。

本份代码由第258题的代码修改而来。

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

const int P=10000;
mpz_class N(string("1")+string(P,'0'));
ll mod = 100000007;

const int M=23;
// a代表递推式的前M个值
ll a[M];
string string_a[M]= {
"1",
"229",
"117805",
"64647289",
"35669566217",
"19690797527709",
"10870506600976757",
"6001202979497804657",
"3313042830624031354513",
"1829008840116358153050197",
"1009728374600381843221483965",
"557433823481589253332775648233",
"307738670509229621147710358375321",
"169891178715542584369273129260748045",
"93790658670253542024618689133882565125",
"51778366130057389441239986148841747669217",
"28584927722109981792301610403923348017948449",
"15780685138381102545287108197623881881376915397",
"8711934690116480171969789787256390490181022415693",
"4809538076408327645969201260680362259835079086427481",
"2655168723276120197512956906659822833388644760430125609",
"1465820799640802552047402979496052449322258430218930512765",
"809225642733724788155919446555896648357335949987871250500245"
};
// p代表Hamilton-Cayley而来的矩阵多项式的系数。
ll p[M]= {
1,
-679,
76177,
-3519127,
85911555,
-1235863045,
11123194131,
-65256474997,
257866595482,
-705239311926,
1363115167354,
-1888426032982,
1888426032982,
-1363115167354,
705239311926,
-257866595482,
65256474997,
-11123194131,
1235863045,
-85911555,
3519127,
-76177,
679
};
// b代表A^(2^i)时所有的A^(2^i)=sum from j=0 to M-1 A^j*b_j时系数b_j的值。
ll b[M];
// v代表解决方案所指的v。
ll v[M];
void mul_ploy(ll a[M],ll b[M],ll p[M]) {
//做的是乘法,对应的就是多项式系数的卷积。
ll c[M << 1];
memset(c, 0, sizeof(c));
for (int i = 0; i < M; i++)
for (int j = 0; j < M; j++)
c[i + j] = (c[i + j] + a[i] * b[j]) % mod;
//利用递推式F(A)的值,将>=M的所有临时系数转化成<M的系数。
for (int i = M * 2 - 1; i >= M; i--)
for (int j = 0; j < M; j++)
c[i - M + j] = (c[i - M + j] + c[i] * p[j]) % mod;
for (int i = 0; i < M; i++)
a[i] = c[i];
}
int main() {
ll ans = 0;
if ((N & 1) == 0) {
N >>= 1;
for (int i = 0; i < M; i++) {
for (char ch:string_a[i])
a[i] = (a[i] * 10 + ch - '0') % mod;
}
for (int i = 0; i < M; i++)
p[i] = (p[i] % mod + mod) % mod;
b[1] = 1;
v[0] = 1;
for (; N; N >>= 1) {
if ((N & 1) == 1) mul_ploy(v, b, p);
mul_ploy(b, b, p);
}
for (int i = 0; i < M; i++)
ans = (ans + v[i] * a[i]) % mod;
}
printf("%lld\n", ans);
}