Project Euler 324

发表于 更新于 分类于 阅读次数： 本文字数： 792 阅读时长 ≈ 1 分钟

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);
}